Program f04zcfe ! F04ZCF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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