! D02RAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02rafe_mod ! D02RAF 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 :: two = 2.0_nag_wp Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp Integer, Parameter :: iset = 1, n = 3, nin = 5, nout = 6 Contains Subroutine fcn(x,eps,y,f,n) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: eps, x Integer, Intent (In) :: n ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: f(n) Real (Kind=nag_wp), Intent (In) :: y(n) ! .. Executable Statements .. f(1) = y(2) f(2) = y(3) f(3) = -y(1)*y(3) - two*(one-y(2)*y(2))*eps Return End Subroutine fcn Subroutine g(eps,ya,yb,bc,n) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: eps Integer, Intent (In) :: n ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: bc(n) Real (Kind=nag_wp), Intent (In) :: ya(n), yb(n) ! .. Executable Statements .. bc(1) = ya(1) bc(2) = ya(2) bc(3) = yb(2) - one Return End Subroutine g Subroutine jaceps(x,eps,y,f,n) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: eps, x Integer, Intent (In) :: n ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: f(n) Real (Kind=nag_wp), Intent (In) :: y(n) ! .. Executable Statements .. f(1:2) = zero f(3) = -two*(one-y(2)*y(2)) Return End Subroutine jaceps Subroutine jacgep(eps,ya,yb,bcep,n) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: eps Integer, Intent (In) :: n ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: bcep(n) Real (Kind=nag_wp), Intent (In) :: ya(n), yb(n) ! .. Executable Statements .. bcep(1:n) = zero Return End Subroutine jacgep Subroutine jacobf(x,eps,y,f,n) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: eps, x Integer, Intent (In) :: n ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: f(n,n) Real (Kind=nag_wp), Intent (In) :: y(n) ! .. Executable Statements .. f(1:n,1:n) = zero f(1,2) = one f(2,3) = one f(3,1) = -y(3) f(3,2) = two*two*y(2)*eps f(3,3) = -y(1) Return End Subroutine jacobf Subroutine jacobg(eps,ya,yb,aj,bj,n) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: eps Integer, Intent (In) :: n ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: aj(n,n), bj(n,n) Real (Kind=nag_wp), Intent (In) :: ya(n), yb(n) ! .. Executable Statements .. aj(1:n,1:n) = zero bj(1:n,1:n) = zero aj(1,1) = one aj(2,2) = one bj(3,2) = one Return End Subroutine jacobg End Module d02rafe_mod Program d02rafe ! D02RAF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02raf, nag_wp, x04abf Use d02rafe_mod, Only: fcn, g, iset, jaceps, jacgep, jacobf, jacobg, n, & nin, nout ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: deleps, tol Integer :: ifail, ijac, init, j, ldy, & liwork, lwork, mnp, np, numbeg, & nummix, outchn ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: abt(:), work(:), x(:), y(:,:) Integer, Allocatable :: iwork(:) ! .. Executable Statements .. Write (nout,*) 'D02RAF Example Program Results' ! Skip heading in data file Read (nin,*) Read (nin,*) mnp, np ldy = n liwork = mnp*(2*n+1) + n lwork = mnp*(3*n*n+6*n+2) + 4*n*n + 3*n Allocate (abt(n),work(lwork),x(mnp),y(ldy,mnp),iwork(liwork)) outchn = nout Write (nout,*) Call x04abf(iset,outchn) Read (nin,*) tol, deleps Read (nin,*) init, ijac, numbeg, nummix Read (nin,*) x(1), x(np) ! ifail: behaviour on error exit ! =1 for quiet-soft exit ! * Set IFAIL to 111 to obtain monitoring information * ifail = 1 Call d02raf(n,mnp,np,numbeg,nummix,tol,init,x,y,ldy,abt,fcn,g,ijac, & jacobf,jacobg,deleps,jaceps,jacgep,work,lwork,iwork,liwork,ifail) If (ifail==0 .Or. ifail==4) Then Write (nout,*) 'Calculation using analytic Jacobians' If (ifail==4) Write (nout,99996) 'On exit from D02RAF IFAIL = 4' Write (nout,*) Write (nout,99999) 'Solution on final mesh of ', np, ' points' Write (nout,*) ' X(I) Y1(I) Y2(I) Y3(I)' Write (nout,99998)(x(j),y(1:n,j),j=1,np) Write (nout,*) Write (nout,*) 'Maximum estimated error by components' Write (nout,99997) abt(1:n) Else Write (nout,99996) ' ** D02RAF returned with IFAIL = ', ifail End If 99999 Format (1X,A,I2,A) 99998 Format (1X,F10.3,3F13.4) 99997 Format (11X,1P,3E13.2) 99996 Format (1X,A,I5) End Program d02rafe