* F12ABF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER LICOMM, NIN, NOUT PARAMETER (LICOMM=140,NIN=5,NOUT=6) INTEGER MAXN, MAXNCV, LDV PARAMETER (MAXN=256,MAXNCV=30,LDV=MAXN) INTEGER LCOMM PARAMETER (LCOMM=3*MAXN+3*MAXNCV*MAXNCV+6*MAXNCV+60) INTEGER IMON PARAMETER (IMON=0) DOUBLE PRECISION ONE, TWO PARAMETER (ONE=1.0D+0,TWO=2.0D+0) * .. Local Scalars .. DOUBLE PRECISION H, RHO, S, S1, S2, S3, SIGMAI, SIGMAR INTEGER I, IFAIL, IFAIL1, INFO, IREVCM, J, N, NCONV, NCV, + NEV, NITER, NSHIFT, NX * .. Local Arrays .. DOUBLE PRECISION COMM(LCOMM), D(MAXNCV,3), DD(MAXN), DL(MAXN), + DU(MAXN), DU2(MAXN), MX(MAXN), RESID(MAXN), + V(LDV,MAXNCV), X(MAXN) INTEGER ICOMM(LICOMM), IPIV(MAXN) * .. External Functions .. DOUBLE PRECISION DNRM2 EXTERNAL DNRM2 * .. External Subroutines .. EXTERNAL DGTTRF, DGTTRS, F12AAF, F12ABF, F12ACF, F12ADF, + F12AEF * .. Intrinsic Functions .. INTRINSIC DBLE * .. Executable Statements .. WRITE (NOUT,*) 'F12ABF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) NX, NEV, NCV, RHO, SIGMAR, SIGMAI N = NX*NX IF (N.LT.1 .OR. N.GT.MAXN) THEN WRITE (NOUT,99999) 'N is out of range: N = ', N ELSE IF (NCV.GT.MAXNCV) THEN WRITE (NOUT,99999) 'NCV is out of range: NCV = ', NCV ELSE IFAIL = 1 CALL F12AAF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL) * IF (IFAIL.EQ.0) THEN * Set the mode. IFAIL = 0 CALL F12ADF('SHIFTED INVERSE REAL',ICOMM,COMM,IFAIL) * Construct C = A - SIGMA*I, and factorize using DGTTRF/F07CDF. H = ONE/DBLE(N+1) S = RHO*H/TWO S1 = -ONE - S S2 = TWO - SIGMAR S3 = -ONE + S DO 20 J = 1, N - 1 DL(J) = S1 DD(J) = S2 DU(J) = S3 20 CONTINUE DD(N) = S2 * CALL DGTTRF(N,DL,DD,DU,DU2,IPIV,INFO) * IREVCM = 0 IFAIL = -1 40 CONTINUE CALL F12ABF(IREVCM,RESID,V,LDV,X,MX,NSHIFT,COMM,ICOMM,IFAIL) IF (IREVCM.NE.5) THEN IF (IREVCM.EQ.-1 .OR. IREVCM.EQ.1) THEN * Perform x <--- OP*x = inv[A-SIGMA*I]*x using DGTTRS/F07CEF. CALL DGTTRS('N',N,1,DL,DD,DU,DU2,IPIV,X,N,INFO) ELSE IF (IREVCM.EQ.4 .AND. IMON.NE.0) THEN * Output monitoring information CALL F12AEF(NITER,NCONV,D,D(1,2),D(1,3),ICOMM,COMM) WRITE (6,99998) NITER, NCONV, DNRM2(NEV,D(1,3),1) END IF GO TO 40 END IF IF (IFAIL.EQ.0) THEN * Post-Process using F12ACF to compute eigenvalues/vectors. IFAIL1 = 0 CALL F12ACF(NCONV,D,D(1,2),V,LDV,SIGMAR,SIGMAI,RESID,V, + LDV,COMM,ICOMM,IFAIL1) * Print computed eigenvalues. WRITE (NOUT,99996) NCONV DO 60 I = 1, NCONV WRITE (NOUT,99995) I, D(I,1), D(I,2) 60 CONTINUE END IF ELSE WRITE (NOUT,99997) IFAIL END IF END IF * 99999 FORMAT (1X,A,I5) 99998 FORMAT (1X,'Iteration',1X,I3,', No. converged =',1X,I3,', norm o', + 'f estimates =',E16.8) 99997 FORMAT (1X,' ** F12ABF returned with IFAIL = ',I5) 99996 FORMAT (1X,/' The ',I4,' Ritz values of closest to unity are:',/) 99995 FORMAT (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )') END