Program d03ubfe
! D03UBF Example Program Text
! Mark 26.1 Release. NAG Copyright 2017.
! .. 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 Procedures ..
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