!------------------------------------------------------------------ ! ! ABSTRACT: ! ! The Fortran 90 code IRKC is intended for the time integration of ! systems of partial differential equations (PDEs) of diffusion- ! reaction type for which the reaction Jacobian has real (negative) ! eigenvalues. It is based on a family of implicit-explicit ! Runge-Kutta-Chebyshev methods which are unconditionally stable for ! reaction terms and which impose a stability constraint associated ! with the diffusion terms that is quadratic in the number of stages. ! Special properties of the family make it possible for the code to ! select at each step the most efficient stable method as well as ! the most efficient step size. Moreover, they make it possible to ! apply the methods using just a few vectors of storage. ! A further step towards minimal storage requirements and optimal ! efficiency is achieved by exploiting the fact that the implicit ! terms, originating from the stiff reactions, are not coupled over ! the spatial grid points. Hence, the systems to be solved have a ! small dimension (viz., equal to the number of PDEs). ! These characteristics of the code make it especially attractive ! for problems in several spatial variables. ! IRKC is a successor to the RKC code (J. Comput. Appl. Math. 88 ! (1997), pp. 315-326) that solves similar problems without ! stiff reaction terms. !------------------------------------------------------------------- ! ! USAGE: ! ! IRKC integrates equations of the form ! ! y'(t) = F_E(t,y) + F_I(t,y), (1) ! ! when the reactions described by the F_I term cause the ODEs to be ! very stiff and the reaction Jacobian F'_I possesses a spectrum ! that is both real and negative. ! In IRKC, the term F_E is handled by a variant of the explicit ! methods of RKC and the term F_I is handled implicitly. Such ! schemes are known as IMEX methods. Implicit treatment of large ! systems of ODEs can have serious implications for storage, so we ! make assumptions about the structure of F_I that allow us to ! solve a useful class of problems with minimal storage. We have in ! mind systems of ODEs that arise from semi-discretization of ! diffusion-reaction PDEs. We make no assumption about how the ! discretization is done, but we do assume that the implicit terms ! at one grid point are not coupled to those at other grid points. ! This kind of decoupling arises naturally with finite differences, ! finite volumes, and discontinuous Galerkin discretizations. To ! take advantage of the great reductions in storage possible ! because of this assumption and the numerical method employed, ! the user must present the problem to the solver in an appropriate ! way. If there are NPDES partial differential equations, each call ! to the subroutine that evaluates F_I is made with a vector y of ! NPDES components that correspond to a single grid point. On ! demand the subroutine is also to evaluate the Jacobian of F_I, a ! matrix that is only NPDES by NPDES. ! ---------------------------------------------------------------- ! ! SOFTWARE ISSUES: ! ! We have redesigned the user interface of RKC and exploited the ! capabilities of Fortran 90 to make solving ODEs easier despite ! solving a larger class of problems. We begin here with an outline ! of the collection of programs that make up the package IRKC ! and then provide some details. Fortran 90 makes it possible to ! encapsulate all information about the solution as a single ! structure of a derived type. This approach and the dynamic storage ! facilities of Fortran 90 make it possible to relieve the user of ! tedious details of allocating storage. The solution structure can ! be called anything, but let us suppose that it is called SOL. ! It must be declared as type IRKC_SOL. This structure and the ! computation itself are initialized by a call to a function ! IRKC_SET. Among other things, the user must specify the initial ! point and a final point in this call. We exploit the possibility ! in Fortran 90 of optional arguments so that the most common ! options can be set by default. The integration is done by the ! subroutine IRKC. The default in IRKC_SET is for the integration ! to proceed a step at a time, but optionally the solver can return ! just the solution at the final point. After each step an ! auxiliary function IRKC_VAL can be used to approximate the ! solution anywhere in the span of the step. Statistics are ! available directly in fields of the solution structure, but they ! can all be displayed conveniently by calling the auxiliary ! subroutine IRKC_STATS. !------------------------------------------------------------------ ! ! FORM OF THE SOLUTION: ! ! The solution structure SOL is initialized by IRKC_SET. It is then ! an input argument of the solver IRKC. If the integration is being ! done a step at a time (the default), the SOL returned at one step ! is input to IRKC for the next step. The solution structure has a ! good many fields. The most interesting fields are the current ! value of the independent variable, SOL%T, and the current ! approximate solution, SOL%Y, a vector of NEQN components. The ! logical quantity SOL%DONE is monitored to learn when the ! integration is complete. ! The methods of IRKC provide a solution between mesh points that ! is evaluated in an auxiliary function. An approximate solution ! YOUT at a point TOUT in the span of the last step is obtained by ! YOUT = IRKC_VAL(SOL,TOUT). Some of the data held in fields of SOL ! for this purpose might be of direct interest, viz., the current ! approximation to the first derivative of the solution, SOL%YP, ! and the size of the last step taken, SOL%HLAST. ! ! A convenient way to see all the statistics is to CALL ! IRKC_STATS(SOL). Individual statistics are available as the ! integer fields ! ! SOL%NFE : number of evaluations of F_E, ! SOL%NFI : number of evaluations of F_I, ! SOL%NSTEPS : number of step, ! SOL%NACCPT : number of accepted step, ! SOL%NREJCT : number of rejected steps, ! SOL%NFESIG : number of function evaluations used in estimating ! the spectral radius of F'_E, ! SOL%MAXM : maximum number of stages used. !----------------------------------------------------------------- ! ! SPECIFICATION OF THE TASK: ! ! Only a few quantities are required to specify the task and ! initialize the integration. This is done with a call ! ! SOL = IRKC_SET(T0,Y0,TEND). ! ! This says that the integration is to start at T0 with a vector Y0 ! of NEQN components as initial value and go to TEND. It also says ! that the solution structure is to be called SOL. These arguments ! must appear, and in the order shown, but the remaining, optional ! arguments can follow in any order because they are specified using ! keywords. The solver must be told how many PDEs there are. This is ! 1 by default and any other value must be supplied with the keyword ! NPDES. For example, if there are 3 PDEs, the call above is changed ! to ! ! SOL = IRKC_SET(T0,Y0,TEND,NPDES=3). ! ! The most commonly used options are the error tolerances. RKC ! provides for a scalar relative error tolerance and a vector of ! absolute error tolerances. To leading order the storage required ! is a multiple of NEQN. A guiding principle of RKC, and especially ! IRKC, is minimize this multiple. Quite often users have scaled the ! variables so that a scalar absolute error tolerance is ! appropriate. If not, they can always rescale the variables so as ! to do away with the need for an array of NEQN absolute error ! tolerances. Accordingly, we have chosen to implement only a scalar ! absolute error tolerance in IRKC. This decision makes the solver ! easier to use and simplifies the program. The default relative ! error tolerance is 10^{-2} and the default absolute error ! tolerance is 10^{-3}. The call just illustrated causes the ! solver to use these default values. Other values are imposed with ! the keywords RE and AE. For example, to use a relative error ! tolerance of 10^{-3} and the default absolute error tolerance, ! the call is changed to ! ! SOL = IRKC_SET(T0,Y0,TEND,NPDES=3,RE=1D-3). ! ! By default IRKC takes one step at a time, but it can be instructed ! to return only after reaching TEND by giving the keyword ONE_STEP ! the value .FALSE. The methods for handling the explicit term F_E ! involve an approximate bound on the spectral radius of the Jacobian ! F'_E. This can be supplied with a function in a manner described ! below, but the default is to have the solver compute a bound. ! The cost of doing this can be reduced substantially when the ! Jacobian is constant by giving the keyword CONSTANT_J the value ! .TRUE. Sometimes it is useful to impose a maximum step size. ! This is done with the keyword HMAX. IRKC selects an initial step ! size automatically, but a value can be supplied with the keyword H0. ! ! It is sometimes useful to reset quantities and continue integrating. ! This is done by replacing the initial data T0,Y0 in the call list of ! IRKC_SET with the solution structure SOL. If the required argument ! TEND or any of the optional arguments are changed, the new values ! replace those stored in SOL and SOL%DONE is changed to .FALSE. ! Of course, NPDES cannot be changed. For example, we can instruct the ! solver to quit returning at every step with ! ! SOL = IRKC_SET(SOL,TEND,ONE_STEP=.FALSE.) ! ! and then call the solver again to continue on to TEND. ! ! Unless instructed to the contrary, the solver will print a message ! and stop when a fatal error occurs. This response to a fatal error ! can be changed with the optional argument STOP_ON_ERROR. Setting it ! to .FALSE. will cause the solver to return with a positive value of ! the integer SOL%ERR_FLAG that indicates the nature of the error. ! In this way a user can prevent an unwelcome termination of the run. ! Of course, if this report of a fatal error is ignored and the solver ! is called again with a positive value of SOL%ERR_FLAG, it will print ! a message and stop. !---------------------------------------------------------------------- ! ! DEFINING THE EQUATIONS: ! ! The differential equations are supplied as subroutines to IRKC. A ! typical invocation of IRKC looks like ! ! CALL IRKC(SOL,F_E,F_I). ! ! All three arguments in this call are required. The solution ! structure SOL is used for both input and output. F_E is the name ! of a subroutine for evaluating the explicit part of the right-hand ! side of the system of ODEs in (1). It is supplied just as for RKC, ! i.e., there is a subroutine of the form F_E(NEQN,T,Y,DY) that ! accepts a vector Y of length NEQN, evaluates F_E(t,y), and ! returns it as a vector DY of length NEQN. ! As explained previously, the term to be handled implicitly must ! be defined in a particular way that we now discuss more fully. ! ! Let NPDES be the number of PDEs. The ODEs must be coded in blocks ! that correspond to grid points. Specifically, ! for I = 1, NPDES+1, 2*NPDES+1, ..., NEQN-NPDES+1, the equations with ! indices I, I+1, ..., I+NPDES-1 must correspond to a single grid point. ! Accordingly there is to be a subroutine of the form ! F_I(GRID_POINT,NPDES,T,YG,DYG,WANT_JAC,JAC). For an index ! I = GRID_POINT, it evaluates F_I(T,YG) for YG = Y(I:I+NPDES-1). ! This value is returned as a vector DYG of NPDES components. The ! solver also requires the Jacobian of this function, a matrix JAC of ! NPDES by NPDES entries. It is convenient to have this computation in ! the subroutine where F_I is evaluated, but it does not have to be ! done every time that F_I is evaluated. The logical variable WANT_JAC ! is used to inform the subroutine when JAC is to be formed along with ! DYG. ! ! By default the solver approximates the spectral radius of the ! Jacobian F'_E using a nonlinear power method. This is convenient ! and often works well, but it can fail. When it is not too much ! trouble to supply a function that returns an upper bound for this ! spectral radius of roughly the correct size, this should be done ! because the integration is then faster, more reliable, and uses ! less storage. This function must have the form SR(NEQN,T,Y). Its ! name is supplied to the solver as an optional fourth argument, but ! in the circumstances there is no point to using a keyword, so the ! call has the form ! ! CALL IRKC(SOL,F_E,F_I,SR). !------------------------------------------------------------------- ! ! Authors: L.F. Shampine ! Mathematics Department ! Southern Methodist University ! Dallas, Texas 75275-0156 ! USA ! e-mail: lshampin@mail.smu.edu ! ! B.P. Sommeijer and J.G. Verwer ! Centre for Mathematics and Computer Science (CWI) ! Kruislaan 413 ! 1098 SJ Amsterdam ! The Netherlands ! e-mail: bsom@cwi.nl, Jan.Verwer@cwi.nl ! ! Details of the methods used and the performance of IRKC can be ! found in ! ! L.F. Shampine. B.P. Sommeijer and J.G. Verwer ! IRKC: an IMEX Solver for Stiff Diffusion-Reaction PDEs. ! ! J. Comput. Appl. Math. 196 (2006), pp. 485 - 497. ! Technical Report MAS-E0513, CWI, Amsterdam, 2005 ! (http://db.cwi.nl/rapporten/index.php?jaar=2005&dept=13) ! ! This source code for IRKC and some examples can also be obtained ! by anonymous ftp from the address ftp://ftp.cwi.nl/pub/bsom/irkc. !--------------------------------------------------------------------- ! MODULE IRKC_M IMPLICIT NONE ! Declare everything in the module to be PRIVATE unless specifically ! declared as PUBLIC. PRIVATE TYPE, PUBLIC :: IRKC_SOL INTEGER :: NEQN,NPDES DOUBLE PRECISION :: T,TEND,HLAST,RTOL,ATOL,HMAX,H0 DOUBLE PRECISION, DIMENSION(:), POINTER :: Y,YLAST,YP,YPLAST INTEGER :: NFE,NFI,NSTEPS,NACCPT,NREJCT,NFESIG,MAXM,ERR_FLAG LOGICAL :: DONE,ONE_STEP,CONSTANT_J,STOP_ON_ERROR END TYPE IRKC_SOL ! A structure SOL of type IRKC_SOL contains information ! about the numerical solution. Of most direct interest ! are the fields that hold the current value t of the ! independent variable and the NEQN components of the ! solution y(t). ! ! SOL%T -- t ! ! SOL%Y(*) -- y(t) (vector of NEQN components) ! ! The logical quantity SOL%DONE indicates whether ! the integration is complete. ! ! After a step by IRKC to SOL%T of size SOL%HLAST, the ! solution structure SOL contains all the information ! needed by the IRKC_VAL function to evaluate an ! approximate solution YARG(*) for any argument ARG in ! the interval [SOL%T-SOL%HLAST, SOL%T] that has the ! same order of accuracy as SOL%Y(*). ! ! Statistics about the integration are available as ! the integer fields ! ! SOL%NFE -- number of evaluations of F_E ! SOL%NFI -- number of evaluations of F_I ! SOL%NSTEPS -- number of steps ! SOL%NACCPT -- number of accepted steps ! SOL%NREJCT -- number of rejected steps ! SOL%NFESIG -- number of evaluations of F_E used ! in estimating the spectral radius ! SOL%MAXM -- maximum number of stages used ! ! They can be displayed conveniently with the auxiliary ! subroutine IRKC_STATS. INTERFACE IRKC_SET MODULE PROCEDURE SET1,SET2 END INTERFACE ! Global variables within module DOUBLE PRECISION :: UROUND=EPSILON(1D0),SQRTU LOGICAL :: HAVE_SPCRAD,STOP_ON_ERROR INTEGER :: NFE,NFI,NSTEPS,NACCPT,NREJCT,NFESIG,MAXM,MMAX DOUBLE PRECISION :: TDIR,HMAX,HLAST,RTOL,ATOL,HMUS1 ! Externally visible subroutines and functions: PUBLIC :: IRKC, IRKC_SET, IRKC_VAL, IRKC_STATS CONTAINS !******************************************************* !*** PUBLIC SUBROUTINES/FUNCTIONS *** !******************************************************* SUBROUTINE IRKC(SOL,F_E,F_I,SPCRAD) TYPE(IRKC_SOL) :: SOL DOUBLE PRECISION, OPTIONAL :: SPCRAD EXTERNAL F_E,F_I,SPCRAD INTEGER :: NEQN,NPDES DOUBLE PRECISION :: T,TEND DOUBLE PRECISION, DIMENSION(SOL%NEQN) :: YNEW LOGICAL :: RHO_CONVG_FAIL,IT_CONVG_FAIL,ONE_STEP,CONSTANT_J,& LAST_STEP INTEGER :: M,I,IER,ERR_FLAG DOUBLE PRECISION :: JACNRM,EST,ERR,FAC,H,HMIN,TEMP1,TEMP2,TEMP3 ! In the one-step mode, a few quantities need to be kept ! from one call to the next: INTEGER, SAVE :: NSTSIG DOUBLE PRECISION, SAVE :: HOLD,ERROLD,ABSH,SPRAD DOUBLE PRECISION, ALLOCATABLE, SAVE :: EV(:) LOGICAL, SAVE :: NEWSPC,JACATT IF (SOL%DONE) THEN PRINT *,' Cannot call IRKC after integration is done.' STOP END IF ERR_FLAG = SOL%ERR_FLAG IF ( (SOL%NFE > 0) .AND. (ERR_FLAG > 0) ) THEN PRINT *,' Cannot call IRKC after a fatal error.' STOP END IF ! Local variables for quantities in SOL NEQN = SOL%NEQN NPDES = SOL%NPDES T = SOL%T TEND = SOL%TEND TDIR = SIGN(1D0,TEND - T) HMAX = SOL%HMAX RTOL = SOL%RTOL ATOL = SOL%ATOL ONE_STEP = SOL%ONE_STEP CONSTANT_J = SOL%CONSTANT_J STOP_ON_ERROR = SOL%STOP_ON_ERROR HLAST = SOL%HLAST NFE = SOL%NFE NFI = SOL%NFI NSTEPS = SOL%NSTEPS NACCPT = SOL%NACCPT NREJCT = SOL%NREJCT NFESIG = SOL%NFESIG MAXM = SOL%MAXM ! Load a work array. YNEW = SOL%Y ! Initialize on the first call. IF (NFE == 0) THEN IF ((RTOL > 0.1D0) .OR. (RTOL < 10D0*UROUND)) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(1) ELSE SOL%ERR_FLAG = 1 RETURN END IF END IF IF (ATOL <= 0D0 ) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(2) ELSE SOL%ERR_FLAG = 2 RETURN END IF END IF HAVE_SPCRAD = PRESENT(SPCRAD) IF (.NOT. HAVE_SPCRAD) THEN ALLOCATE(EV(NEQN),STAT=IER) IF (IER /= 0) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(7) ELSE SOL%ERR_FLAG = 7 RETURN END IF END IF END IF SQRTU = SQRT(UROUND) MMAX = NINT(SQRT(RTOL/(10D0*UROUND))) MMAX = MAX(MMAX,2) JACATT = .FALSE. NSTSIG = 0 NEWSPC = .TRUE. CALL F_E(NEQN,T,YNEW,SOL%YP) NFE = NFE + 1 ! Get maximum norm of Jacobian of F_I for computation ! of initial step size. CALL ADD_F_I(NEQN,T,YNEW,SOL%YP,NPDES,F_I,JACNRM) HMIN = 10D0*UROUND*MAX(ABS(T),ABS(TEND)) END IF ! Start of loop for taking a step: TAKE_A_STEP: DO ! Estimate the spectral radius of the Jacobian ! when NEWSPC = .TRUE. A convergence failure ! in RHO is reported by RHO_CONVG_FAIL. IF (NEWSPC) THEN IF (HAVE_SPCRAD) THEN SPRAD = SPCRAD(NEQN,T,YNEW) ELSE CALL RHO(F_E,NEQN,T,YNEW,SOL%YP,SOL%YLAST,SOL%YPLAST,EV,& SPRAD,RHO_CONVG_FAIL) IF (RHO_CONVG_FAIL) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(4) ELSE SOL%ERR_FLAG = 4 RETURN END IF END IF ! The value SOL%YP returned by RHO corresponds to F_E only. CALL ADD_F_I(NEQN,T,YNEW,SOL%YP,NPDES,F_I) END IF JACATT = .TRUE. END IF ! Compute an initial step size. IF (NSTEPS == 0) THEN IF (SOL%H0 > 0D0) THEN ABSH = MIN(SOL%H0,HMAX) ELSE ABSH = HMAX IF (MAX(SPRAD,JACNRM)*ABSH > 1D0) THEN ABSH = 1D0/MAX(SPRAD,JACNRM) END IF END IF ABSH = MAX(ABSH,HMIN) SOL%YLAST = YNEW + ABSH*SOL%YP CALL F_E(NEQN,T+ABSH,SOL%YLAST,SOL%YPLAST) NFE = NFE + 1 CALL ADD_F_I(NEQN,T+ABSH,SOL%YLAST,SOL%YPLAST,NPDES,F_I) EST = RMSNORM( (SOL%YPLAST - SOL%YP) / (ATOL + RTOL*ABS(YNEW)) ) EST = ABSH*SQRT(EST/NEQN) IF (0.1D0*ABSH < HMAX*SQRT(EST)) THEN ABSH = MAX(0.1D0*ABSH/SQRT(EST),HMIN) ELSE ABSH = HMAX END IF END IF ! Adjust step size and determine number of stages M. LAST_STEP = .FALSE. IF (1.1D0*ABSH >= ABS(TEND - T)) THEN ABSH = ABS(TEND - T) LAST_STEP = .TRUE. END IF M = 1 + INT(SQRT(1.54D0*ABSH*SPRAD + 1D0)) ! Limit M to MMAX to control the growth of roundoff error. IF (M > MMAX) THEN M = MMAX ABSH = (M**2 - 1)/(1.54D0*SPRAD) LAST_STEP = .FALSE. END IF MAXM = MAX(M,MAXM) ! A tentative solution at T+H is returned in ! SOL%Y and its slope is evaluated in SOL%YLAST. H = TDIR*ABSH HMIN = 10D0*UROUND*MAX(ABS(T),ABS(T+H)) CALL STEP(NEQN,F_E,NPDES,F_I,T,YNEW,& H,M,SOL%Y,SOL%YLAST,SOL%YPLAST,& IT_CONVG_FAIL,ERR_FLAG) IF (ERR_FLAG > 0) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(ERR_FLAG) ELSE SOL%ERR_FLAG = ERR_FLAG RETURN END IF END IF IF (IT_CONVG_FAIL) THEN ERR = 10D0 ! Dummy value to trigger failed step. ELSE ! Here YNEW is y(T) and SOL%Y is the tentative y(T+H). ! Reuse of storage has y'(T) stored in SOL%YP and ! y'(T+H) is to be stored in SOL%YLAST. CALL F_E(NEQN,T+H,SOL%Y,SOL%YLAST) NFE = NFE + 1 CALL ERR_EST(NEQN,NPDES,T,H,YNEW,& SOL%YP,SOL%Y,SOL%YLAST,F_I,ERR,ERR_FLAG) IF (ERR_FLAG > 0) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(ERR_FLAG) ELSE SOL%ERR_FLAG = ERR_FLAG RETURN END IF END IF END IF NSTEPS = NSTEPS + 1 IF (ERR > 1D0) THEN ! Step is rejected. NREJCT = NREJCT + 1 IF (IT_CONVG_FAIL) THEN ABSH = ABSH/2D0 ELSE ABSH = 0.8D0*ABSH/SQRT(ERR) END IF IF (ABSH < HMIN) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(3) ELSE SOL%ERR_FLAG = 3 RETURN END IF ELSE NEWSPC = .NOT. JACATT CYCLE TAKE_A_STEP END IF END IF ! Step is accepted. NACCPT = NACCPT + 1 T = T + H JACATT = CONSTANT_J NSTSIG = MOD(NSTSIG+1,25) NEWSPC = .FALSE. IF (HAVE_SPCRAD .OR. (NSTSIG == 0)) THEN NEWSPC = .NOT. JACATT END IF ! Update stored quantities. HLAST = H DO I = 1,NEQN TEMP3 = SOL%YLAST(I) SOL%YLAST(I) = YNEW(I) YNEW(I) = SOL%Y(I) SOL%YPLAST(I) = SOL%YP(I) SOL%YP(I) = TEMP3 END DO ! Compute new step size. FAC = 10D0 IF (NACCPT == 1) THEN TEMP2 = SQRT(ERR) IF (0.8D0 < FAC*TEMP2) THEN FAC = 0.8D0/TEMP2 END IF ELSE TEMP1 = 0.8D0*ABSH*SQRT(ERROLD) TEMP2 = ABS(HOLD)*ERR IF (TEMP1 < FAC*TEMP2) THEN FAC = TEMP1/TEMP2 END IF END IF ABSH = MAX(0.1D0,FAC)*ABSH ABSH = MAX(HMIN,MIN(HMAX,ABSH)) ERROLD = ERR HOLD = H H = TDIR*ABSH IF (LAST_STEP .OR. ONE_STEP) THEN SOL%DONE = LAST_STEP EXIT TAKE_A_STEP END IF ! End of loop for taking a step: END DO TAKE_A_STEP ! Load working variables into SOL for return SOL%T = T SOL%HLAST = HLAST SOL%NFE = NFE SOL%NFI = NFI SOL%NSTEPS = NSTEPS SOL%NACCPT = NACCPT SOL%NREJCT = NREJCT SOL%NFESIG = NFESIG SOL%MAXM = MAXM IF (SOL%DONE .AND. (.NOT. HAVE_SPCRAD) ) THEN DEALLOCATE(EV,STAT=IER) IF (IER /= 0) THEN IF (STOP_ON_ERROR) THEN CALL REPORT_ERRORS(7) ELSE SOL%ERR_FLAG = 7 RETURN END IF END IF END IF END SUBROUTINE IRKC !--------------BEGIN GENERIC FUNCTION IRKC_SET-------------- FUNCTION SET1(T0,Y0,TEND,RE,AE,ONE_STEP,CONSTANT_J,& STOP_ON_ERROR,NPDES,HMAX,H0) RESULT(SOL) DOUBLE PRECISION :: T0,Y0(:),TEND INTEGER, OPTIONAL :: NPDES DOUBLE PRECISION, OPTIONAL :: RE,AE,HMAX,H0 LOGICAL, OPTIONAL :: ONE_STEP,CONSTANT_J,STOP_ON_ERROR TYPE(IRKC_SOL) :: SOL INTEGER :: NEQN,IER SOL%T = T0 SOL%TEND = TEND NEQN = SIZE(Y0) SOL%NEQN = NEQN ALLOCATE(SOL%Y(NEQN),SOL%YP(NEQN),SOL%YLAST(NEQN),& SOL%YPLAST(NEQN),STAT=IER) IF (IER /= 0) THEN PRINT *,'A storage allocation error occurred in IRKC_SET.' STOP END IF SOL%Y = Y0 SOL%YP = 0D0 SOL%YLAST = 0D0 SOL%YPLAST = 0D0 IF (PRESENT(RE)) THEN SOL%RTOL = RE ELSE SOL%RTOL = 1D-2 END IF IF (PRESENT(AE)) THEN SOL%ATOL = AE ELSE SOL%ATOL = 1D-3 END IF IF (PRESENT(ONE_STEP)) THEN SOL%ONE_STEP = ONE_STEP ELSE SOL%ONE_STEP = .TRUE. END IF IF (PRESENT(CONSTANT_J)) THEN SOL%CONSTANT_J = CONSTANT_J ELSE SOL%CONSTANT_J = .FALSE. END IF IF (PRESENT(STOP_ON_ERROR)) THEN SOL%STOP_ON_ERROR = STOP_ON_ERROR ELSE SOL%STOP_ON_ERROR = .TRUE. END IF SOL%ERR_FLAG = 0 IF (PRESENT(NPDES)) THEN SOL%NPDES = NPDES ELSE SOL%NPDES = 1 END IF IF (PRESENT(HMAX)) THEN SOL%HMAX = HMAX ELSE SOL%HMAX = ABS(TEND - T0) END IF IF (PRESENT(H0)) THEN SOL%H0 = ABS(H0) ELSE SOL%H0 = -1D0 END IF SOL%NFE = 0 SOL%NFI = 0 SOL%NSTEPS = 0 SOL%NACCPT = 0 SOL%NREJCT = 0 SOL%NFESIG = 0 SOL%MAXM = 0 SOL%HLAST = 0D0 SOL%DONE = .FALSE. END FUNCTION SET1 FUNCTION SET2(SOLIN,TEND,RE,AE,ONE_STEP,CONSTANT_J,& STOP_ON_ERROR,NPDES,HMAX,H0) RESULT(SOL) TYPE(IRKC_SOL) :: SOLIN,SOL DOUBLE PRECISION :: TEND DOUBLE PRECISION, OPTIONAL :: RE,AE,HMAX,H0 LOGICAL, OPTIONAL :: ONE_STEP,CONSTANT_J,STOP_ON_ERROR INTEGER, OPTIONAL :: NPDES IF (PRESENT(NPDES)) THEN IF (NPDES /= SOLIN%NPDES) THEN PRINT *,' Cannot change NPDES.' STOP END IF END IF ! Is this a new TEND? IF (ABS(TEND - SOLIN%TEND) > 0D0) THEN ! Is this a change of direction? IF (SIGN(1D0,TEND - SOLIN%TEND)*SOLIN%HLAST < 0) THEN PRINT *,' Cannot change direction without restarting.' STOP END IF SOLIN%TEND = TEND ! Reset in SOLIN for later copy to SOL. END IF SOL = SOLIN IF (PRESENT(RE)) SOL%RTOL = RE IF (PRESENT(AE)) SOL%ATOL = AE IF (PRESENT(ONE_STEP)) SOL%ONE_STEP = ONE_STEP IF (PRESENT(CONSTANT_J)) SOL%CONSTANT_J = CONSTANT_J IF (PRESENT(STOP_ON_ERROR)) SOL%STOP_ON_ERROR = STOP_ON_ERROR IF (PRESENT(HMAX)) SOL%HMAX = HMAX IF (PRESENT(H0)) SOL%H0 = ABS(H0) SOL%DONE = .FALSE. END FUNCTION SET2 !----------------END GENERIC FUNCTION IRKC_SET-------------- FUNCTION IRKC_VAL(SOL,ARG) RESULT(YARG) ! IRKC_VAL is used to compute approximate solutions at specific t ! and to compute cheaply the large number of approximations that ! may be needed for plotting or locating when events occur. ! After a step by IRKC to SOL%T of size SOL%HLAST, the solution ! structure SOL contains all the information needed by IRKC_VAL to ! evaluate an approximate solution YARG(*) for any argument ARG ! in the interval [SOL%T-SOL%HLAST, SOL%T] that has the same order ! of accuracy as SOL%Y(*). TYPE(IRKC_SOL) :: SOL DOUBLE PRECISION :: ARG DOUBLE PRECISION, DIMENSION(SOL%NEQN) :: YARG DOUBLE PRECISION :: TLAST,A1,A2,B1,B2,S TLAST = SOL%T - SOL%HLAST S = (ARG - TLAST) / SOL%HLAST A1 = (1D0 + 2D0*S)*(S - 1D0)**2 A2 = (3D0 - 2D0*S)*S**2 B1 = SOL%HLAST*S*(S - 1D0)**2 B2 = SOL%HLAST*(S - 1D0)*S**2 YARG = A1*SOL%YLAST + A2*SOL%Y + B1*SOL%YPLAST + B2*SOL%YP END FUNCTION IRKC_VAL SUBROUTINE IRKC_STATS(SOL) TYPE(IRKC_SOL) :: SOL INTEGER :: ITEMP ITEMP = DBLE(SOL%NFI)/( DBLE(SOL%NEQN)/DBLE(SOL%NPDES) ) PRINT *,' ' PRINT *,' Integration statistics: ' PRINT *,' ' PRINT *,' Evaluations of F_E = ', & SOL%NFE PRINT *,' Evaluations of F_I / (NEQN/NPDES) = ', & ITEMP PRINT *,' Number of steps = ', & SOL%NSTEPS PRINT *,' Number of accepted steps = ', & SOL%NACCPT PRINT *,' Number of rejected steps = ', & SOL%NREJCT PRINT *,' Evaluations of F_E for spectral radius = ', & SOL%NFESIG PRINT *,' Maximum number of stages used = ', & SOL%MAXM PRINT *,' ' END SUBROUTINE IRKC_STATS !******************************************************* !*** PRIVATE SUBROUTINES/FUNCTIONS *** !******************************************************* SUBROUTINE ADD_F_I(NEQN,T,Y,YP,NPDES,F_I,JACNRM) ! Evaluate the terms to be handled implicitly ! and add in blocks of size NPDEs. Optionally ! compute the maximum norm of the Jacobian. INTEGER :: NEQN,NPDES,I DOUBLE PRECISION :: T,Y(NEQN),YP(NEQN) DOUBLE PRECISION, OPTIONAL :: JACNRM DOUBLE PRECISION :: TEMP(NPDES),JAC(NPDES,NPDES) EXTERNAL F_I IF (.NOT. PRESENT(JACNRM)) THEN DO I = 1,NEQN,NPDES CALL F_I(I,NPDES,T,Y(I:I+NPDES-1),TEMP,.FALSE.,JAC) NFI = NFI + 1 YP(I:I+NPDES-1) = YP(I:I+NPDES-1) + TEMP END DO ELSE JACNRM = 0D0 DO I = 1,NEQN,NPDES CALL F_I(I,NPDES,T,Y(I:I+NPDES-1),TEMP,.TRUE.,JAC) NFI = NFI + 1 YP(I:I+NPDES-1) = YP(I:I+NPDES-1) + TEMP ! Compute maximum norm of the matrix JAC. JACNRM = MAX(JACNRM,MAXVAL( SUM(ABS(JAC),DIM=2) )) END DO END IF END SUBROUTINE ADD_F_I SUBROUTINE STEP(NEQN,F_E,NPDES,F_I,T,YN,H,M,Y,YJM1,YJM2,& IT_CONVG_FAIL,ERR_FLAG) ! Take an implicit step of size H from T to T+H ! to get Y(*). INTEGER :: NEQN,M,NPDES,ERR_FLAG DOUBLE PRECISION :: T,H DOUBLE PRECISION, DIMENSION(NEQN) :: YN,Y,YJM1,YJM2 LOGICAL :: IT_CONVG_FAIL EXTERNAL F_E,F_I INTEGER :: J,I DOUBLE PRECISION :: AJM1,ARG,BJ,BJM1,BJM2,DZJ,DZJM1,& DZJM2,D2ZJ,D2ZJM1,D2ZJM2,MU,MUS,NU,TEMP1,TEMP2,THJ,& THJM1,THJM2,W0,W1,ZJ,ZJM1,ZJM2,GAMMA,MUS1 DOUBLE PRECISION :: RHS(NEQN),DYE_0(NEQN),DYE_JM1(NEQN),& DYI_0(NEQN),DYI(NEQN,3),JAC(NPDES,NPDES) ! Swap columns of DYI by swapping these pointers: INTEGER :: PY,PY_JM1,PY_JM2,ITEMP PY = 1 PY_JM1 = 2 PY_JM2 = 3 ! Flag to report convergence failure. IT_CONVG_FAIL = .FALSE. W0 = 1D0 + 2D0/(13D0*M**2) TEMP1 = W0**2 - 1D0 TEMP2 = SQRT(TEMP1) ARG = M*LOG(W0 + TEMP2) W1 = SINH(ARG)*TEMP1 / (COSH(ARG)*M*TEMP2 - W0*SINH(ARG)) BJM1 = 1D0/W0 BJM2 = 1D0/(2D0*W0)**2 ! Evaluate the first stage. ! FN of explicit step is recomputed because ! in calling program it includes the implicit ! part which we need here in its own right. CALL F_E(NEQN,T,YN,DYE_0) NFE = NFE + 1 DO I = 1,NEQN,NPDES CALL F_I(I,NPDES,T,YN(I:I+NPDES-1),& DYI_0(I:I+NPDES-1),.FALSE.,JAC) NFI = NFI + 1 END DO DYI(:,PY_JM2) = DYI_0 ! Initialize for j = 2 YJM2 = YN MUS = W1*BJM1 MUS1 = MUS HMUS1 = H*MUS1 ! Global variable RHS = YN + HMUS1*DYE_0 YJM1 = YN CALL STAGE(NEQN,NPDES,F_I,T+HMUS1,RHS,YJM1,& DYI(:,PY_JM1),IT_CONVG_FAIL,ERR_FLAG) IF (IT_CONVG_FAIL .OR. ERR_FLAG > 0) RETURN THJM2 = 0D0 THJM1 = MUS ZJM1 = W0 ZJM2 = 1D0 DZJM1 = 1D0 DZJM2 = 0D0 D2ZJM1 = 0D0 D2ZJM2 = 0D0 ! Evaluate stages j = 2,...,m. DO J = 2, M ZJ = 2D0*W0*ZJM1 - ZJM2 DZJ = 2D0*W0*DZJM1 - DZJM2 + 2D0*ZJM1 D2ZJ = 2D0*W0*D2ZJM1 - D2ZJM2 + 4D0*DZJM1 BJ = D2ZJ/DZJ**2 AJM1 = 1D0 - ZJM1*BJM1 MU = 2D0*W0*BJ/BJM1 NU = - BJ/BJM2 MUS = MU*W1/W0 GAMMA = - AJM1*MUS CALL F_E(NEQN,T + H*THJM1,YJM1,DYE_JM1) NFE = NFE + 1 RHS = (1D0 - MU - NU)*YN + MU*YJM1 + NU*YJM2 + & H*MUS*DYE_JM1 + H*GAMMA*DYE_0 + & (GAMMA - (1D0 - MU - NU)*MUS1)*H*DYI_0 & - NU*HMUS1*DYI(:,PY_JM2) Y = YJM1 ! Guess for next stage THJ = MU*THJM1 + NU*THJM2 + MUS*(1D0 - AJM1) CALL STAGE(NEQN,NPDES,F_I,T+H*THJ,& RHS,Y,DYI(:,PY),IT_CONVG_FAIL,ERR_FLAG) IF (IT_CONVG_FAIL .OR. ERR_FLAG > 0) RETURN ! Shift the data for the next stage. IF (J < M) THEN YJM2 = YJM1 YJM1 = Y ! Swap column pointers. ITEMP = PY_JM2 PY_JM2 = PY_JM1 PY_JM1 = PY PY = ITEMP THJM2 = THJM1 THJM1 = THJ BJM2 = BJM1 BJM1 = BJ ZJM2 = ZJM1 ZJM1 = ZJ DZJM2 = DZJM1 DZJM1 = DZJ D2ZJM2 = D2ZJM1 D2ZJM1 = D2ZJ END IF END DO END SUBROUTINE STEP SUBROUTINE STAGE(NEQN,NPDES,F_I,TNJ,RHS,YS,DYS,& IT_CONVG_FAIL,ERR_FLAG) INTEGER :: NEQN,NPDES,ERR_FLAG DOUBLE PRECISION :: TNJ,RHS(NEQN),YS(NEQN),DYS(NEQN) LOGICAL :: IT_CONVG_FAIL EXTERNAL F_I INTEGER, PARAMETER :: ITMAX=10 INTEGER :: I,M,INFO,IPVT(NPDES) DOUBLE PRECISION :: ERR,W(NPDES),RES(NPDES),WT(NPDES),& ITMAT(NPDES,NPDES) IT_CONVG_FAIL = .FALSE. ! Compute implicit stage a block at at time DO I = 1,NEQN,NPDES ! Initialize guess for block W = YS(I:I+NPDES-1) ! Form and factor iteration matrix. CALL F_I(I,NPDES,TNJ,W,DYS(I:I+NPDES-1),.TRUE.,ITMAT) NFI = NFI + 1 ! HMUS1 = H*MUS1 is a global variable. ITMAT = EYE(NPDES) - HMUS1*ITMAT CALL FACTOR(ITMAT,NPDES,INFO,IPVT) IF (INFO > 0) THEN ERR_FLAG = 5 RETURN END IF DO M = 1,ITMAX ! Compute the residual of current iterate W. RES = RHS(I:I+NPDES-1) + HMUS1*DYS(I:I+NPDES-1) - W ! Compute the correction. RES = SOLVE(ITMAT,NPDES,IPVT,RES) WT = ABS(W) W = W + RES WT = MAX(WT,ABS(W)) ! Prepare for another iteration; ! place result for slope in DYS. CALL F_I(I,NPDES,TNJ,W,DYS(I:I+NPDES-1),.FALSE.,ITMAT) NFI = NFI + 1 ! Test for convergence ERR = RMSNORM( RES/(ATOL + RTOL*WT) ) IF (ERR <= 0.5D0) THEN EXIT ELSEIF (M == ITMAX) THEN IT_CONVG_FAIL = .TRUE. RETURN END IF END DO ! Place result for solution in YS. YS(I:I+NPDES-1) = W END DO END SUBROUTINE STAGE SUBROUTINE ERR_EST(NEQN,NPDES,T,H,YNEW,YP,Y,YLAST,& F_I,ERR,ERR_FLAG) ! Here YNEW is y(T) and Y is the tentative y(T+H). Reuse of ! storage has y'(T) stored in YP and y'(T+H) is to be returned ! in YLAST. INTEGER :: NEQN,NPDES,ERR_FLAG,I,INFO DOUBLE PRECISION :: T,H,ERR DOUBLE PRECISION, DIMENSION(NEQN) :: YNEW,YP,Y,YLAST EXTERNAL F_I DOUBLE PRECISION :: TEMP(NPDES),EST(NPDES),FMAT(NPDES,NPDES) INTEGER :: IPVT(NPDES) ! The error estimate is formed, filtered, weighted, and ! the squares of the components are summed in blocks. ! After all blocks are processed, the RMS norm is formed. ! Form y'(T+H) in YLAST. The input YLAST contains F_E. ! Here F_I is evaluated in blocks and added. The input ! YP holds y'(T). However, we also need F_I at y(T), ! so it is evaluated in blocks and used in the error ! estimate in blocks. ERR = 0D0 DO I = 1,NEQN,NPDES CALL F_I(I,NPDES,T+H,Y(I:I+NPDES-1),TEMP,.FALSE.,FMAT) NFI = NFI + 1 ! Form a block of y'(T+H) in YLAST. YLAST(I:I+NPDES-1) = YLAST(I:I+NPDES-1) + TEMP EST = 0.5D0*H*(YLAST(I:I+NPDES-1) - YP(I:I+NPDES-1)) & + HMUS1*TEMP ! Get Jacobian along with F_I. CALL F_I(I,NPDES,T,YNEW(I:I+NPDES-1),TEMP,.TRUE.,FMAT) NFI = NFI + 1 ! Complete computation of the unfiltered estimate. EST = EST - HMUS1*TEMP ! Form and factor the filter. FMAT = EYE(NPDES) - H*FMAT CALL FACTOR(FMAT,NPDES,INFO,IPVT) IF (INFO > 0) THEN ERR_FLAG = 6 RETURN END IF ! Filter the estimate. EST = SOLVE(FMAT,NPDES,IPVT,EST) ! Weight the components of the estimated error. TEMP = EST / ( ATOL + & RTOL*MAX(ABS(Y(I:I+NPDES-1)),ABS(YNEW(I:I+NPDES-1))) ) ! Sum the squares of the weighted error estimates. ERR = ERR + DOT_PRODUCT(TEMP,TEMP) END DO ! Form the RMS norm ERR = SQRT(ERR/NEQN) END SUBROUTINE ERR_EST SUBROUTINE RHO(F_E,NEQN,T,YN,FN,V,FV,EV,SPCRAD,RHO_CONVG_FAIL) ! RHO attempts to compute a close upper bound, SPCRAD, on ! the spectral radius of the Jacobian matrix using a nonlinear ! power method. A convergence failure is reported by ! RHO_CONVG_FAIL = .TRUE. EXTERNAL F_E INTEGER :: NEQN DOUBLE PRECISION :: T,SPCRAD DOUBLE PRECISION, DIMENSION(NEQN) :: YN,FN,V,FV,EV LOGICAL :: RHO_CONVG_FAIL INTEGER, PARAMETER :: ITMAX=50 INTEGER :: ITER,INDEX DOUBLE PRECISION :: YNRM,SIGMA,SIGMAL,DYNRM,DFNRM,VNRM,SMALL ! HMAX is a global variable. SPCRAD smaller than 1/HMAX are ! not interesting because they do not constrain the step size. SMALL = 1D0/HMAX RHO_CONVG_FAIL = .FALSE. ! The initial slope is used as guess when NSTEPS = 0 and ! thereafter the last computed eigenvector. Some care ! is needed to deal with special cases. Approximations to ! the eigenvector are normalized so that their Euclidean ! norm has the constant value DYNRM. ! It is the spectral radius of the Jacobian of F_E that is ! computed, but the input value FN includes the effect of ! the implicit term, so FN must be recomputed here. CALL F_E(NEQN,T,YN,FN) NFESIG = NFESIG + 1 YNRM = ENORM(YN) IF (NSTEPS == 0) THEN V = FN ELSE V = EV END IF VNRM = ENORM(V) IF (ABS(YNRM) > 0D0 .AND. ABS(VNRM) > 0D0) THEN DYNRM = YNRM*SQRTU V = YN + V*(DYNRM/VNRM) ELSEIF (ABS(YNRM) > 0D0) THEN DYNRM = YNRM*SQRTU V = YN*(1D0 + SQRTU) ELSEIF (ABS(VNRM) > 0D0) THEN DYNRM = UROUND V = V*(DYNRM/VNRM) ELSE DYNRM = UROUND V = DYNRM END IF ! Now iterate with a nonlinear power method. SIGMA = 0D0 DO ITER = 1,ITMAX CALL F_E(NEQN,T,V,FV) NFESIG = NFESIG + 1 DFNRM = ENORM(FV - FN) SIGMAL = SIGMA SIGMA = DFNRM/DYNRM ! SPCRAD is a little bigger than the estimate SIGMA of the ! spectral radius, so is more likely to be an upper bound. SPCRAD = 1.2D0*SIGMA IF (ITER >= 2 .AND. & ABS(SIGMA - SIGMAL) <= MAX(SIGMA,SMALL)*0.01D0) THEN EV = V - YN RETURN END IF ! The next V(*) is the change in f ! scaled so that norm(V - YN) = DYNRM. IF (ABS(DFNRM) > 0D0) THEN V = YN + (FV - FN)*(DYNRM/DFNRM) ELSE ! The new V(*) degenerated to YN(*)--"randomly" perturb ! current approximation to the eigenvector by changing ! the sign of one component. INDEX = 1 + MOD(ITER,NEQN) V(INDEX) = YN(INDEX) - (V(INDEX) - YN(INDEX)) END IF END DO ! Set flag to report a convergence failure. RHO_CONVG_FAIL = .TRUE. END SUBROUTINE RHO SUBROUTINE REPORT_ERRORS(ERR_FLAG) INTEGER :: ERR_FLAG SELECT CASE (ERR_FLAG) CASE (1) PRINT *,' Must have 10*UROUND <= RE <= 0.1.' CASE (2) PRINT *,' AE must be positive.' CASE (3) PRINT *,' Unable to achieve the desired accuracy with the',& ' precision available.' CASE (4) PRINT *,' Convergence failure in estimating the spectral radius.' CASE (5) PRINT *,' Singular iteration matrix.' CASE (6) PRINT *,' Singular filter matrix.' CASE (7) PRINT *,' A storage allocation error occurred in IRKC.' END SELECT STOP END SUBROUTINE REPORT_ERRORS FUNCTION ENORM(V) ! Euclidean norm of a vector. DOUBLE PRECISION :: V(:),ENORM ENORM = SQRT( DOT_PRODUCT(V,V) ) END FUNCTION ENORM FUNCTION RMSNORM(V) ! RMS norm of a vector. DOUBLE PRECISION :: V(:),RMSNORM RMSNORM = SQRT( DOT_PRODUCT(V,V) / SIZE(V) ) END FUNCTION RMSNORM FUNCTION EYE(N) RESULT(IDENTITY) INTEGER :: N,I,J DOUBLE PRECISION :: IDENTITY(N,N) IDENTITY = RESHAPE( (/ ((0D0,J=1,I-1),1D0,(0D0,J=I+1,N), & I=1,N) /), (/ N,N /) ) END FUNCTION EYE SUBROUTINE FACTOR(A,NEQ,FLAG,PIVOTS) ! FACTOR decomposes the matrix A using Gaussian elimination. It ! is used in conjunction with SOLVE to solve A*X = B, a system ! of NEQ linear equations in NEQ unknowns. It is sometimes ! convenient for the matrix to be defined as the portion ! A(1:NEQ,1:NEQ) of a larger matrix, and similarly for the ! vector B. FACTOR/SOLVE are coded so that this is possible. ! ! Input arguments: ! A = real matrix A(1:NEQ,1:NEQ) to be factored. ! NEQ = the number of equations to be solved, an integer. ! ! Output arguments: ! A = contains the upper triangular matrix U in its upper ! portion (by rows) and a permuted version of a lower ! triangular matrix (I - L). The factorization is ! such that (permutation matrix)*A = L*U. ! FLAG = an integer that reports success or the reason for ! failure. FLAG = 0 indicates success. If FLAG > 0, ! a zero pivot occurred at equation number FLAG and ! the computation was terminated. FLAG = -1 means ! there was an input error: either NEQ <= 0 or A does ! not have at least NEQ rows and columns. ! PIVOTS = an integer vector of at least NEQ entries that ! records row interchanges. The entry PIVOTS(NEQ) = ! (-1)**(number of row interchanges). ! ! When FLAG > 0, the determinant of A is 0 and when FLAG = 0, ! det(A) = pivots(NEQ) * A(1,1) * ... * A(NEQ,NEQ). INTEGER :: NEQ, FLAG, PIVOTS(NEQ) DOUBLE PRECISION :: A(:,:) INTEGER :: I, K, M, OCCURRED(1) DOUBLE PRECISION :: T, TEMP(NEQ) IF ((NEQ <= 0) .OR. & (SIZE(A,DIM=1) < NEQ) .OR. (SIZE(A,DIM=2) < NEQ)) THEN FLAG = -1 RETURN END IF FLAG = 0 PIVOTS(NEQ) = 1 IF (NEQ == 1) THEN IF (ABS(A(1,1)) <= 0D0) THEN FLAG = 1 RETURN END IF END IF ! Gaussian elimination with partial pivoting. DO K = 1,NEQ-1 ! Determine the row M containing the largest element in ! magnitude to be used as a pivot. OCCURRED = MAXLOC(ABS(A(K:NEQ,K))) M = OCCURRED(1) + K - 1 ! If all possible pivots are 0, A is numerically singular. IF (ABS(A(M,K)) <= 0D0) THEN FLAG = K RETURN END IF PIVOTS(K) = M IF (M /= K) THEN ! Interchange the current row K with the pivot row M. TEMP(K:NEQ) = A(K,K:NEQ) A(K,K:NEQ) = A(M,K:NEQ) A(M,K:NEQ) = TEMP(K:NEQ) PIVOTS(NEQ) = - PIVOTS(NEQ) END IF ! Eliminate subdiagonal entries of column K. DO I = K+1,NEQ T = A(I,K)/A(K,K) A(I,K) = - T IF (ABS(T) > 0D0) THEN A(I,K+1:NEQ) = A(I,K+1:NEQ) - T*A(K,K+1:NEQ) END IF END DO END DO IF (ABS(A(NEQ,NEQ)) <= 0D0) THEN FLAG = NEQ RETURN END IF END SUBROUTINE FACTOR FUNCTION SOLVE(A,NEQ,PIVOTS,B) RESULT(X) ! SOLVE solves A*X = B, a system of NEQ linear equations in NEQ ! unknowns using the decomposition obtained from a successful call ! to FACTOR. ! ! Input arguments: ! A = output of FACTOR. Triangular decomposition ! of the coefficient matrix. ! NEQ = the number of equations to be solved. ! PIVOTS = output of FACTOR. Record of row interchanges. ! B = right hand side vector B, a real vector of at ! least NEQ entries. ! ! Output argument: ! X = solution vector of the same size and type as B. INTEGER :: NEQ,PIVOTS(:) DOUBLE PRECISION :: A(:,:),B(:) DOUBLE PRECISION, DIMENSION(SIZE(B)) :: X INTEGER :: I, K, M DOUBLE PRECISION :: TEMP X = B IF (NEQ == 1) THEN X(1) = X(1)/A(1,1) ELSE ! Forward elimination. DO K = 1,NEQ-1 M = PIVOTS(K) TEMP = X(M) X(M) = X(K) X(K) = TEMP X(K+1:NEQ) = X(K+1:NEQ) + A(K+1:NEQ,K)*X(K) END DO ! Back substitution. X(NEQ) = X(NEQ)/A(NEQ,NEQ) DO I = NEQ-1,1,-1 X(I) = (X(I) - DOT_PRODUCT(A(I,I+1:NEQ),X(I+1:NEQ)))/A(I,I) END DO END IF END FUNCTION SOLVE END MODULE IRKC_M