Program d03edfe
! D03EDF Example Program Text
! Mark 27.0 Release. NAG Copyright 2019.
! .. Use Statements ..
Use nag_library, Only: d03edf, nag_wp, x04abf
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: four = 4.0_nag_wp
Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp
Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp
Integer, Parameter :: iset = 1, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp) :: acc, alpha, hx, hy
Integer :: i, ifail, iout, ix, iy1, iy2, j, k, &
lda, levels, maxit, ngx, ngy, numit, &
outchn
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:,:), rhs(:), u(:), ub(:), us(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*) 'D03EDF Example Program Results'
! Skip heading in data file
Read (nin,*)
Read (nin,*) levels
ngx = 2**levels + 1
ngy = ngx
lda = 4*(ngx+1)*(ngy+1)/3
Allocate (a(lda,7),rhs(lda),u(lda),ub(ngx*ngy),us(lda))
outchn = nout
Write (nout,*)
Read (nin,*) alpha, acc
Read (nin,*) maxit
Call x04abf(iset,outchn)
! ** Set iout > 2 to obtain intermediate output **
iout = 0
hx = one/real(ngx+1,kind=nag_wp)
hy = one/real(ngy+1,kind=nag_wp)
! Set up operator, right-hand side and initial guess for
! step-lengths HX and HY
a(1:ngx*ngy,1) = one - half*alpha
a(1:ngx*ngy,2) = half*alpha
a(1:ngx*ngy,3) = one - half*alpha
a(1:ngx*ngy,4) = -four + alpha
a(1:ngx*ngy,5) = one - half*alpha
a(1:ngx*ngy,6) = half*alpha
a(1:ngx*ngy,7) = one - half*alpha
rhs(1:ngx*ngy) = -four*hx*hy
ub(1:ngx*ngy) = zero
! Correction for the boundary conditions
! Horizontal boundaries --
! Boundary condition on Y=0 -- U=0
a(2:ngx-1,1:2) = zero
! Boundary condition on Y=1 -- U=0
ix = (ngy-1)*ngx
a(ix+1:ix+ngx-1,6:7) = zero
! Vertical boundaries --
iy1 = 1
iy2 = ngx
Do j = 2, ngy - 1
! Boundary condition on X=0 -- U=0
iy1 = iy1 + ngx
a(iy1,3) = zero
a(iy1,6) = zero
! Boundary condition on X=1 -- U=1
iy2 = iy2 + ngx
rhs(iy2) = rhs(iy2) - a(iy2,5) - a(iy2,2)
a(iy2,2) = zero
a(iy2,5) = zero
End Do
! Now the four corners --
! Bottom left corner
a(1,1:3) = zero
a(1,6) = zero
! Top left corner
a(ix+1,3) = zero
a(ix+1,6:7) = zero
! Bottom right corner
! Use average value at discontinuity ( = 0.5 )
k = ngx
rhs(k) = rhs(k) - a(k,2)*half - a(k,5)
a(k,1:2) = zero
a(k,5) = zero
! Top right corner
k = ngx*ngy
rhs(k) = rhs(k) - a(k,2) - a(k,5)
a(k,2) = zero
a(k,5:7) = zero
! Solve the equations
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03edf(ngx,ngy,lda,a,rhs,ub,maxit,acc,us,u,iout,numit,ifail)
Write (nout,99999) ngx, ngy, acc, maxit
Write (nout,*)
Write (nout,99998) 'Residual norm =', us(1)
Write (nout,99997) 'Number of iterations =', numit
Write (nout,*)
Write (nout,*) 'Solution'
Write (nout,*)
Write (nout,99996) ' I/J', (i,i=1,ngx)
Do j = 1, ngy
Write (nout,99995) j, (u(i+(j-1)*ngx),i=1,ngx)
End Do
99999 Format (1X,'NGX = ',I3,' NGY = ',I3,' ACC =',1P,E10.2,' MAXIT',' = ', &
I3)
99998 Format (1X,A,1P,E12.2)
99997 Format (1X,A,I5)
99996 Format (1X,A,10I7,:)
99995 Format (1X,I3,2X,10F7.3,:)
End Program d03edfe