* F04ZCF Example Program Text * Mark 17 Revised. NAG Copyright 1995. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX, KLMAX, KUMAX, LDA, LDX, NRHS PARAMETER (NMAX=8,KLMAX=NMAX-1,KUMAX=NMAX-1, + LDA=2*KLMAX+KUMAX+1,LDX=NMAX,NRHS=1) * .. Local Scalars .. DOUBLE PRECISION ANORM, COND, ESTNRM INTEGER I, ICASE, IFAIL, INFO, J, K, KL, KU, N * .. Local Arrays .. COMPLEX *16 A(LDA,NMAX), WORK(NMAX), X(LDX,NRHS) DOUBLE PRECISION RWORK(1) INTEGER IPIV(NMAX) * .. External Functions .. DOUBLE PRECISION F06UBF EXTERNAL F06UBF * .. External Subroutines .. EXTERNAL F04ZCF, ZGBTRF, ZGBTRS * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. Executable Statements .. WRITE (NOUT,*) 'F04ZCF Example Program Results' * Skip heading in data file READ (NIN,*) READ (NIN,*) N, KL, KU IF (N.LE.NMAX .AND. KL.LE.KLMAX .AND. KU.LE.KUMAX) THEN 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). We do not form the * inverse explicitly. * Factorise A into P*L*U. * CALL ZGBTRF(N,N,KL,KU,A,LDA,IPIV,INFO) * ICASE = 0 IFAIL = 0 20 CALL F04ZCF(ICASE,N,X,ESTNRM,WORK,IFAIL) * IF (ICASE.EQ.0) THEN WRITE (NOUT,99999) 'Estimated norm of inverse(A) =', ESTNRM ELSE IF (ICASE.EQ.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.EQ.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. GO TO 20 END IF COND = ANORM*ESTNRM WRITE (NOUT,99998) 'Estimated condition number of A =', COND END IF STOP * 99999 FORMAT (1X,A,F8.4) 99998 FORMAT (1X,A,F6.1) END