! D03EEF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (Out) :: alpha, beta, delta, epslon, & gamma, phi, psi Real (Kind=nag_wp), Intent (In) :: x, y ! .. Intrinsic Procedures .. 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) ! .. 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 Procedures .. 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) ! .. Function Return Value .. Real (Kind=nag_wp) :: fexact ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: x, y ! .. Intrinsic Procedures .. 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 Procedures .. 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