C Work performed under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number 
C W-7405-Eng-48.
C
C Copyright 1995 the Regents of the University of California.
C All rights reserved.
C
C-----------------------------------------------------------------------
C
C    Preconditioner Tools for Reaction-Transport Problems
C    Part I: Block-Diagonal Reaction-Based Factor without Grouping
C                        14 September 1995
C
C The following three subroutines -- DMSET2, DRBDJA, DRBDPS --
C are provided to assist in the generation and solution of
C preconditioner matrices for problems arising from reaction-transport
C systems, as solved with DASPK.  More specifically, they are intended
C as tools for preconditioners that include a contribution from the
C reaction terms of the system.  These are intended as auxiliary
C routines for the user-supplied routines JAC and PSOL called by
C DASPK when the Krylov method is selected.
C
C These routines are intended for a DAE system obtained from a system
C of reaction-transport PDEs, in which some of the PDE variables obey
C evolution equations, and the rest obey algebraic (time-independent)
C equations.  See Ref. 2, Section 4.  It is assumed that the 
C right-hand sides of all the equations have the form of a sum of a
C reaction term R and a transport term S, that the transport term
C is discretized by finite differences, and that in the spatial
C discretization the PDE variables at each spatial point are kept
C together.  Thus the DAE system function, in terms of a dependent
C variable vector u, has the form
C     G(t,u,u') = I_d u' - R(t,u) - S(t,u) ,  where
C     I_d = identity matrix with zeros in the positions corresponding
C           to the algebraic components, ones in those for the
C           evolution (differential) components.
C     R(t,u) = the reaction terms (spatial coupling absent), and
C     S(t,u) = the spatial transport terms.
C
C As shown in [2], two possible preconditioners for such a system are:
C (a) P_R = c I_d - dR/du, based on the reaction term R alone, and
C (b) P_SR = (I - (1/c) dS/du) (c I_d - dR/du), the product of two
C     factors (in either order), one being P_R and the other being
C     based on the spatial term S alone.
C Here c is the scalar CJ that is input to the JAC and PSOL routines
C provided by the user (1/c is proportional to the step size H).
C
C The routines given here can be used to provide the reaction-based
C factor P_R.  More precisely, they provide an approximation A_R to 
C P_R.  The matrix P_R is block-diagonal, with each block corresponding
C to one spatial point.  In A_R, we compute each block by difference
C quotient approximations, by way of calls to a user-supplied routine,
C subroutine RBLOCK, that evaluates the reaction terms at a single
C spatial point.  A_R has one such block for each spatial point in
C the mesh.  (For a more economical approximation, see Part II,
C on block-grouping in A_R.)
C
C The routines given here are specialized to the case of a 2-D problem
C on a rectangular mesh in the x-y plane.  However, they can be easily
C modified for a different problem geometry.  It is also assumed
C that the PDE variables are ordered so that the differential
C variables appear first, followed by the algebraic variables.
C
C To make use of these routines in a DASPK solution, the user must
C provide:
C (a) a calling program that sets the DASPK input parameters, and calls
C     DMSET2 to set mesh data and mesh-related DASPK inputs;
C (b) a JAC routine, as prescribed by the DASPK instructions, which
C     calls DRBDJA, and does any other Jacobian-related preprocessing
C     needed for preconditioning; and
C (c) a PSOL routine, as prescribed by the DASPK instructions, which
C     calls DRBDPS for the solution of systems A_R x = b, and does
C     any other linear system solving required by the preconditioner.
C Detailed descriptions and instructions are given below.
C
C In addition, the use of these routines requires:
C  * the LINPACK routines DGEFA and DGESL for dense linear sytems, and
C  * the machine constant routine D1MACH for the machine unit roundoff.
C
C (a) The calling program.
C The calling program sets the DASPK inputs and makes calls to DASPK.
C Here the DASPK inputs include
C   INFO(12) = 1 [to signal the Krylov method]
C   INFO(15) = 1 [to signal the presence of a JAC routine]
C
C Also, the use of the DRBDJA/DRBDPS routines in conjunction with
C DASPK requires that certain mesh-related data be set.  This can be
C done with the call
C     CALL DMSET2 (MX, MY, NS, NSD, LID, IWORK)
C The input arguments to DMSET2 are:
C   MX and MY = the mesh dimensions.
C   NS  = number of PDE variables.
C   NSD = number of differential PDE variables.
C   LID = offset in IWORK for array showing the differential and
C         algebraic components on input to DASPK, required if either
C         INFO(11) = 1 or INFO(16) = 1.  Set LID = 0 otherwise.
C         If this array is required, set LID = 40 or 40 + NEQ,
C         depending on the value of the constraint option INFO(10).
C DMSET2 loads mesh data in a COMMON block /DRPRE1/ used by the
C DRBDJA/DRBDPS routines.
C
C DMSET2 also loads the preconditioner work lengths into
C IWORK(27) and IWORK(28), and if LID > 0 it sets the ID array
C in IWORK showing the differential and algebraic components.
C 
C (b) The JAC routine.
C The user-supplied JAC routine called by DASPK with the Krylov
C method specified, is to generate and preprocess Jacobian-related
C data as needed for later solution of the preconditioner system
C P x = b.  Assuming that P is to be an approximation of either P_R
C or P_SR, the JAC routine should call DRBDJA for the approximation
C A_R to P_R.  Subroutine DRBDJA generates A_R using difference
C quotients.  It then performs an LU decomposition of each block,
C using the LINPACK routine DGEFA.
C
C In terms of the arguments passed to JAC by DASPK, the call to
C DRBDJA should have the form
C     CALL DRBDJA (T, U, R0, RBLOCK, WK, REWT, CJ, WP, IWP, IER)
C where we use U instead of Y for the dependent variable array.
C The argument R0 is an array assumed to contain the current value
C of the R vector, at the current values (T,U).  This can be done, for
C example, by taking R0 to be RPAR, and loading RPAR with the
C vector R in the last call to the RES routine; in that case, the
C calling program must declare RPAR to have length at least NEQ.
C Alternatively, insert a call to RBLOCK (see below) within the
C loop over mesh points in DRBDJA.
C
C To use DRBDJA, the user must provide the following subroutine,
C which DRBDJA calls to obtain individual blocks of R:
C      SUBROUTINE RBLOCK (T, JX, JY, UXY, RXY, senpara)
C The input arguments to RBLOCK are:
C   T     = current time.
C   JX,JY = spatial indices in x- and y-directions.
C   UXY   = block of NS dependent variables at spatial point (JX,JY).
C RBLOCK is to load block (JX,JY) of R(t,u) into the array RXY.
C
C (c) The PSOL routine.
C The user-supplied PSOL routine must solve the linear system P x = b,
C where P is the preconditioner matrix.  For this, the PSOL routine
C should call DRBDPS for the solution of A_R.  Subroutine DRBDPS
C solves a linear system A_R x = b, using the LINPACK backsolve
C routine DGESL.  In terms of the arguments passed to PSOL by DASPK,
C the call to DRBDPS should have the form
C     CALL DRBDPS (B, WP, IWP)
C DRBDPS overwrites the B array (containing b) with the solution x.
C
C-----------------------------------------------------------------------
C
C References
C [1] Peter N. Brown and Alan C. Hindmarsh,
C     Reduced Storage Matrix Methods in Stiff ODE Systems,
C     J. Appl. Math. & Comp., 31 (1989), pp. 40-91.
C [2] Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold,
C     Using Krylov Methods in the Solution of Large-Scale Differential-
C     Algebraic Systems, SIAM J. Sci. Comput., 15 (1994), pp. 1467-1488.
C [3] Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold,
C     Consistent Initial Condition Calculation for Differential-
C     Algebraic Systems, LLNL Report UCRL-JC-122175, August 1995;
C     submitted to SIAM J. Sci. Comp.
C-----------------------------------------------------------------------


      SUBROUTINE DMSET2 (MX, MY, NS, NSD, LID, IWORK)
C***BEGIN PROLOGUE  DMSET2
C***DATE WRITTEN   950830   (YYMMDD)
C
C***AUTHORS  A. C. Hindmarsh
C            Lawrence Livermore National Laboratory
C            L-316, P.O. Box 808
C            Livermore, CA 94551
C
C***DESCRIPTION
C
C-----------------------------------------------------------------------
C This routine sets mesh parameters needed to use the routines
C DRBDJA and DRBDPS, assuming a 2-D rectangular problem.
C Given the mesh parameters, it loads the COMMON block /DRPRE1/,
C and the lengths LENWP and LENIWP in IWORK.
C Then if LID > 0, it also sets the ID array in IWORK, indicating
C which components are differential and which are algebraic.
C
C The variables in the COMMON block are defined as follows:
C   SRUR   = SQRT(unit roundoff), used in difference quotients.
C            UROUND = D1MACH(4) generates the unit roundoff.
C   MP     = NS = number of PDE variables, the size of each block in 
C            the block-diagonal preconditioner matrix P_R.
C   MPD    = NSD = number of differential PDE variables.  In the DAE
C            system, the first MPD variables at each spatial point have
C            time  derivatives, and the remaining (MP - MPD) do not.
C   MPSQ   = MP*MP.
C   MESHX  = MX = x mesh size.
C   MESHY  = MY = y mesh size (the mesh is MESHX by MESHY).
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   D1MACH
C
C***END PROLOGUE  DMSET2
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION IWORK(*)
      COMMON /DRPRE1/ SRUR, MP, MPD, MPSQ, MESHX, MESHY, MXMP
C
C Load the COMMON block.
      UROUND = D1MACH(4)
      SRUR = SQRT(UROUND)
      MP = NS
      MPD = NSD
      MPSQ = NS*NS
      MESHX = MX
      MESHY = MY
      MXMP = MESHX*MP
C
C Here set the sizes of the preconditioning storage space segments
C in RWORK and IWORK.
      IWORK(27) = MPSQ*MESHX*MESHY
      IWORK(28) = MP*MESHX*MESHY
C
C If LID .GT. 0, set the ID array in IWORK.
      IF (LID .EQ. 0) RETURN
      I0 = LID
      DO 40 JY = 1,MY
        DO 30 JX = 1,MX
          DO 10 I = 1,MPD
 10         IWORK(I0+I) = 1
          DO 20 I = MPD+1,MP
 20         IWORK(I0+I) = -1
          I0 = I0 + MP
 30       CONTINUE
 40     CONTINUE
C
      RETURN
C------------  End of Subroutine DMSET2  -------------------------------
      END

      SUBROUTINE DRBDJA (
     *     T, U, R0, RBLOCK, R1, REWT, CJ, BD, IPBD, IER, senpara)
C***BEGIN PROLOGUE  DRBDJA
C***DATE WRITTEN   950914   (YYMMDD)
C
C***AUTHORS  A. C. Hindmarsh
C            Lawrence Livermore National Laboratory
C            L-316, P.O. Box 808
C            Livermore, CA 94551
C
C***DESCRIPTION
C
C-----------------------------------------------------------------------
C This routine generates and preprocesses a block-diagonal
C preconditioner matrix, based on the part of the Jacobian corresponding
C to the reaction terms R of the problem.
C It generates a matrix of the form CJ * I_d - dR/du.
C It calls DGEFA to do LU decomposition of each diagonal block.
C The computation of the diagonal blocks uses the mesh information in
C the COMMON block /DRPRE1/.  One block per spatial point is computed.
C The Jacobian elements are generated by difference quotients.
C This routine calls a user routine of the form
C      SUBROUTINE RBLOCK (T, JX, JY, UXY, RXY, senpara)
C which is to set RXY to block (JX,JY) of R, as a function of the
C current time T and block UXY of current dependent variable vector U.
C The array R0 is assumed to contain the current value of R at (T,U).
C-----------------------------------------------------------------------
C On input:
C   T      = current value of independent variable.
C   U      = current dependent variable array.
C   R0  = array of current values of the vector R at (T,U)
C   RBLOCK = name of external routine that computes a single block of R.
C   R1     = array of length NEQ for work space.
C   REWT   = reciprocal error weights.
C   CJ     = scalar used in forming the system Jacobian.
C
C On output:
C   BD     = array containing the LU factors of the diagonal blocks.
C   IPBD   = integer array of pivots for the LU factorizations.
C   IER    = integer error flag.  If no error occurred, IER = 0.
C            If a zero pivot was found at stage k in one of the LU
C            factorizations, this routine returns IER = k > 0.
C Here BD is the RWORK segment WP, and IPBD is the IWORK segment IWP.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   RBLOCK, DGEFA
C
C***END PROLOGUE  DRBDJA
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL RBLOCK
      DIMENSION U(*), R0(*), R1(*), REWT(*), BD(*),IPBD(*),senpara(*)
     *     
C
      COMMON /DRPRE1/ SRUR, MP, MPD, MPSQ, MESHX, MESHY, MXMP
C
C Make MP calls to RBLOCK to approximate each diagonal block of dR/du. 
      DFAC = 1.0D-2
      IBD = 0
      J0 = 0
      DO 40 JY = 1,MESHY
        DO 30 JX = 1,MESHX
C If R0 has not been set previously as an array of length NEQ, it can
C be set here, as an array of length MP, with the call
C         CALL RBLOCK (T, JX, JY, U(J0+1), R0)
C In this case, change R0(J0+I) below to R0(I).
          DO 20 JS = 1,MP
            J = J0 + JS
            UJ = U(J)
            DEL = MAX(SRUR*ABS(UJ),DFAC/REWT(J))
            U(J) = U(J) + DEL
            FAC = -1.0D0/DEL
            CALL RBLOCK (T, JX, JY, U(J0+1), R1, senpara)
            DO 10 I = 1,MP
 10           BD(IBD+I) = (R1(I) - R0(J0+I))*FAC
            U(J) = UJ
            IBD = IBD + MP
 20         CONTINUE
          J0 = J0 + MP
 30       CONTINUE
 40     CONTINUE
C
C Add matrix CJ * I_d, and do LU decomposition on blocks. --------------
      IBD = 1
      IIP = 1
      DO 80 J = 1,MESHX*MESHY
        IDIAG = IBD
        DO 70 I = 1,MP
          IF (I .LE. MPD) BD(IDIAG) = BD(IDIAG) + CJ
 70       IDIAG = IDIAG + (MP + 1)
        CALL DGEFA (BD(IBD), MP, MP, IPBD(IIP), IER)
        IF (IER .NE. 0) GO TO 90
        IBD = IBD + MPSQ
        IIP = IIP + MP
 80     CONTINUE
 90   RETURN
C------------  End of Subroutine DRBDJA  -------------------------------
      END

      SUBROUTINE  DRBDPS (B, BD, IPBD)
C***BEGIN PROLOGUE  DRBDPS
C***DATE WRITTEN   950914   (YYMMDD)
C
C***AUTHORS  A. C. Hindmarsh
C            Lawrence Livermore National Laboratory
C            L-316, P.O. Box 808
C            Livermore, CA 94551
C
C***DESCRIPTION
C
C-----------------------------------------------------------------------
C This routine solves a linear system A_R x = b, using the LU factors
C of the diagonal blocks computed in DRBDJA, and mesh parameters
C in the COMMON block /DRPRE1/.
C Here BD is the RWORK segment WP, and IPBD is the IWORK segment IWP.
C The right-hand side vector b, contained in B on entry, is overwritten
C with the solution vector x on return.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   DGESL
C
C***END PROLOGUE  DRBDPS
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION B(*), BD(*), IPBD(*)
C
      COMMON /DRPRE1/ SRUR, MP, MPD, MPSQ, MESHX, MESHY, MXMP
C
      IER = 0
      IB = 1
      IBD = 1
      DO 20 JY = 1,MESHY
        DO 10 JX = 1,MESHX
          CALL DGESL (BD(IBD), MP, MP, IPBD(IB), B(IB), 0)
          IB = IB + MP
          IBD = IBD + MPSQ
 10       CONTINUE
 20     CONTINUE
C
      RETURN
C------------  End of Subroutine DRBDPS  -------------------------------
      END
