PROGRAM d03ubfe ! D03UBF Example Program Text ! Mark 23 Release. NAG Copyright 2011. ! .. Use Statements .. USE nag_library, ONLY : d03ubf, 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) :: adel, aparam, ares, delmax, delmn, & resmax, resmn, root2, x1, x2, y1, & y2, yy, z1, z2 INTEGER :: i, ifail, it, j, k, lda, n1, n2, n3, & nits, sda ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:,:), b(:,:,:), c(:,:,:), & d(:,:,:), e(:,:,:), f(:,:,:), & g(:,:,:), q(:,:,:), r(:,:,:), & t(:,:,:), wrksp1(:,:,:), & wrksp2(:,:,:), wrksp3(:,:,:), x(:), & y(:), z(:) ! .. Intrinsic Functions .. INTRINSIC abs, cos, exp, max, real, sqrt ! .. Executable Statements .. WRITE (nout,*) 'D03UBF Example Program Results' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) n1, n2, n3, nits lda = n1 sda = n2 ALLOCATE (a(lda,sda,n3),b(lda,sda,n3),c(lda,sda,n3),d(lda,sda,n3), & e(lda,sda,n3),f(lda,sda,n3),g(lda,sda,n3),q(lda,sda,n3), & r(lda,sda,n3),t(lda,sda,n3),wrksp1(lda,sda,n3),wrksp2(lda,sda,n3), & wrksp3(lda,sda,n3),x(n1),y(n2),z(n3)) READ (nin,*) x(1:n1) READ (nin,*) y(1:n2) READ (nin,*) z(1:n3) root2 = sqrt(two) aparam = one ! 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 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 ! Specification for internal nodes DO k = 2, n3 - 1 a(2:n1-1,2:n2-1,k) = two/((z(k)-z(k-1))*(z(k+1)-z(k-1))) g(2:n1-1,2:n2-1,k) = two/((z(k+1)-z(k))*(z(k+1)-z(k-1))) END DO DO j = 2, n2 - 1 b(2:n1-1,j,2:n3-1) = two/((y(j)-y(j-1))*(y(j+1)-y(j-1))) f(2:n1-1,j,2:n3-1) = two/((y(j+1)-y(j))*(y(j+1)-y(j-1))) END DO DO i = 2, n1 - 1 c(i,2:n2-1,2:n3-1) = two/((x(i)-x(i-1))*(x(i+1)-x(i-1))) e(i,2:n2-1,2:n3-1) = two/((x(i+1)-x(i))*(x(i+1)-x(i-1))) END DO d(1:n1,1:n2,1:n3) = -a(1:n1,1:n2,1:n3) - b(1:n1,1:n2,1:n3) - & c(1:n1,1:n2,1:n3) - e(1:n1,1:n2,1:n3) - f(1:n1,1:n2,1:n3) - & g(1:n1,1:n2,1:n3) ! 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 ! Iterative loop DO it = 1, nits resmax = zero resmn = zero DO k = 1, n3 DO j = 1, n2 DO i = 1, n1 IF (d(i,j,k)/=zero) THEN ! Seven point molecule formula r(i,j,k) = q(i,j,k) - a(i,j,k)*t(i,j,k-1) - & b(i,j,k)*t(i,j-1,k) - c(i,j,k)*t(i-1,j,k) - & d(i,j,k)*t(i,j,k) - e(i,j,k)*t(i+1,j,k) - & f(i,j,k)*t(i,j+1,k) - g(i,j,k)*t(i,j,k+1) ELSE ! Explicit equation r(i,j,k) = q(i,j,k) - t(i,j,k) END IF ares = abs(r(i,j,k)) resmax = max(resmax,ares) resmn = resmn + ares END DO END DO END DO resmn = resmn/(real(n1*n2*n3,kind=nag_wp)) ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 CALL d03ubf(n1,n2,n3,lda,sda,a,b,c,d,e,f,g,aparam,it,r,wrksp1, & wrksp2,wrksp3,ifail) IF (it==1) THEN WRITE (nout,99997) 'Iteration', 'Residual', 'Change' WRITE (nout,99996) 'No', 'Max.', 'Mean', 'Max.', 'Mean' END IF ! Update the dependent variable delmax = zero delmn = zero DO k = 1, n3 DO j = 1, n2 DO i = 1, n1 t(i,j,k) = t(i,j,k) + r(i,j,k) adel = abs(r(i,j,k)) delmax = max(delmax,adel) delmn = delmn + adel END DO END DO END DO delmn = delmn/(real(n1*n2*n3,kind=nag_wp)) WRITE (nout,99999) it, resmax, resmn, delmax, delmn ! Convergence tests here if required END DO ! End of iterative loop WRITE (nout,*) WRITE (nout,*) 'Table of calculated function values' WRITE (nout,99995) WRITE (nout,*) WRITE (nout,99998) ((k,j,(i,t(i,j,k),i=1,n1),j=1,n2),k=1,n3) 99999 FORMAT (1X,I5,4(2X,E11.4)) 99998 FORMAT ((1X,I1,I3,1X,4(1X,I3,2X,F8.3))) 99997 FORMAT (1X,A,6X,A,19X,A) 99996 FORMAT (2X,A,7X,A,8X,A,11X,A,6X,A/) 99995 FORMAT (1X,'K J',2X,4(1X,'(I T )')) END PROGRAM d03ubfe