Program d03uafe ! D03UAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: d03uaf, 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 Integer :: i, ifail, it, j, lda, n1, n2, nits ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), b(:,:), c(:,:), d(:,:), & e(:,:), q(:,:), r(:,:), t(:,:), & wrksp1(:,:), wrksp2(:,:), x(:), y(:) ! .. Intrinsic Procedures .. Intrinsic :: abs, cos, exp, max, real ! .. Executable Statements .. Write (nout,*) 'D03UAF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) Read (nin,*) n1, n2, nits lda = n1 Allocate (a(lda,n2),b(lda,n2),c(lda,n2),d(lda,n2),e(lda,n2),q(lda,n2), & r(lda,n2),t(lda,n2),wrksp1(lda,n2),wrksp2(lda,n2),x(n1),y(n2)) Read (nin,*) x(1:n1) Read (nin,*) y(1:n2) aparam = one ! Set up difference equation coefficients, source terms and ! initial S a(1:n1,1:n2) = zero b(1:n1,1:n2) = zero d(1:n1,1:n2) = zero e(1:n1,1:n2) = zero q(1:n1,1:n2) = zero t(1:n1,1:n2) = zero ! Specification for internal nodes Do j = 2, n2 - 1 a(2:n1-1,j) = two/((y(j)-y(j-1))*(y(j+1)-y(j-1))) e(2:n1-1,j) = two/((y(j+1)-y(j))*(y(j+1)-y(j-1))) End Do Do i = 2, n1 - 1 b(i,2:n2-1) = two/((x(i)-x(i-1))*(x(i+1)-x(i-1))) d(i,2:n2-1) = two/((x(i+1)-x(i))*(x(i+1)-x(i-1))) End Do c(1:n1,1:n2) = -a(1:n1,1:n2) - b(1:n1,1:n2) - d(1:n1,1:n2) - & e(1:n1,1:n2) ! Specification for boundary nodes Do j = 1, n2 q(1,j) = exp((x(1)+one)/y(n2))*cos(y(j)/y(n2)) q(n1,j) = exp((x(n1)+one)/y(n2))*cos(y(j)/y(n2)) End Do Do i = 1, n1 q(i,1) = exp((x(i)+one)/y(n2))*cos(y(1)/y(n2)) q(i,n2) = exp((x(i)+one)/y(n2))*cos(y(n2)/y(n2)) End Do ! Iterative loop Do it = 1, nits ! Calculate the residuals resmax = zero resmn = zero Do j = 1, n2 Do i = 1, n1 If (c(i,j)/=zero) Then ! Five point molecule formula r(i,j) = q(i,j) - a(i,j)*t(i,j-1) - b(i,j)*t(i-1,j) - & c(i,j)*t(i,j) - d(i,j)*t(i+1,j) - e(i,j)*t(i,j+1) Else ! Explicit equation r(i,j) = q(i,j) - t(i,j) End If ares = abs(r(i,j)) resmax = max(resmax,ares) resmn = resmn + ares End Do End Do resmn = resmn/(real(n1*n2,kind=nag_wp)) ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d03uaf(n1,n2,lda,a,b,c,d,e,aparam,it,r,wrksp1,wrksp2,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 j = 1, n2 Do i = 1, n1 t(i,j) = t(i,j) + r(i,j) adel = abs(r(i,j)) delmax = max(delmax,adel) delmn = delmn + adel End Do End Do delmn = delmn/(real(n1*n2,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,*) Write (nout,99995) 'I', 1, (i,i=2,6) Write (nout,*) ' J' Do j = 1, n2 Write (nout,99998) j, (t(i,j),i=1,n1) End Do 99999 Format (1X,I3,4(2X,E11.4)) 99998 Format (1X,I2,1X,6(F9.3,2X)) 99997 Format (1X,A,6X,A,19X,A) 99996 Format (3X,A,7X,A,8X,A,11X,A,6X,A/) 99995 Format (4X,A,4X,I1,5I11) End Program d03uafe