! D03EEF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d03eefe_mod ! D03EEF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: zero = 0.0_nag_wp INTEGER, PARAMETER :: nin = 5, nout = 6 CONTAINS SUBROUTINE pdef(x,y,alpha,beta,gamma,delta,epslon,phi,psi) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: alpha, beta, delta, epslon, & gamma, phi, psi REAL (KIND=nag_wp), INTENT (IN) :: x, y ! .. Intrinsic Functions .. INTRINSIC cos, sin ! .. Executable Statements .. alpha = one beta = zero gamma = one delta = 50.0_nag_wp epslon = 50.0_nag_wp phi = zero psi = sin(x)*((-alpha-gamma+phi)*sin(y)+epslon*cos(y)) + & cos(x)*(delta*sin(y)+beta*cos(y)) RETURN END SUBROUTINE pdef SUBROUTINE bndy(x,y,a,b,c,ibnd) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: bottom = 0, left = 3, & right = 1, top = 2 ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: a, b, c REAL (KIND=nag_wp), INTENT (IN) :: x, y INTEGER, INTENT (IN) :: ibnd ! .. Intrinsic Functions .. INTRINSIC sin ! .. Executable Statements .. IF (ibnd==top .OR. ibnd==right) THEN ! Solution prescribed a = one b = zero c = sin(x)*sin(y) ELSE IF (ibnd==bottom) THEN ! Derivative prescribed a = zero b = one c = -sin(x) ELSE IF (ibnd==left) THEN ! Derivative prescribed a = zero b = one c = -sin(y) END IF RETURN END SUBROUTINE bndy FUNCTION fexact(x,y) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: fexact ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x, y ! .. Intrinsic Functions .. INTRINSIC sin ! .. Executable Statements .. fexact = sin(x)*sin(y) RETURN END FUNCTION fexact END MODULE d03eefe_mod PROGRAM d03eefe ! D03EEF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d03edf, d03eef, nag_wp USE d03eefe_mod, ONLY : bndy, fexact, nin, nout, pdef, zero ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: acc, hx, hy, rmserr, xmax, xmin, & xx, ymax, ymin, yy INTEGER :: i, icase, ifail, iout, ix, j, & lda, levels, maxit, ngx, ngxy, & ngy, numit CHARACTER (7) :: scheme ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), rhs(:), u(:), ub(:), & us(:), x(:), y(:) ! .. Intrinsic Functions .. INTRINSIC real, sqrt ! .. Executable Statements .. WRITE (nout,*) 'D03EEF Example Program Results' WRITE (nout,*) FLUSH (nout) ! Skip heading in data file READ (nin,*) READ (nin,*) levels ngx = 2**levels + 1 ngy = ngx lda = 4*(ngx+1)*(ngy+1)/3 ngxy = ngx*ngy ALLOCATE (a(lda,7),rhs(lda),u(lda),ub(ngxy),us(lda),x(ngxy),y(ngxy)) READ (nin,*) xmin, xmax READ (nin,*) ymin, ymax hx = (xmax-xmin)/real(ngx-1,kind=nag_wp) DO i = 1, ngx xx = xmin + real(i-1,kind=nag_wp)*hx x(i:ngxy:ngx) = xx END DO hy = (ymax-ymin)/real(ngy-1,kind=nag_wp) DO j = 1, ngy yy = ymin + real(j-1,kind=nag_wp)*hy y((j-1)*ngx+1:j*ngx) = yy END DO ! ** set iout > 2 to obtain intermediate output from D03EDF ** iout = 0 READ (nin,*) acc READ (nin,*) maxit CASES: DO icase = 1, 2 SELECT CASE (icase) CASE (1) ! Central differences scheme = 'Central' CASE (2) ! Upwind differences scheme = 'Upwind' END SELECT ! Discretize the equations ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = -1 CALL d03eef(xmin,xmax,ymin,ymax,pdef,bndy,ngx,ngy,lda,a,rhs,scheme, & ifail) IF (ifail<0) THEN WRITE (nout,99995) ifail EXIT CASES END IF ! Set the initial guess to zero ub(1:ngxy) = zero ! Solve the equations ifail = 0 CALL d03edf(ngx,ngy,lda,a,rhs,ub,maxit,acc,us,u,iout,numit,ifail) ! Print out the solution WRITE (nout,*) WRITE (nout,*) 'Exact solution above computed solution' WRITE (nout,*) WRITE (nout,99998) ' I/J', (i,i=1,ngx) rmserr = zero DO j = ngy, 1, -1 ix = (j-1)*ngx WRITE (nout,*) WRITE (nout,99999) j, (fexact(x(ix+i),y(ix+i)),i=1,ngx) WRITE (nout,99999) j, u(ix+1:ix+ngx) DO i = 1, ngx rmserr = rmserr + (fexact(x(ix+i),y(ix+i))-u(ix+i))**2 END DO END DO rmserr = sqrt(rmserr/real(ngxy,kind=nag_wp)) WRITE (nout,*) WRITE (nout,99997) 'Number of Iterations = ', numit WRITE (nout,99996) 'RMS Error = ', rmserr END DO CASES 99999 FORMAT (1X,I3,2X,10F7.3:/(6X,10F7.3)) 99998 FORMAT (1X,A,10I7:/(6X,10I7)) 99997 FORMAT (1X,A,I3) 99996 FORMAT (1X,A,1P,E10.2) 99995 FORMAT (1X,' ** D03EEF returned with IFAIL = ',I5) END PROGRAM d03eefe