irkc_m.f90

!------------------------------------------------------------------
!
! 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,SPRAD
    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
      ! print *,'SPRAD = ',SPRAD

      !  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
    SOL%SPRAD  = SPRAD

    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