* D03FAF Example Program Text * Mark 14 Release. NAG Copyright 1989. * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) INTEGER L, M, N, MAXLM, LDF, LDF2, LWRK PARAMETER (L=16,M=32,N=20,MAXLM=32,LDF=L+1,LDF2=M+1, + LWRK=2*(N+1)*MAXLM+3*L+3*M+4*N+6000) * .. Local Scalars .. DOUBLE PRECISION DX, DY, DZ, ERROR, LAMBDA, PERTRB, PI, T, XF, XS, + YF, YS, ZF, ZS INTEGER I, IFAIL, J, K, LBDCND, MBDCND, NBDCND * .. Local Arrays .. DOUBLE PRECISION BDXF(LDF2,N+1), BDXS(LDF2,N+1), BDYF(LDF,N+1), + BDYS(LDF,N+1), BDZF(LDF,M+1), BDZS(LDF,M+1), + F(LDF,LDF2,N+1), W(LWRK), X(L+1), Y(M+1), Z(N+1) * .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF * .. External Subroutines .. EXTERNAL D03FAF * .. Intrinsic Functions .. INTRINSIC ABS, COS, DBLE, SIN * .. Executable Statements .. WRITE (NOUT,*) 'D03FAF Example Program Results' LAMBDA = -2.0D0 XS = 0.0D0 XF = 1.0D0 LBDCND = 1 YS = 0.0D0 PI = X01AAF(PI) YF = 2.0D0*PI MBDCND = 0 ZS = 0.0D0 ZF = PI/2.0D0 NBDCND = 2 * * Define the grid points for later use. * DX = (XF-XS)/DBLE(L) DO 20 I = 1, L + 1 X(I) = XS + DBLE(I-1)*DX 20 CONTINUE DY = (YF-YS)/DBLE(M) DO 40 J = 1, M + 1 Y(J) = YS + DBLE(J-1)*DY 40 CONTINUE DZ = (ZF-ZS)/DBLE(N) DO 60 K = 1, N + 1 Z(K) = ZS + DBLE(K-1)*DZ 60 CONTINUE * * Define the array of derivative boundary values. * DO 100 J = 1, M + 1 DO 80 I = 1, L + 1 BDZF(I,J) = -X(I)**4*SIN(Y(J)) 80 CONTINUE 100 CONTINUE * * Note that for this example all other boundary arrays are * dummy variables. * * We define the function boundary values in the F array. * DO 140 K = 1, N + 1 DO 120 J = 1, M + 1 F(1,J,K) = 0.0D0 F(L+1,J,K) = SIN(Y(J))*COS(Z(K)) 120 CONTINUE 140 CONTINUE DO 180 J = 1, M + 1 DO 160 I = 1, L + 1 F(I,J,1) = X(I)**4*SIN(Y(J)) 160 CONTINUE 180 CONTINUE * * Define the values of the right hand side of the Helmholtz * equation. * DO 240 K = 2, N + 1 DO 220 J = 1, M + 1 DO 200 I = 2, L F(I,J,K) = 4.0D0*X(I)**2*(3.0D0-X(I)**2)*SIN(Y(J)) + *COS(Z(K)) 200 CONTINUE 220 CONTINUE 240 CONTINUE * * Call D03FAF to generate and solve the finite difference equation. * IFAIL = 1 * CALL D03FAF(XS,XF,L,LBDCND,BDXS,BDXF,YS,YF,M,MBDCND,BDYS,BDYF,ZS, + ZF,N,NBDCND,BDZS,BDZF,LAMBDA,LDF,LDF2,F,PERTRB,W,LWRK, + IFAIL) IF (IFAIL.EQ.0) THEN * * Compute discretization error. The exact solution to the * problem is * * U(X,Y,Z) = X**4*SIN(Y)*COS(Z) * ERROR = 0.0D0 DO 300 K = 1, N + 1 DO 280 J = 1, M + 1 DO 260 I = 1, L + 1 T = ABS(F(I,J,K)-X(I)**4*SIN(Y(J))*COS(Z(K))) IF (T.GT.ERROR) ERROR = T 260 CONTINUE 280 CONTINUE 300 CONTINUE WRITE (NOUT,*) WRITE (NOUT,99999) ERROR ELSE WRITE (NOUT,99998) IFAIL END IF * 99999 FORMAT (1X,'Maximum component of discretization error =',1P,E13.6) 99998 FORMAT (1X,/1X,' ** D03FAF returned with IFAIL = ',I5) END