Program f04zdfe

!     F04ZDF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: f04zdf, f06uaf, nag_wp, zgetrf, zgetrs
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: cond, nrma, nrminv
      Integer                          :: i, ifail, irevcm, lda, ldx, ldy, m,  &
                                          n, seed, t
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: a(:,:), work(:), x(:,:), y(:,:)
      Real (Kind=nag_wp), Allocatable  :: rwork(:)
      Real (Kind=nag_wp)               :: workr(1)
      Integer, Allocatable             :: ipiv(:), iwork(:)
!     .. Executable Statements ..
      Write (nout,*) 'F04ZDF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) m, n, t

      lda = m
      ldx = n
      ldy = m
      Allocate (a(lda,n))
      Allocate (x(ldx,t))
      Allocate (y(ldy,t))
      Allocate (work(m*t))
      Allocate (rwork(2*n))
      Allocate (iwork(2*n+5*t+20))
      Allocate (ipiv(n))

!     Read A from data file
      Read (nin,*)(a(i,1:n),i=1,m)

!     Compute 1-norm of A
      nrma = f06uaf('1',m,n,a,lda,workr)
      Write (nout,99999) 'The norm of A is: ', nrma

!     Estimate the norm of A^(-1) without forming A^(-1)
      irevcm = 0
      ifail = 0

      seed = 652

!     Perform an LU factorization so that A=LU where L and U are lower
!     and upper triangular.

!     The NAG name equivalent of zgetrf is f07arf
      Call zgetrf(m,n,a,lda,ipiv,ifail)

loop: Do
        Call f04zdf(irevcm,m,n,x,ldx,y,ldy,nrminv,t,seed,work,rwork,iwork,     &
          ifail)
        If (irevcm/=0) Then
          If (irevcm==1) Then
!           Compute Y = inv(A)*X

!           The NAG name equivalent of zgetrf is f07asf
            Call zgetrs('N',n,t,a,lda,ipiv,x,ldx,ifail)
!           X was overwritten by ZGETRS, so set Y=X
            y(1:n,1:t) = x(1:n,1:t)
          Else
!           Compute X = herm(inv(A))*Y

!           The NAG name equivalent of zgetrf is f07asf
            Call zgetrs('C',n,t,a,lda,ipiv,y,ldy,ifail)
!           Y was overwritten by ZGETRS, so set X=Y
            x(1:n,1:t) = y(1:n,1:t)
          End If
        Else

          Write (nout,99999) 'The estimated norm of inverse(A) is: ', nrminv

!         Compute and print the estimated condition number
          cond = nrminv*nrma
          Write (nout,*)
          Write (nout,99999) 'Estimated condition number of A: ', cond
          Write (nout,*)
          Exit loop
        End If
      End Do loop
99999 Format (1X,A,F6.2)

    End Program f04zdfe