Program d03ecfe ! D03ECF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: d03ecf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: aparam, conchn, conres, root2, x1, & x2, y1, y2, yy, z1, z2 Integer :: i, ifail, itcoun, itmax, itused, & ixn, iyn, izn, j, k, lda, n1, n2, & n3, ndir, sda ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:,:), b(:,:,:), c(:,:,:), & chngs(:), d(:,:,:), e(:,:,:), & f(:,:,:), g(:,:,:), q(:,:,:), & resids(:), t(:,:,:), wrksp1(:,:,:), & wrksp2(:,:,:), wrksp3(:,:,:), & wrksp4(:,:,:), x(:), y(:), z(:) ! .. Intrinsic Procedures .. Intrinsic :: cos, exp, sqrt ! .. Executable Statements .. Write (nout,*) 'D03ECF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) Read (nin,*) n1, n2, n3, itmax lda = n1 sda = n2 Allocate (a(lda,sda,n3),b(lda,sda,n3),c(lda,sda,n3),chngs(itmax), & d(lda,sda,n3),e(lda,sda,n3),f(lda,sda,n3),g(lda,sda,n3),q(lda,sda,n3), & resids(itmax),t(lda,sda,n3),wrksp1(lda,sda,n3),wrksp2(lda,sda,n3), & wrksp3(lda,sda,n3),wrksp4(lda,sda,n3),x(n1),y(n2),z(n3)) Read (nin,*) x(1:n1) Read (nin,*) y(1:n2) Read (nin,*) z(1:n3) Read (nin,*) conres, conchn Read (nin,*) ndir root2 = sqrt(two) aparam = one itcoun = 0 ! Set up difference equation coefficients, source terms and ! initial approximation. a(1:n1,1:n2,1:n3) = zero b(1:n1,1:n2,1:n3) = zero c(1:n1,1:n2,1:n3) = zero d(1:n1,1:n2,1:n3) = zero e(1:n1,1:n2,1:n3) = zero f(1:n1,1:n2,1:n3) = zero g(1:n1,1:n2,1:n3) = zero q(1:n1,1:n2,1:n3) = zero t(1:n1,1:n2,1:n3) = zero ! Non-zero Specification for internal nodes Do k = 2, n3 - 1 Do j = 2, n2 - 1 Do i = 2, n1 - 1 a(i,j,k) = two/((z(k)-z(k-1))*(z(k+1)-z(k-1))) g(i,j,k) = two/((z(k+1)-z(k))*(z(k+1)-z(k-1))) b(i,j,k) = two/((y(j)-y(j-1))*(y(j+1)-y(j-1))) f(i,j,k) = two/((y(j+1)-y(j))*(y(j+1)-y(j-1))) c(i,j,k) = two/((x(i)-x(i-1))*(x(i+1)-x(i-1))) e(i,j,k) = two/((x(i+1)-x(i))*(x(i+1)-x(i-1))) d(i,j,k) = -a(i,j,k) - b(i,j,k) - c(i,j,k) - e(i,j,k) - f(i,j,k) - & g(i,j,k) End Do End Do End Do ! Non-zero specification for boundary nodes yy = one/y(n2) x1 = (x(1)+one)*yy x2 = (x(n1)+one)*yy Do j = 1, n2 y1 = root2*y(j)*yy q(1,j,1:n3) = exp(x1)*cos(y1)*exp((-z(1:n3)-one)*yy) q(n1,j,1:n3) = exp(x2)*cos(y1)*exp((-z(1:n3)-one)*yy) End Do y1 = root2*y(1)*yy y2 = root2*y(n2)*yy Do i = 1, n1 x1 = (x(i)+one)*yy q(i,1,1:n3) = exp(x1)*cos(y1)*exp((-z(1:n3)-one)*yy) q(i,n2,1:n3) = exp(x1)*cos(y2)*exp((-z(1:n3)-one)*yy) End Do z1 = (-z(1)-one)*yy z2 = (-z(n3)-one)*yy Do i = 1, n1 x1 = (x(i)+one)*yy q(i,1:n2,1) = exp(x1)*cos(root2*y(1:n2)*yy)*exp(z1) q(i,1:n2,n3) = exp(x1)*cos(root2*y(1:n2)*yy)*exp(z2) End Do ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d03ecf(n1,n2,n3,lda,sda,a,b,c,d,e,f,g,q,t,aparam,itmax,itcoun, & itused,ndir,ixn,iyn,izn,conres,conchn,resids,chngs,wrksp1,wrksp2, & wrksp3,wrksp4,ifail) Write (nout,*) 'Iteration Maximum Maximum' Write (nout,*) ' number residual change' If (itused/=0) Then Write (nout,99999)(i,resids(i),chngs(i),i=1,itused) End If Write (nout,*) Write (nout,*) 'Table of calculated function values' Write (nout,*) Write (nout,99998) Do k = 1, n3 Do j = 1, n2 Write (nout,99997) k, j, (i,t(i,j,k),i=1,n1) End Do End Do 99999 Format (2X,I3,9X,E11.4,4X,E11.4) 99998 Format (1X,'K J ',4(' (I T )')) 99997 Format (1X,I1,2X,I1,1X,4(1X,I3,2X,F8.3)) End Program d03ecfe