    Program f04ycfe

!     F04YCF Example Program Text

!     Mark 24 Release. NAG Copyright 2012.

!     .. 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, dgetrf, dgetrs, f04ycf, nag_wp
!       .. Parameters ..
        Real (Kind=nag_wp), Parameter    :: zero = 0.0E+0_nag_wp
        Integer, Parameter               :: nin = 5
!       .. Local Scalars ..
        Real (Kind=nag_wp)               :: anorm, cond, estnrm
        Integer                          :: i, icase, ifail, info, j, lda, n
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable  :: a(:,:), p(:), work(:), x(:)
        Integer, Allocatable             :: ipiv(:), iwork(:)
!       .. Intrinsic Procedures ..
        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),ipiv(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
!       LU Factorize A
!       The NAG name equivalent of dgetrf is f07adf
        Call dgetrf(n,n,a,lda,ipiv,info)

        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(A)*X 
!             The NAG name equivalent of dgetrs is f07aef
              Call dgetrs('N',n,1,a,lda,ipiv,x,n,info)
            Else If (icase==2) Then
!             Return the vector inv(A)'*X 
              Call dgetrs('T',n,1,a,lda,ipiv,x,n,info)
            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
!       .. 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 Procedures ..
        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
