PROGRAM f04ycfe

!      F04YCF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nout = 6
!      .. Executable Statements ..
       WRITE (nout,*) 'F04YCF Example Program Results'

       CALL ex1

       CALL ex2

    CONTAINS
       SUBROUTINE ex1

!         .. Use Statements ..
          USE nag_library, ONLY : dasum, dtrsv, f03aff, f04ycf, nag_wp
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          REAL (KIND=nag_wp), PARAMETER   :: zero = 0.0E+0_nag_wp
          INTEGER, PARAMETER              :: nin = 5
!         .. Local Scalars ..
          REAL (KIND=nag_wp)              :: anorm, cond, d1, eps, estnrm
          INTEGER                         :: i, icase, id, ifail, j, lda, n
!         .. Local Arrays ..
          REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), p(:), work(:), x(:)
          INTEGER, ALLOCATABLE            :: iwork(:)
!         .. Intrinsic Functions ..
          INTRINSIC                          max
!         .. Executable Statements ..
          WRITE (nout,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Example 1'
          WRITE (nout,*)
!         Skip heading in data file
          READ (nin,*)
          READ (nin,*)
          READ (nin,*)
          READ (nin,*) n
          lda = n
          ALLOCATE (a(lda,n),p(n),work(n),x(n),iwork(n))
          READ (nin,*) (a(i,1:n),i=1,n)
!         First compute the norm of A. 
!         DASUM (f06ekf) returns the sum of the absolute values of a column of A.
          anorm = zero
          DO j = 1, n
             anorm = max(anorm,dasum(n,a(1,j),1))
          END DO
          WRITE (nout,99999) 'Computed norm of A =', anorm

!         Next estimate the norm of inverse(A). We do not form the
!         inverse explicitly.

!         ifail: behaviour on error exit
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
!         Factorise A as P*A = L*U using F03AFF.
          CALL f03aff(n,eps,a,lda,d1,id,p,ifail)

          icase = 0
LOOP:     DO
             ifail = 0
             CALL f04ycf(icase,n,x,estnrm,work,iwork,ifail)

             IF (icase/=0) THEN
                IF (icase==1) THEN
!                  Return the vector inv(P*A)*X by solving the equations
!                  L*U*Y = X, overwriting Y on X. First solve L*Z = X for Z.
!                  The NAG name equivalent of dtrsv is f06pjf
                   CALL dtrsv('Lower','No Transpose','Non-Unit',n,a,lda,x,1)
!                  Then solve U*Y = Z for Y.
                   CALL dtrsv('Upper','No Transpose','Unit',n,a,lda,x,1)
                ELSE IF (icase==2) THEN
!                  Return the vector inv(P*A)'*X by solving U'*L'*Y = X,
!                  overwriting Y on X. First solve U'*Z = X for Z.
                   CALL dtrsv('Upper','Transpose','Unit',n,a,lda,x,1)
!                  Then solve L'*Y = Z for Y.
                   CALL dtrsv('Lower','Transpose','Non-Unit',n,a,lda,x,1)
                END IF
!            Continue until icase is returned as 0.
             ELSE
                WRITE (nout,99999) 'Estimated norm of inverse(A) =', estnrm
                cond = anorm*estnrm
                WRITE (nout,99998) 'Estimated condition number of A =', cond
                WRITE (nout,*)
                EXIT LOOP
             END IF
          END DO LOOP

99999     FORMAT (1X,A,F8.4)
99998     FORMAT (1X,A,F5.1)
       END SUBROUTINE ex1
       SUBROUTINE ex2

!         .. Use Statements ..
          USE nag_library, ONLY : f01brf, f04axf, f04ycf, nag_wp
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          REAL (KIND=nag_wp), PARAMETER   :: tenth = 0.1E+0_nag_wp
          REAL (KIND=nag_wp), PARAMETER   :: zero = 0.0E+0_nag_wp
          INTEGER, PARAMETER              :: nin = 5
!         .. Local Scalars ..
          REAL (KIND=nag_wp)              :: anorm, cond, estnrm, resid, sum, u
          INTEGER                         :: i, icase, ifail, j, licn, lirn,   &
                                             n, nz
          LOGICAL                         :: grow, lblock
!         .. Local Arrays ..
          REAL (KIND=nag_wp), ALLOCATABLE :: a(:), w(:), work1(:), x(:)
          INTEGER, ALLOCATABLE            :: icn(:), ikeep(:), irn(:), iw(:),  &
                                             iwork(:)
          INTEGER                         :: idisp(10)
          LOGICAL                         :: abort(4)
!         .. Intrinsic Functions ..
          INTRINSIC                          abs, max
!         .. Executable Statements ..
          WRITE (nout,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Example 2'
          WRITE (nout,*)
!         Skip heading in data file
          READ (nin,*)
          READ (nin,*)
!         Input N, the order of matrix A, and NZ, the number of non-zero
!         elements of A.
          READ (nin,*) n, nz
          licn = 4*nz
          lirn = 2*nz
          ALLOCATE (a(licn),w(n),work1(n),x(n),icn(licn),ikeep(5*n),irn(lirn), &
             iw(8*n),iwork(n))
!         Input the elements of A, along with row and column information.
          READ (nin,*) (a(i),irn(i),icn(i),i=1,nz)
!         First compute the norm of A.
          anorm = zero
          DO i = 1, n
             sum = zero
             DO j = 1, nz
                IF (icn(j)==i) sum = sum + abs(a(j))
             END DO
             anorm = max(anorm,sum)
          END DO
          WRITE (nout,99999) 'Computed norm of A =', anorm
!         Next estimate the norm of inverse(A). We do not form the
!         inverse explicitly.
!         Factorise A into L*U using F01BRF.
          u = tenth
          lblock = .TRUE.
          grow = .TRUE.
          abort(1) = .TRUE.
          abort(2) = .TRUE.
          abort(3) = .FALSE.
          abort(4) = .TRUE.

!         ifail: behaviour on error exit
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
          CALL f01brf(n,nz,a,licn,irn,lirn,icn,u,ikeep,iw,w,lblock,grow,abort, &
             idisp,ifail)

          icase = 0
LOOP:     DO
             ifail = 0
             CALL f04ycf(icase,n,x,estnrm,work1,iwork,ifail)

             IF (icase/=0) THEN
!               Return X := inv(A)*X or X = inv(A)'*X, depending on the
!               value of ICASE, by solving A*Y = X or A'*Y = X,
!               overwriting Y on X.
                CALL f04axf(n,a,licn,icn,ikeep,x,w,icase,idisp,resid)
!            Continue until icase is returned as 0.
             ELSE
                WRITE (nout,99999) 'Estimated norm of inverse(A) =', estnrm
                cond = anorm*estnrm
                WRITE (nout,99998) 'Estimated condition number of A =', cond
                EXIT LOOP
             END IF
          END DO LOOP

99999     FORMAT (1X,A,F8.4)
99998     FORMAT (1X,A,F5.1)
       END SUBROUTINE ex2
    END PROGRAM f04ycfe