* F12AEF 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) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D+0) * .. Local Scalars .. COMPLEX *16 C1, C2, C3, CSIG DOUBLE PRECISION DENI, DENR, NUMI, NUMR, SIGMAI, SIGMAR INTEGER IFAIL, IFAIL1, INFO, IREVCM, J, N, NCONV, NCV, + NEV, NITER, NSHIFT LOGICAL FIRST * .. Local Arrays .. COMPLEX *16 CDD(MAXN), CDL(MAXN), CDU(MAXN), CDU2(MAXN), + CTEMP(MAXN) DOUBLE PRECISION AX(MAXN), COMM(LCOMM), D(MAXNCV,3), MX(MAXN), + RESID(MAXN), V(LDV,MAXNCV), X(MAXN) INTEGER ICOMM(LICOMM), IPIV(MAXN) * .. External Functions .. DOUBLE PRECISION DDOT, DNRM2, F06BNF EXTERNAL DDOT, DNRM2, F06BNF * .. External Subroutines .. EXTERNAL AV, DCOPY, F12AAF, F12ABF, F12ACF, F12ADF, + F12AEF, MV, ZGTTRF, ZGTTRS * .. Intrinsic Functions .. INTRINSIC CMPLX, DBLE * .. Executable Statements .. WRITE (NOUT,*) 'F12AEF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) N, NEV, NCV, SIGMAR, SIGMAI 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 REAL',ICOMM,COMM,IFAIL) * Set problem type CALL F12ADF('GENERALIZED',ICOMM,COMM,IFAIL) * Solve A*x = lambda*B*x in shift-invert mode. * The shift, sigma, is a complex number (sigmar, sigmai). * OP = Real_Part{inv[A-(SIGMAR,SIGMAI)*M]*M and B = M. CSIG = CMPLX(SIGMAR,SIGMAI,KIND=KIND(DENI)) C1 = CMPLX(-2.0D0,0.0D0,KIND=KIND(DENI)) - CSIG C2 = CMPLX(2.0D0,0.0D0,KIND=KIND(DENI)) - 4.0D0*CSIG C3 = CMPLX(3.0D0,0.0D0,KIND=KIND(DENI)) - CSIG * DO 20 J = 1, N - 1 CDL(J) = C1 CDD(J) = C2 CDU(J) = C3 20 CONTINUE CDD(N) = C2 * CALL ZGTTRF(N,CDL,CDD,CDU,CDU2,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) THEN * Perform x <--- OP*x = inv[A-SIGMA*M]*M*x CALL MV(N,X) DO 60 J = 1, N CTEMP(J) = CMPLX(X(J),KIND=KIND(DENI)) 60 CONTINUE CALL ZGTTRS('N',N,1,CDL,CDD,CDU,CDU2,IPIV,CTEMP,N, + INFO) DO 80 J = 1, N X(J) = DBLE(CTEMP(J)) 80 CONTINUE ELSE IF (IREVCM.EQ.1) THEN * Perform x <--- OP*x = inv[A-SIGMA*M]*M*x, * M*X stored in MX. DO 100 J = 1, N CTEMP(J) = CMPLX(MX(J),KIND=KIND(DENI)) 100 CONTINUE CALL ZGTTRS('N',N,1,CDL,CDD,CDU,CDU2,IPIV,CTEMP,N, + INFO) DO 120 J = 1, N X(J) = DBLE(CTEMP(J)) 120 CONTINUE ELSE IF (IREVCM.EQ.2) THEN * Perform y <--- M*x CALL MV(N,X) ELSE IF (IREVCM.EQ.4) 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) FIRST = .TRUE. DO 140 J = 1, NCONV * Use Rayleigh Quotient to recover eigenvalues of the * original problem. IF (D(J,2).EQ.ZERO) THEN * Ritz value is real. CALL AV(N,V(1,J),AX) NUMR = DDOT(N,V(1,J),1,AX,1) CALL DCOPY(N,V(1,J),1,MX,1) CALL MV(N,MX) DENR = DDOT(N,V(1,J),1,MX,1) D(J,1) = NUMR/DENR ELSE IF (FIRST) THEN * Ritz value is complex. * Compute x'(Ax) CALL AV(N,V(1,J),AX) NUMR = DDOT(N,V(1,J),1,AX,1) NUMI = DDOT(N,V(1,J+1),1,AX,1) CALL AV(N,V(1,J+1),AX) NUMR = NUMR + DDOT(N,V(1,J+1),1,AX,1) NUMI = -NUMI + DDOT(N,V(1,J),1,AX,1) * Compute x'(Mx) CALL DCOPY(N,V(1,J),1,MX,1) CALL MV(N,MX) DENR = DDOT(N,V(1,J),1,MX,1) DENI = DDOT(N,V(1,J+1),1,MX,1) CALL DCOPY(N,V(1,J+1),1,MX,1) CALL MV(N,MX) DENR = DENR + DDOT(N,V(1,J+1),1,MX,1) DENI = -DENI + DDOT(N,V(1,J),1,MX,1) * d=x'(Ax)/x'(Mx) D(J,1) = (NUMR*DENR+NUMI*DENI)/F06BNF(DENR,DENI) D(J,2) = (NUMI*DENR-NUMR*DENI)/F06BNF(DENR,DENI) FIRST = .FALSE. ELSE * Second of complex conjugate pair. D(J,1) = D(J-1,1) D(J,2) = -D(J-1,2) FIRST = .TRUE. END IF 140 CONTINUE * Print computed eigenvalues. WRITE (NOUT,99996) NCONV, SIGMAR, SIGMAI DO 160 J = 1, NCONV WRITE (NOUT,99995) J, D(J,1), D(J,2) 160 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 =',E12.4) 99997 FORMAT (1X,' ** F12AAF returned with IFAIL = ',I5) 99996 FORMAT (1X,/' The ',I4,' generalized Ritz values closest to ','(', + F8.4,', ',F8.4,')',' are:',/) 99995 FORMAT (1X,I8,5X,'(',F7.4,',',F7.4,')') END * SUBROUTINE MV(N,V) * Compute the in-place matrix vector multiplication X<---M*X, * where M is mass matrix formed by using piecewise linear elements * on [0,1]. * * .. Parameters .. DOUBLE PRECISION FOUR PARAMETER (FOUR=4.0D+0) * .. Scalar Arguments .. INTEGER N * .. Array Arguments .. DOUBLE PRECISION V(N) * .. Local Scalars .. DOUBLE PRECISION VM1, VV INTEGER J * .. Executable Statements .. VM1 = V(1) V(1) = FOUR*V(1) + V(2) DO 20 J = 2, N - 1 VV = V(J) V(J) = VM1 + FOUR*VV + V(J+1) VM1 = VV 20 CONTINUE V(N) = VM1 + FOUR*V(N) RETURN END * SUBROUTINE AV(N,V,W) * .. Parameters .. DOUBLE PRECISION THREE, TWO PARAMETER (THREE=3.0D+0,TWO=2.0D+0) * .. Scalar Arguments .. INTEGER N * .. Array Arguments .. DOUBLE PRECISION V(N), W(N) * .. Local Scalars .. INTEGER J * .. Executable Statements .. W(1) = TWO*V(1) + THREE*V(2) DO 20 J = 2, N - 1 W(J) = -TWO*V(J-1) + TWO*V(J) + THREE*V(J+1) 20 CONTINUE W(N) = -TWO*V(N-1) + TWO*V(N) RETURN END