Program d03edfe ! D03EDF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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