Program f04zcfe

!     F04ZCF Example Program Text

!     Mark 25 Release. NAG Copyright 2014.

!     .. Use Statements ..
      Use nag_library, Only: f04zcf, f06ubf, nag_wp, zgbtrf, zgbtrs
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: anorm, cond, estnrm
      Integer                          :: i, icase, ifail, info, j, k, kl, ku, &
                                          lda, ldx, n, nrhs
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: a(:,:), work(:), x(:,:)
      Real (Kind=nag_wp)               :: rwork(1)
      Integer, Allocatable             :: ipiv(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, min
!     .. Executable Statements ..
      Write (nout,*) 'F04ZCF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, kl, ku, nrhs
      lda = 2*kl + ku + 1
      ldx = n
      Allocate (a(lda,n),work(n),x(ldx,nrhs),ipiv(n))
      k = kl + ku + 1
      Read (nin,*)((a(k+i-j,j),j=max(i-kl,1),min(i+ku,n)),i=1,n)

!     First compute the 1-norm of A.
      anorm = f06ubf('1-norm',n,kl,ku,a(kl+1,1),lda,rwork)
      Write (nout,*)
      Write (nout,99999) 'Computed norm of A =', anorm

!     Next estimate the 1-norm of inverse(A).
!     Factorise A into P*L*U.
!     The NAG name equivalent of zgbtrf is f07brf
      Call zgbtrf(n,n,kl,ku,a,lda,ipiv,info)

      icase = 0

loop: Do
!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call f04zcf(icase,n,x,estnrm,work,ifail)

        If (icase/=0) Then
!         The NAG name equivalent of the backsolve routine zgbtrs is f07bsf
          If (icase==1) Then
!           Return X := inv(A)*X by solving A*Y = X, overwriting
!           Y on X.
            Call zgbtrs('No transpose',n,kl,ku,nrhs,a,lda,ipiv,x,ldx,info)
          Else If (icase==2) Then
!           Return X := conjg(inv(A)')*X by solving conjg(A')*Y
!           = X, overwriting Y on X.
            Call zgbtrs('Conjugate transpose',n,kl,ku,nrhs,a,lda,ipiv,x,ldx, &
              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
          Exit loop
        End If
      End Do loop

99999 Format (1X,A,F8.4)
99998 Format (1X,A,F6.1)
    End Program f04zcfe