Program d03fafe ! D03FAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: d03faf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: dx, dy, dz, error, lambda, pertrb, & t, xf, xs, yf, ys, yy, zf, zs, zz Integer :: i, ifail, j, k, l, lbdcnd, ldf, & ldf2, lwrk, m, maxlm, mbdcnd, n, & nbdcnd ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: bdxf(:,:), bdxs(:,:), bdyf(:,:), & bdys(:,:), bdzf(:,:), bdzs(:,:), & f(:,:,:), w(:), x(:), y(:), z(:) ! .. Intrinsic Procedures .. Intrinsic :: abs, cos, real, sin ! .. Executable Statements .. Write (nout,*) 'D03FAF Example Program Results' ! Skip heading in data file Read (nin,*) Read (nin,*) l, m, n, maxlm ldf = l + 1 ldf2 = m + 1 lwrk = 2*(n+1)*maxlm + 3*l + 3*m + 4*n + 6000 Allocate (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)) Read (nin,*) lambda Read (nin,*) xs, xf Read (nin,*) ys, yf Read (nin,*) zs, zf Read (nin,*) lbdcnd, mbdcnd, nbdcnd ! Define the grid points for later use. dx = (xf-xs)/real(l,kind=nag_wp) Do i = 1, l + 1 x(i) = xs + real(i-1,kind=nag_wp)*dx End Do dy = (yf-ys)/real(m,kind=nag_wp) Do j = 1, m + 1 y(j) = ys + real(j-1,kind=nag_wp)*dy End Do dz = (zf-zs)/real(n,kind=nag_wp) Do k = 1, n + 1 z(k) = zs + real(k-1,kind=nag_wp)*dz End Do ! Define the array of derivative boundary values. Do j = 1, m + 1 yy = sin(y(j)) bdzf(1:l+1,j) = -yy*x(1:l+1)**4 End Do ! Note that for this example all other boundary arrays are ! dummy variables. ! We define the function boundary values in the F array. f(1,1:m+1,1:n+1) = 0.0_nag_wp Do k = 1, n + 1 zz = cos(z(k)) Do j = 1, m + 1 f(l+1,j,k) = zz*sin(y(j)) End Do End Do f(1:l+1,1:m+1,1) = -bdzf(1:l+1,1:m+1) ! Define the values of the right hand side of the Helmholtz ! equation. Do k = 2, n + 1 zz = 4.0_nag_wp*cos(z(k)) Do j = 1, m + 1 yy = sin(y(j))*zz Do i = 2, l f(i,j,k) = yy*x(i)**2*(3.0_nag_wp-x(i)**2) End Do End Do End Do ! Call D03FAF to generate and solve the finite difference equation. ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 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) ! Compute discretization error. The exact solution to the ! problem is ! U(X,Y,Z) = X**4*SIN(Y)*COS(Z) error = 0.0_nag_wp Do k = 1, n + 1 zz = cos(z(k)) Do j = 1, m + 1 yy = sin(y(j))*zz Do i = 1, l + 1 t = abs(f(i,j,k)-yy*x(i)**4) If (t>error) error = t End Do End Do End Do Write (nout,*) Write (nout,99999) error 99999 Format (1X,'Maximum component of discretization error =',1P,E13.6) End Program d03fafe