! D02GAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02gafe_mod ! Data for D02GAF example program ! .. 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 :: iset = 1, n = 3, nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: beta Contains Subroutine fcn(x,y,f) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: x ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: f(*) Real (Kind=nag_wp), Intent (In) :: y(*) ! .. Executable Statements .. f(1) = y(2) f(2) = y(3) f(3) = -y(1)*y(3) - beta*(1.0E0_nag_wp-y(2)*y(2)) Return End Subroutine fcn End Module d02gafe_mod Program d02gafe ! D02GAF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02gaf, nag_wp, x04abf Use d02gafe_mod, Only: beta, fcn, iset, n, nin, nout, one, zero ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: a, b, h, tol Integer :: i, ifail, j, k, liw, lw, mnp, & np, outchn ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: u(:,:), v(:,:), w(:), x(:), y(:,:) Integer, Allocatable :: iw(:) ! .. Intrinsic Procedures .. Intrinsic :: real ! .. Executable Statements .. Write (nout,*) 'D02GAF Example Program Results' ! Skip heading in data file Read (nin,*) ! n: number of differential equations ! mnp: maximum permitted number of mesh points. Read (nin,*) mnp liw = mnp*(2*n+1) + n*n + 4*n + 2 lw = mnp*(3*n*n+6*n+2) + 4*n*n + 4*n Allocate (iw(liw),u(n,2),v(n,2),w(lw),x(mnp),y(n,mnp)) ! tol: positive absolute error tolerance ! np : determines whether a default or user-supplied mesh is used. ! a : left-hand boundary point, b: right-hand boundary point. Read (nin,*) tol Read (nin,*) np Read (nin,*) a, b outchn = nout Call x04abf(iset,outchn) beta = zero u(1:n,1:2) = zero v(1:n,1:2) = zero v(1,2) = one v(3,1) = one v(3,2) = one u(2,2) = one u(1,2) = b x(1) = a h = (b-a)/real(np-1,kind=nag_wp) Do i = 2, np - 1 x(i) = x(i-1) + h End Do x(np) = b loop: Do k = 1, 2 Select Case (k) Case (1) beta = zero Case (2) beta = 0.2_nag_wp End Select ! ifail: behaviour on error exit ! =1 for quiet-soft exit ! * Set ifail to 111 to obtain monitoring information * ifail = 1 Call d02gaf(u,v,n,a,b,tol,fcn,mnp,x,y,np,w,lw,iw,liw,ifail) If (ifail>=0) Write (nout,99999) 'Problem with BETA = ', beta If (ifail==0 .Or. ifail==3) Then Write (nout,*) If (ifail==3) Write (nout,*) ' IFAIL = 3' Write (nout,99998) np Write (nout,99997) Write (nout,99996)(x(i),(y(j,i),j=1,n),i=1,np) beta = beta + 0.2E0_nag_wp Else Write (nout,99995) ifail Exit loop End If End Do loop 99999 Format (/1X,A,F7.2) 99998 Format (1X,'Solution on final mesh of ',I2,' points') 99997 Format (1X,' X(I) Y1(I) Y2(I) Y3(I)') 99996 Format (1X,F11.3,3F13.4) 99995 Format (1X/1X,' ** D02GAF returned with IFAIL = ',I5) End Program d02gafe