text:
*DECK DIR
      SUBROUTINE DIR(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSOLVE,
     $     ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ,
     $     RWORK, IWORK)
C***BEGIN PROLOGUE  DIR
C***DATE WRITTEN   890404   (YYMMDD)
C***REVISION DATE  890404   (YYMMDD)
C***CATEGORY NO.  D2A4, D2B4
C***KEYWORDS  LIBRARY=SLATEC(SLAP),
C             TYPE=DOUBLE PRECISION(DIR-D),
C             Linear system, Sparse, Iterative Precondition
C***AUTHOR  Greenbaum, Anne, Courant Institute
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-300
C             Livermore, CA 94550 (415) 423-3141
C             seager@lll-crg.llnl.gov
C***PURPOSE  Preconditioned Iterative Refinement sparse Ax = b solver.
C            Routine to solve a general linear system  Ax = b  using
C            iterative refinement with a matrix splitting.
C***DESCRIPTION
C *Usage:
C     INTEGER  N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
C     INTEGER  ITER, IERR, IUNIT, IWORK(USER DEFINABLE)
C     DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N), DZ(N)
C     DOUBLE PRECISION RWORK(USER DEFINABLE)
C     EXTERNAL MATVEC, MSOLVE
C
C     CALL DIR(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSLOVE, ITOL,
C    $     TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK, IWORK)
C
C *Arguments:
C N      :IN       Integer.
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right-hand side vector.
C X      :INOUT    Double Precision X(N).
C         On input X is your initial guess for solution vector.
C         On output X is the final approximate solution.
C NELT   :IN       Integer.
C         Number of Non-Zeros stored in A.
C IA     :IN       Integer IA(NELT).
C JA     :IN       Integer JA(NELT).
C A      :IN       Integer A(NELT).
C         These arrays contain the matrix data structure for A.
C         It could take any form.  See "Description", below
C         for more late breaking details...
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all nonzero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the upper
C         or lower triangle of the matrix is stored.
C MATVEC :EXT      External.
C         Name of a routine which performs the matrix vector multiply
C         Y = A*X given A and X.  The name of the MATVEC routine must
C         be declared external in the calling program.  The calling
C         sequence to MATVEC is:
C             CALL MATVEC( N, X, Y, NELT, IA, JA, A, ISYM )
C         Where N is the number of unknowns, Y is the product A*X
C         upon return, X is an input vector, NELT is the number of
C         non-zeros in the SLAP IA, JA, A storage for the matrix A.
C         ISYM is a flag which, if non-zero, denotes that A is
C         symmetric and only the lower or upper triangle is stored.
C MSOLVE :EXT      External.
C         Name of a routine which solves a linear system MZ = R for
C         Z given R with the preconditioning matrix M (M is supplied via
C         RWORK and IWORK arrays).  The name of the MSOLVE routine must
C         be declared external in the calling program.  The calling
C         sequence to MSOLVE is:
C             CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
C         Where N is the number of unknowns, R is the right-hand side
C         vector, and Z is the solution upon return.  IA, JA, A and
C         ISYM are defined as above.  RWORK is a double precision array
C         that can be used to pass necessary preconditioning information
C         and/or workspace to MSOLVE.  IWORK is an integer work array
C         for the same purpose as RWORK.
C ITOL   :IN       Integer.
C         Flag to indicate type of convergence criterion.
C         If ITOL=1, iteration stops when the 2-norm of the residual
C         divided by the 2-norm of the right-hand side is less than TOL.
C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
C         residual divided by the 2-norm of M-inv times the right hand
C         side is less than TOL, where M-inv is the inverse of the
C         diagonal of A.
C         ITOL=11 is often useful for checking and comparing different
C         routines.  For this case, the user must supply the "exact"
C         solution or a very accurate approximation (one with an error
C         much less than TOL) through a common block,
C                     COMMON /SOLBLK/ SOLN(1)
C         if ITOL=11, iteration stops when the 2-norm of the difference
C         between the iterative approximation and the user-supplied
C         solution divided by the 2-norm of the user-supplied solution
C         is less than TOL.  Note that this requires the user to set up
C         the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine.
C         The routine with this declaration should be loaded before the
C         stop test so that the correct length is used by the loader.
C         This procedure is not standard Fortran and may not work
C         correctly on your system (although it has worked on every
C         system the authors have tried).  If ITOL is not 11 then this
C         common block is indeed standard Fortran.
C TOL    :IN       Double Precision.
C         Convergence criterion, as described above.
C ITMAX  :IN       Integer.
C         Maximum number of iterations.
C ITER   :OUT      Integer.
C         Number of iterations required to reach convergence, or
C         ITMAX+1 if convergence criterion could not be achieved in
C         ITMAX iterations.
C ERR    :OUT      Double Precision.
C         Error estimate of error in final approximate solution, as
C         defined by ITOL.
C IERR   :OUT      Integer.
C         Return error flag.
C           IERR = 0 => All went well.
C           IERR = 1 => Insufficient storage allocated
C                       for WORK or IWORK.
C           IERR = 2 => Method failed to converge in
C                       ITMAX steps.
C           IERR = 3 => Error in user input.  Check input
C                       value of N, ITOL.
C           IERR = 4 => User error tolerance set too tight.
C                       Reset to 500.0*D1MACH(3).  Iteration proceeded.
C           IERR = 5 => Preconditioning matrix, M,  is not
C                       Positive Definite.  $(r,z) < 0.0$.
C           IERR = 6 => Matrix A is not Positive Definite.
C                       $(p,Ap) < 0.0$.
C IUNIT  :IN       Integer.
C         Unit number on which to write the error at each iteration,
C         if this is desired for monitoring convergence.  If unit
C         number is 0, no writing will occur.
C R      :WORK     Double Precision R(N).
C Z      :WORK     Double Precision Z(N).
C DZ     :WORK     Double Precision DZ(N).
C RWORK  :WORK     Double Precision RWORK(USER DEFINABLE).
C         Double Precision array that can be used by  MSOLVE.
C IWORK  :WORK     Integer IWORK(USER DEFINABLE).
C         Integer array that can be used by  MSOLVE.
C
C *Description:
C       The basic algorithm for iterative refinement (also known as
C       iterative improvement) is:
C
C            n+1    n    -1       n
C           X    = X  + M  (B - AX  ).
C
C           -1   -1
C       If M =  A then   this is  the standard iterative  refinement
C       algorithm and the "subtraction" in  the residual calculation
C       should be done in double precision (which it is  not in this
C       routine).  If M = DIAG(A), the diagonal of A, then iterative
C       refinement is known  as Jacobi's  method.   The SLAP routine
C       DSJAC  implements this iterative strategy.   If  M = L,  the
C       lower  triangle of A,  then iterative refinement is known as
C       Gauss-Seidel.   The    SLAP  routine  DSGS  implements  this
C       iterative strategy.
C
C       This routine does  not care  what matrix data   structure is
C       used for  A and M.  It simply   calls  the MATVEC and MSOLVE
C       routines, with  the arguments as  described above.  The user
C       could write any type of structure and the appropriate MATVEC
C       and MSOLVE routines.  It is assumed  that A is stored in the
C       IA, JA, A  arrays in some fashion and  that M (or INV(M)) is
C       stored  in  IWORK  and  RWORK)  in  some fashion.   The SLAP
C       routines DSJAC and DSGS are examples of this procedure.
C
C       Two  examples  of  matrix  data structures  are the: 1) SLAP
C       Triad  format and 2) SLAP Column format.
C
C       =================== S L A P Triad format ===================
C
C       In  this   format only the  non-zeros are  stored.  They may
C       appear  in *ANY* order.   The user  supplies three arrays of
C       length NELT, where  NELT  is the number  of non-zeros in the
C       matrix:  (IA(NELT), JA(NELT),  A(NELT)).  For each  non-zero
C       the  user puts   the row  and  column index   of that matrix
C       element in the IA and JA arrays.  The  value of the non-zero
C       matrix  element is  placed in  the corresponding location of
C       the A  array.  This is  an extremely easy data  structure to
C       generate.  On  the other hand it  is  not too  efficient  on
C       vector  computers   for the  iterative  solution  of  linear
C       systems.  Hence, SLAP  changes this input  data structure to
C       the SLAP   Column  format for the  iteration (but   does not
C       change it back).
C
C       Here is an example of the  SLAP Triad   storage format for a
C       5x5 Matrix.  Recall that the entries may appear in any order.
C
C
C           5x5 Matrix       SLAP Triad format for 5x5 matrix on left.
C                              1  2  3  4  5  6  7  8  9 10 11
C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C       =================== S L A P Column format ==================
C       This routine requires  that the  matrix  A be  stored in the
C       SLAP Column format.  In this format the non-zeros are stored
C       counting down columns (except for  the diagonal entry, which
C       must appear first in each  "column")  and are stored  in the
C       double precision array A.   In other words,  for each column
C       in the matrix put the diagonal entry in  A.  Then put in the
C       other non-zero  elements going down  the column (except  the
C       diagonal) in order.   The  IA array holds the  row index for
C       each non-zero.  The JA array holds the offsets  into the IA,
C       A arrays  for  the  beginning  of each   column.   That  is,
C       IA(JA(ICOL)),  A(JA(ICOL)) points   to the beginning  of the
C       ICOL-th   column    in    IA and   A.      IA(JA(ICOL+1)-1),
C       A(JA(ICOL+1)-1) points to  the  end of the   ICOL-th column.
C       Note that we always have  JA(N+1) = NELT+1,  where N is  the
C       number of columns in  the matrix and NELT  is the number  of
C       non-zeros in the matrix.
C
C       Here is an example of the  SLAP Column  storage format for a
C       5x5 Matrix (in the A and IA arrays '|'  denotes the end of a
C       column):
C
C       5x5 Matrix      SLAP Column format for 5x5 matrix on left.
C       1  2  3    4  5    6  7    8    9 10 11
C       |11 12  0  0 15|   A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
C       |21 22  0  0  0|  IA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  JA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C *Examples:
C       See the SLAP routines DSJAC, DSGS
C
C *Precision:           Double Precision
C *See Also:
C       DSJAC, DSGS
C***REFERENCES  1. Gene Golub \& Charles Van Loan, "Matrix
C                 Computations", John Hopkins University Press; 3
C                 (1983) IBSN 0-8018-3010-9.
C***ROUTINES CALLED  MATVEC, MSOLVE, ISDIR.
C***END PROLOGUE  DIR
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER N, NELT, IA(NELT), JA(NELT), ISYM
      INTEGER ITOL, ITMAX, ITER, IERR, IUNIT, IWORK(*)
      DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N)
      DOUBLE PRECISION DZ(N), RWORK(*)
      EXTERNAL MSOLVE, MATVEC, ISDIR
C
C         Check some of the input data.
C***FIRST EXECUTABLE STATEMENT  DIR
      ITER = 0
      IERR = 0
      IF( N.LT.1 ) THEN
         IERR = 3
         RETURN
      ENDIF
      TOLMIN = 500.0*D1MACH(3)
      IF( TOL.LT.TOLMIN ) THEN
         TOL = TOLMIN
         IERR = 4
      ENDIF
C
C         Calculate initial residual and pseudo-residual, and check
C         stopping criterion.
      CALL MATVEC(N, X, R, NELT, IA, JA, A, ISYM)
      DO 10 I = 1, N
         R(I) = B(I) - R(I)
 10   CONTINUE
      CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
C
      IF( ISDIR(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL,
     $     ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK,
     $     IWORK, BNRM, SOLNRM) .NE. 0 ) GO TO 200
      IF( IERR.NE.0 ) RETURN
C
C         ***** iteration loop *****
C
      DO 100 K=1,ITMAX
         ITER = K
C
C         Calculate new iterate x, new residual r, and new
C         pseudo-resid z.
         DO 20 I = 1, N
            X(I) = X(I) + Z(I)
 20      CONTINUE
         CALL MATVEC(N, X, R, NELT, IA, JA, A, ISYM)
         DO 30 I = 1, N
            R(I) = B(I) - R(I)
 30      CONTINUE
         CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
C
C         check stopping criterion.
         IF( ISDIR(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL,
     $        ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK,
     $        IWORK, BNRM, SOLNRM) .NE. 0 ) GO TO 200
C
 100  CONTINUE
C
C         *****   end of loop  *****
C         Stopping criterion not satisfied.
      ITER = ITMAX + 1
      IERR = 2
C
 200  RETURN
C------------- LAST LINE OF DIR FOLLOWS -------------------------------
      END
*DECK DSJAC
      SUBROUTINE DSJAC(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
     $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW )
C***BEGIN PROLOGUE  DSJAC
C***DATE WRITTEN   890404   (YYMMDD)
C***REVISION DATE  890404   (YYMMDD)
C***CATEGORY NO.  D2A4, D2B4
C***KEYWORDS  LIBRARY=SLATEC(SLAP),
C             TYPE=DOUBLE PRECISION(DSJAC-D),
C             Linear system, Sparse, Iterative Precondition
C***AUTHOR  Greenbaum, Anne, Courant Institute
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-300
C             Livermore, CA 94550 (415) 423-3141
C             seager@lll-crg.llnl.gov
C***PURPOSE  Jacobi's method iterative sparse Ax = b solver.
C            Routine to solve a general linear system  Ax = b  using
C            Jacobi iteration.
C***DESCRIPTION
C *Usage:
C     INTEGER  N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
C     INTEGER  ITER, IERR, IUNIT, LENW, IWORK(LENIW), LENIW
C     DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(LENW)
C
C     CALL DSJAC(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
C    $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW )
C
C *Arguments:
C N      :IN       Integer.
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right-hand side vector.
C X      :INOUT    Double Precision X(N).
C         On input X is your initial guess for solution vector.
C         On output X is the final approximate solution.
C NELT   :IN       Integer.
C         Number of Non-Zeros stored in A.
C IA     :INOUT    Integer IA(NELT).
C JA     :INOUT    Integer JA(NELT).
C A      :INOUT    Integer A(NELT).
C         These arrays should hold the matrix A in either the SLAP
C         Triad format or the SLAP Column format.  See "Description",
C         below.  If the SLAP Triad format is chosen it is changed
C         internally to the SLAP Column format.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all nonzero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the upper
C         or lower triangle of the matrix is stored.
C ITOL   :IN       Integer.
C         Flag to indicate type of convergence criterion.
C         If ITOL=1, iteration stops when the 2-norm of the residual
C         divided by the 2-norm of the right-hand side is less than TOL.
C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
C         residual divided by the 2-norm of M-inv times the right hand
C         side is less than TOL, where M-inv is the inverse of the
C         diagonal of A.
C         ITOL=11 is often useful for checking and comparing different
C         routines.  For this case, the user must supply the "exact"
C         solution or a very accurate approximation (one with an error
C         much less than TOL) through a common block,
C                     COMMON /SOLBLK/ SOLN(1)
C         if ITOL=11, iteration stops when the 2-norm of the difference
C         between the iterative approximation and the user-supplied
C         solution divided by the 2-norm of the user-supplied solution
C         is less than TOL.  Note that this requires the user to set up
C         the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine.
C         The routine with this declaration should be loaded before the
C         stop test so that the correct length is used by the loader.
C         This procedure is not standard Fortran and may not work
C         correctly on your system (although it has worked on every
C         system the authors have tried).  If ITOL is not 11 then this
C         common block is indeed standard Fortran.
C TOL    :IN       Double Precision.
C         Convergence criterion, as described above.
C ITMAX  :IN       Integer.
C         Maximum number of iterations.
C ITER   :OUT      Integer.
C         Number of iterations required to reach convergence, or
C         ITMAX+1 if convergence criterion could not be achieved in
C         ITMAX iterations.
C ERR    :OUT      Double Precision.
C         Error estimate of error in final approximate solution, as
C         defined by ITOL.
C IERR   :OUT      Integer.
C         Return error flag.
C           IERR = 0 => All went well.
C           IERR = 1 => Insufficient storage allocated
C                       for WORK or IWORK.
C           IERR = 2 => Method failed to converge in
C                       ITMAX steps.
C           IERR = 3 => Error in user input.  Check input
C                       value of N, ITOL.
C           IERR = 4 => User error tolerance set too tight.
C                       Reset to 500.0*D1MACH(3).  Iteration proceeded.
C           IERR = 5 => Preconditioning matrix, M,  is not
C                       Positive Definite.  $(r,z) < 0.0$.
C           IERR = 6 => Matrix A is not Positive Definite.
C                       $(p,Ap) < 0.0$.
C IUNIT  :IN       Integer.
C         Unit number on which to write the error at each iteration,
C         if this is desired for monitoring convergence.  If unit
C         number is 0, no writing will occur.
C RWORK  :WORK     Double Precision RWORK(LENW).
C         Double Precision array used for workspace.
C LENW   :IN       Integer.
C         Length of the double precision workspace, RWORK.  LENW >= 4*N.
C IWORK  :WORK     Integer IWORK(LENIW).
C         Used to hold pointers into the double precision workspace,
C         RWORK. Upon return the following locations of IWORK hold
C         information which may be of use to the user:
C         IWORK(9)  Amount of Integer workspace actually used.
C         IWORK(10) Amount of Double Precision workspace actually used.
C LENIW  :IN       Integer.
C         Length of the integer workspace, IWORK.  LENIW >= 10.
C
C *Description:
C       Jacobi's method solves the linear system Ax=b with the
C       basic iterative method (where A = L + D + U):
C
C            n+1    -1       n    n
C           X    = D  (B - LX - UX )
C
C                   n    -1       n
C                = X  + D  (B - AX )
C
C       The Sparse Linear Algebra Package (SLAP) utilizes two matrix
C       data structures: 1) the  SLAP Triad  format or  2)  the SLAP
C       Column format.  The user can hand this routine either of the
C       of these data structures and SLAP  will figure out which one
C       is being used and act accordingly.
C
C       =================== S L A P Triad format ===================
C
C       This routine requires that the  matrix A be   stored in  the
C       SLAP  Triad format.  In  this format only the non-zeros  are
C       stored.  They may appear in  *ANY* order.  The user supplies
C       three arrays of  length NELT, where  NELT is  the number  of
C       non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)).  For
C       each non-zero the user puts the row and column index of that
C       matrix element  in the IA and  JA arrays.  The  value of the
C       non-zero   matrix  element is  placed  in  the corresponding
C       location of the A array.   This is  an  extremely  easy data
C       structure to generate.  On  the  other hand it   is  not too
C       efficient on vector computers for  the iterative solution of
C       linear systems.  Hence,   SLAP changes   this  input    data
C       structure to the SLAP Column format  for  the iteration (but
C       does not change it back).
C
C       Here is an example of the  SLAP Triad   storage format for a
C       5x5 Matrix.  Recall that the entries may appear in any order.
C
C           5x5 Matrix       SLAP Triad format for 5x5 matrix on left.
C                              1  2  3  4  5  6  7  8  9 10 11
C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
C       | 0  0  0 44  0|
C
C       =================== S L A P Column format ==================
C       This routine requires  that the  matrix  A be  stored in the
C       SLAP Column format.  In this format the non-zeros are stored
C       counting down columns (except for  the diagonal entry, which
C       must appear first in each  "column")  and are stored  in the
C       double precision array A.   In other words,  for each column
C       in the matrix put the diagonal entry in  A.  Then put in the
C       other non-zero  elements going down  the column (except  the
C       diagonal) in order.   The  IA array holds the  row index for
C       each non-zero.  The JA array holds the offsets  into the IA,
C       A arrays  for  the  beginning  of each   column.   That  is,
C       IA(JA(ICOL)),  A(JA(ICOL)) points   to the beginning  of the
C       ICOL-th   column    in    IA and   A.      IA(JA(ICOL+1)-1),
C       A(JA(ICOL+1)-1) points to  the  end of the   ICOL-th column.
C       Note that we always have  JA(N+1) = NELT+1,  where N is  the
C       number of columns in  the matrix and NELT  is the number  of
C       non-zeros in the matrix.
C
C       Here is an example of the  SLAP Column  storage format for a
C       5x5 Matrix (in the A and IA arrays '|'  denotes the end of a
C       column):
C
C       5x5 Matrix      SLAP Column format for 5x5 matrix on left.
C       1  2  3    4  5    6  7    8    9 10 11
C       |11 12  0  0 15|   A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
C       |21 22  0  0  0|  IA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  JA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C *Precision:           Double Precision
C *Side Effects:
C       The SLAP Triad format (IA, JA, A) is modified internally to be
C       the SLAP Column format.  See above.
C
C *See Also:
C       DSGS, DIR
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DS2Y, DDCHKW, DSDS, DIR, DSMV, DSDI
C***END PROLOGUE  DSJAC
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
      INTEGER ITER, IUNIT, LENW, IWORK(LENIW), LENIW
      DOUBLE PRECISION B(N), X(N), A(NELT), RWORK(LENW)
      EXTERNAL DSMV, DSDI
      PARAMETER(LOCRB=1,LOCIB=11)
C
C         Compute the inverse of the diagonal of the matrix.  This
C         will be used as the precontioner.
C***FIRST EXECUTABLE STATEMENT  DSJAC
      IERR = 0
      IF( N.LT.1 .OR. NELT.LT.1 ) THEN
         IERR = 3
         RETURN
      ENDIF
      LOCIW = LOCIB
      LOCD = LOCRB
      LOCR = LOCD + N
      LOCZ = LOCR + N
      LOCDZ = LOCZ + N
      LOCW = LOCDZ + N
C
C         Check the workspace allocations.
      CALL DCHKW( 'DSJAC', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR )
      IF( IERR.NE.0 ) RETURN
C
      IWORK(4) = LOCD
      IWORK(9) = LOCIW
      IWORK(10) = LOCW
C
      CALL DS2Y(N, NELT, IA, JA, A, ISYM )
      CALL DSDS(N, NELT, IA, JA, A, ISYM, RWORK(LOCD))
C
C         Set up the work array and perform the iterative refinement.
      CALL DIR(N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSDI, ITOL, TOL,
     $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK(LOCR), RWORK(LOCZ),
     $     RWORK(LOCDZ), RWORK, IWORK )
      RETURN
C------------- LAST LINE OF DSJAC FOLLOWS -----------------------------
      END
*DECK DSGS
      SUBROUTINE DSGS(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
     $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW )
C***BEGIN PROLOGUE  DSGS
C***DATE WRITTEN   871119   (YYMMDD)
C***REVISION DATE  881213   (YYMMDD)
C***CATEGORY NO.  D2A4, D2B4
C***KEYWORDS  LIBRARY=SLATEC(SLAP),
C             TYPE=DOUBLE PRECISION(DSGS-S),
C             Linear system, Sparse, Iterative Precondition
C***AUTHOR  Greenbaum, Anne, Courant Institute
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-300
C             Livermore, CA 94550 (415) 423-3141
C             seager@lll-crg.llnl.gov
C***PURPOSE  Gauss-Seidel method iterative sparse Ax = b solver.
C            Routine to solve a general linear system  Ax = b  using
C            Gauss-Seidel iteration.
C***DESCRIPTION
C *Usage:
C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
C     INTEGER ITER, IERR, IUNIT, LENW, IWORK(NEL+2*N+1), LENIW
C     DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(NEL+3*N)
C
C     CALL DSGS(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
C    $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW )
C
C *Arguments:
C N      :IN       Integer.
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right-hand side vector.
C X      :INOUT    Double Precision X(N).
C         On input X is your initial guess for solution vector.
C         On output X is the final approximate solution.
C NELT   :IN       Integer.
C         Number of Non-Zeros stored in A.
C IA     :INOUT    Integer IA(NELT).
C JA     :INOUT    Integer JA(NELT).
C A      :INOUT    Integer A(NELT).
C         These arrays should hold the matrix A in either the SLAP
C         Triad format or the SLAP Column format.  See "Description",
C         below.  If the SLAP Triad format is chosen it is changed
C         internally to the SLAP Column format.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all nonzero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the lower
C         lower triangle of the matrix is stored.
C ITOL   :IN       Integer.
C         Flag to indicate type of convergence criterion.
C         If ITOL=1, iteration stops when the 2-norm of the residual
C         divided by the 2-norm of the right-hand side is less than TOL.
C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
C         residual divided by the 2-norm of M-inv times the right hand
C         side is less than TOL, where M-inv is the inverse of the
C         diagonal of A.
C         ITOL=11 is often useful for checking and comparing different
C         routines.  For this case, the user must supply the "exact"
C         solution or a very accurate approximation (one with an error
C         much less than TOL) through a common block,
C                     COMMON /SOLBLK/ SOLN(1)
C         if ITOL=11, iteration stops when the 2-norm of the difference
C         between the iterative approximation and the user-supplied
C         solution divided by the 2-norm of the user-supplied solution
C         is less than TOL.  Note that this requires the user to set up
C         the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine.
C         The routine with this declaration should be loaded before the
C         stop test so that the correct length is used by the loader.
C         This procedure is not standard Fortran and may not work
C         correctly on your system (although it has worked on every
C         system the authors have tried).  If ITOL is not 11 then this
C         common block is indeed standard Fortran.
C TOL    :IN       Double Precision.
C         Convergence criterion, as described above.
C ITMAX  :IN       Integer.
C         Maximum number of iterations.
C ITER   :OUT      Integer.
C         Number of iterations required to reach convergence, or
C         ITMAX+1 if convergence criterion could not be achieved in
C         ITMAX iterations.
C ERR    :OUT      Double Precision.
C         Error estimate of error in final approximate solution, as
C         defined by ITOL.
C IERR   :OUT      Integer.
C         Return error flag.
C           IERR = 0 => All went well.
C           IERR = 1 => Insufficient storage allocated
C                       for WORK or IWORK.
C           IERR = 2 => Method failed to converge in
C                       ITMAX steps.
C           IERR = 3 => Error in user input.  Check input
C                       value of N, ITOL.
C           IERR = 4 => User error tolerance set too tight.
C                       Reset to 500.0*D1MACH(3).  Iteration proceeded.
C           IERR = 5 => Preconditioning matrix, M,  is not
C                       Positive Definite.  $(r,z) < 0.0$.
C           IERR = 6 => Matrix A is not Positive Definite.
C                       $(p,Ap) < 0.0$.
C IUNIT  :IN       Integer.
C         Unit number on which to write the error at each iteration,
C         if this is desired for monitoring convergence.  If unit
C         number is 0, no writing will occur.
C RWORK  :WORK     Double Precision RWORK(LENW).
C         Double Precision array used for workspace.  NEL is the number
C         of non-zeros in the lower triangle of the matrix (including
C         the diagonal).
C LENW   :IN       Integer.
C         Length of the double precision workspace, RWORK.
C         LENW >= NEL+3*N.
C IWORK  :WORK     Integer IWORK(LENIW).
C         Integer array used for workspace.  NEL is the number of non-
C         zeros in the lower triangle of the matrix (including the
C         diagonal).
C         Upon return the following locations of IWORK hold information
C         which may be of use to the user:
C         IWORK(9)  Amount of Integer workspace actually used.
C         IWORK(10) Amount of Double Precision workspace actually used.
C LENIW  :IN       Integer.
C         Length of the integer workspace, IWORK.  LENIW >=
C         NEL+N+11.
C
C *Description
C       The Sparse Linear Algebra Package (SLAP) utilizes two matrix
C       data structures: 1) the  SLAP Triad  format or  2)  the SLAP
C       Column format.  The user can hand this routine either of the
C       of these data structures and SLAP  will figure out  which on
C       is being used and act accordingly.
C
C       =================== S L A P Triad format ===================
C
C       This routine requires that the  matrix A be   stored in  the
C       SLAP  Triad format.  In  this format only the non-zeros  are
C       stored.  They may appear in  *ANY* order.  The user supplies
C       three arrays of  length NELT, where  NELT is  the number  of
C       non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)).  For
C       each non-zero the user puts the row and column index of that
C       matrix element  in the IA and  JA arrays.  The  value of the
C       non-zero   matrix  element is  placed  in  the corresponding
C       location of the A array.   This is  an  extremely  easy data
C       structure to generate.  On  the  other hand it   is  not too
C       efficient on vector computers for  the iterative solution of
C       linear systems.  Hence,   SLAP changes   this  input    data
C       structure to the SLAP Column format  for  the iteration (but
C       does not change it back).
C
C       Here is an example of the  SLAP Triad   storage format for a
C       5x5 Matrix.  Recall that the entries may appear in any order.
C
C           5x5 Matrix       SLAP Triad format for 5x5 matrix on left.
C                              1  2  3  4  5  6  7  8  9 10 11
C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C       =================== S L A P Column format ==================
C       This routine requires  that the  matrix  A be  stored in the
C       SLAP Column format.  In this format the non-zeros are stored
C       counting down columns (except for  the diagonal entry, which
C       must appear first in each  "column")  and are stored  in the
C       double precision array A.   In other words,  for each column
C       in the matrix put the diagonal entry in  A.  Then put in the
C       other non-zero  elements going down  the column (except  the
C       diagonal) in order.   The  IA array holds the  row index for
C       each non-zero.  The JA array holds the offsets  into the IA,
C       A arrays  for  the  beginning  of each   column.   That  is,
C       IA(JA(ICOL)),  A(JA(ICOL)) points   to the beginning  of the
C       ICOL-th   column    in    IA and   A.      IA(JA(ICOL+1)-1),
C       A(JA(ICOL+1)-1) points to  the  end of the   ICOL-th column.
C       Note that we always have  JA(N+1) = NELT+1,  where N is  the
C       number of columns in  the matrix and NELT  is the number  of
C       non-zeros in the matrix.
C
C       Here is an example of the  SLAP Column  storage format for a
C       5x5 Matrix (in the A and IA arrays '|'  denotes the end of a
C       column):
C
C       5x5 Matrix      SLAP Column format for 5x5 matrix on left.
C       1  2  3    4  5    6  7    8    9 10 11
C       |11 12  0  0 15|   A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
C       |21 22  0  0  0|  IA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  JA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C *Precision:           Double Precision
C *Side Effects:
C       The SLAP Triad format (IA, JA, A) is modified internally to be
C       the SLAP Column format.  See above.
C
C *See Also:
C       DSJAC, DIR
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DS2Y, DCHKW, DS2LT, SDIR, DSMV, DSLI
C***END PROLOGUE  DSGS
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
      INTEGER ITER, IUNIT, IWORK(10)
      DOUBLE PRECISION B(N), X(N), A(N), TOL, ERR, RWORK(1)
      EXTERNAL DSMV, DSLI
      PARAMETER(LOCRB=1,LOCIB=11)
C
C         Modify the SLAP matrix data structure to YSMP-Column.
C***FIRST EXECUTABLE STATEMENT  DSGS
      IF( N.LT.1 .OR. NELT.LT.1 ) THEN
         IERR = 3
         RETURN
      ENDIF
      CALL DS2Y( N, NELT, IA, JA, A, ISYM )
C
C         Count number of elements in lower triangle of the matrix.
      IF( ISYM.EQ.0 ) THEN
         NEL = 0
         DO 20 ICOL = 1, N
            JBGN = JA(ICOL)
            JEND = JA(ICOL+1)-1
            DO 10 J = JBGN, JEND
               IF( IA(J).GE.ICOL ) NEL = NEL + 1
 10         CONTINUE
 20      CONTINUE
      ELSE
         NEL = JA(N+1)-1
      ENDIF
C
C         Set up the work arrays.  Then store the lower triangle of
C         the matrix.
C
      LOCJEL = LOCIB
      LOCIEL = LOCJEL + N+1
      LOCIW = LOCIEL + NEL
C
      LOCEL = LOCRB
      LOCR = LOCEL + NEL
      LOCZ = LOCR + N
      LOCDZ = LOCZ + N
      LOCW = LOCDZ + N
C
C         Check the workspace allocations.
      CALL DCHKW( 'DSGS', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR )
      IF( IERR.NE.0 ) RETURN
C
      IWORK(1) = NEL
      IWORK(2) = LOCIEL
      IWORK(3) = LOCJEL
      IWORK(4) = LOCEL
      IWORK(9) = LOCIW
      IWORK(10) = LOCW
C
      CALL DS2LT( N, NELT, IA, JA, A, ISYM, NEL, IWORK(LOCIEL),
     $     IWORK(LOCJEL), RWORK(LOCEL) )
C
C         Call iterative refinement routine.
      CALL DIR(N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSLI,
     $     ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK(LOCR),
     $     RWORK(LOCZ), RWORK(LOCDZ), RWORK, IWORK )
C
C         Set the amount of Integer and Double Precision Workspace used.
      IWORK(9) = LOCIW+N+NELT
      IWORK(10) = LOCW+NELT
      RETURN
C------------- LAST LINE OF DSGS FOLLOWS ------------------------------
      END
*DECK DSILUR
      SUBROUTINE DSILUR(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
     $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW )
C***BEGIN PROLOGUE  DSILUR
C***DATE WRITTEN   871119   (YYMMDD)
C***REVISION DATE  881213   (YYMMDD)
C***CATEGORY NO.  D2A4, D2B4
C***KEYWORDS  LIBRARY=SLATEC(SLAP),
C             TYPE=DOUBLE PRECISION(DSILUR-S),
C             Linear system, Sparse, Iterative Precondition
C***AUTHOR  Greenbaum, Anne, Courant Institute
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-300
C             Livermore, CA 94550 (415) 423-3141
C             seager@lll-crg.llnl.gov
C***PURPOSE  Incomplete LU Iterative Refinement sparse Ax = b solver.
C            Routine to solve a general linear system  Ax = b  using
C            the incomplete LU decomposition with iterative refinement.
C***DESCRIPTION
C *Usage:
C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
C     INTEGER ITER, IERR, IUNIT, LENW, IWORK(NEL+NU+4*N+2), LENIW
C     DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(NEL+NU+4*N)
C
C     CALL DSILUR(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
C    $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW)
C
C *Arguments:
C N      :IN       Integer
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right-hand side vector.
C X      :INOUT    Double Precision X(N).
C         On input X is your initial guess for solution vector.
C         On output X is the final approximate solution.
C NELT   :IN       Integer.
C         Number of Non-Zeros stored in A.
C IA     :INOUT    Integer IA(NELT).
C JA     :INOUT    Integer JA(NELT).
C A      :INOUT    Integer A(NELT).
C         These arrays should hold the matrix A in either the SLAP
C         Triad format or the SLAP Column format.  See "Description",
C         below.  If the SLAP Triad format is chosen it is changed
C         internally to the SLAP Column format.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all nonzero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the upper
C         or lower triangle of the matrix is stored.
C ITOL   :IN       Integer.
C         Flag to indicate type of convergence criterion.
C         If ITOL=1, iteration stops when the 2-norm of the residual
C         divided by the 2-norm of the right-hand side is less than TOL.
C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
C         residual divided by the 2-norm of M-inv times the right hand
C         side is less than TOL, where M-inv is the inverse of the
C         diagonal of A.
C         ITOL=11 is often useful for checking and comparing different
C         routines.  For this case, the user must supply the "exact"
C         solution or a very accurate approximation (one with an error
C         much less than TOL) through a common block,
C                     COMMON /SOLBLK/ SOLN(1)
C         if ITOL=11, iteration stops when the 2-norm of the difference
C         between the iterative approximation and the user-supplied
C         solution divided by the 2-norm of the user-supplied solution
C         is less than TOL.  Note that this requires the user to set up
C         the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine.
C         The routine with this declaration should be loaded before the
C         stop test so that the correct length is used by the loader.
C         This procedure is not standard Fortran and may not work
C         correctly on your system (although it has worked on every
C         system the authors have tried).  If ITOL is not 11 then this
C         common block is indeed standard Fortran.
C TOL    :IN       Double Precision.
C         Convergence criterion, as described above.
C ITMAX  :IN       Integer.
C         Maximum number of iterations.
C ITER   :OUT      Integer.
C         Number of iterations required to reach convergence, or
C         ITMAX+1 if convergence criterion could not be achieved in
C         ITMAX iterations.
C ERR    :OUT      Double Precision.
C         Error estimate of error in final approximate solution, as
C         defined by ITOL.
C IERR   :OUT      Integer.
C         Return error flag.
C           IERR = 0 => All went well.
C           IERR = 1 => Insufficient storage allocated
C                       for WORK or IWORK.
C           IERR = 2 => Method failed to converge in
C                       ITMAX steps.
C           IERR = 3 => Error in user input.  Check input
C                       value of N, ITOL.
C           IERR = 4 => User error tolerance set too tight.
C                       Reset to 500.0*D1MACH(3).  Iteration proceeded.
C           IERR = 5 => Preconditioning matrix, M,  is not
C                       Positive Definite.  $(r,z) < 0.0$.
C           IERR = 6 => Matrix A is not Positive Definite.
C                       $(p,Ap) < 0.0$.
C           IERR = 7 => Incomplete factorization broke down
C                       and was fudged.  Resulting preconditioning may
C                       be less than the best.
C IUNIT  :IN       Integer.
C         Unit number on which to write the error at each iteration,
C         if this is desired for monitoring convergence.  If unit
C         number is 0, no writing will occur.
C RWORK  :WORK     Double Precision RWORK(LENW).
C         Double Precision array used for workspace.  NEL is the number
C         of non-zeros in the lower triangle of the matrix (including
C         the diagonal).  NU is the number of nonzeros in the upper
C         triangle of the matrix (including the diagonal).
C LENW   :IN       Integer.
C         Length of the double precision workspace, RWORK.
C         LENW >= NEL+NU+4*N.
C IWORK  :WORK     Integer IWORK(LENIW).
C         Integer array used for workspace.  NEL is the number of non-
C         zeros in the lower triangle of the matrix (including the
C         diagonal).  NU is the number of nonzeros in the upper
C         triangle of the matrix (including the diagonal).
C         Upon return the following locations of IWORK hold information
C         which may be of use to the user:
C         IWORK(9)  Amount of Integer workspace actually used.
C         IWORK(10) Amount of Double Precision workspace actually used.
C LENIW  :IN       Integer.
C         Length of the integer workspace, IWORK.  LENIW >=
C         NEL+NU+4*N+10.
C
C *Description
C       The Sparse Linear Algebra Package (SLAP) utilizes two matrix
C       data structures: 1) the  SLAP Triad  format or  2)  the SLAP
C       Column format.  The user can hand this routine either of the
C       of these data structures and SLAP  will figure out  which on
C       is being used and act accordingly.
C
C       =================== S L A P Triad format ===================
C
C       This routine requires that the  matrix A be   stored in  the
C       SLAP  Triad format.  In  this format only the non-zeros  are
C       stored.  They may appear in  *ANY* order.  The user supplies
C       three arrays of  length NELT, where  NELT is  the number  of
C       non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)).  For
C       each non-zero the user puts the row and column index of that
C       matrix element  in the IA and  JA arrays.  The  value of the
C       non-zero   matrix  element is  placed  in  the corresponding
C       location of the A array.   This is  an  extremely  easy data
C       structure to generate.  On  the  other hand it   is  not too
C       efficient on vector computers for  the iterative solution of
C       linear systems.  Hence,   SLAP changes   this  input    data
C       structure to the SLAP Column format  for  the iteration (but
C       does not change it back).
C
C       Here is an example of the  SLAP Triad   storage format for a
C       5x5 Matrix.  Recall that the entries may appear in any order.
C
C           5x5 Matrix       SLAP Triad format for 5x5 matrix on left.
C                              1  2  3  4  5  6  7  8  9 10 11
C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C       =================== S L A P Column format ==================
C       This routine requires  that the  matrix  A be  stored in the
C       SLAP Column format.  In this format the non-zeros are stored
C       counting down columns (except for  the diagonal entry, which
C       must appear first in each  "column")  and are stored  in the
C       double precision array A.   In other words,  for each column
C       in the matrix put the diagonal entry in  A.  Then put in the
C       other non-zero  elements going down  the column (except  the
C       diagonal) in order.   The  IA array holds the  row index for
C       each non-zero.  The JA array holds the offsets  into the IA,
C       A arrays  for  the  beginning  of each   column.   That  is,
C       IA(JA(ICOL)),  A(JA(ICOL)) points   to the beginning  of the
C       ICOL-th   column    in    IA and   A.      IA(JA(ICOL+1)-1),
C       A(JA(ICOL+1)-1) points to  the  end of the   ICOL-th column.
C       Note that we always have  JA(N+1) = NELT+1,  where N is  the
C       number of columns in  the matrix and NELT  is the number  of
C       non-zeros in the matrix.
C
C       Here is an example of the  SLAP Column  storage format for a
C       5x5 Matrix (in the A and IA arrays '|'  denotes the end of a
C       column):
C
C       5x5 Matrix      SLAP Column format for 5x5 matrix on left.
C       1  2  3    4  5    6  7    8    9 10 11
C       |11 12  0  0 15|   A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
C       |21 22  0  0  0|  IA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  JA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C *Precision:           Double Precision
C *Side Effects:
C       The SLAP Triad format (IA, JA, A) is modified internally to be
C       the SLAP Column format.  See above.
C
C *Portability:
C       DSJAC, DSGS, DIR
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DS2Y, DCHKW, DSILUS, DIR, DSMV, DSLUI
C***END PROLOGUE  DSILUR
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER
      INTEGER IERR, IUNIT, LENW, IWORK(LENIW), LENIW
      DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(LENW)
      PARAMETER (LOCRB=1, LOCIB=11)
C
      EXTERNAL DSMV, DSLUI
C
C         Change the SLAP input matrix IA, JA, A to SLAP-Column format.
C***FIRST EXECUTABLE STATEMENT  DSILUR
      IERR = 0
      IF( N.LT.1 .OR. NELT.LT.1 ) THEN
         IERR = 3
         RETURN
      ENDIF
      CALL DS2Y( N, NELT, IA, JA, A, ISYM )
C
C         Count number of Non-Zero elements in preconditioner ILU
C         matrix.  Then set up the work arrays.
      NL = 0
      NU = 0
      DO 20 ICOL = 1, N
C         Don't count diagonal.
         JBGN = JA(ICOL)+1
         JEND = JA(ICOL+1)-1
         IF( JBGN.LE.JEND ) THEN
CVD$ NOVECTOR
            DO 10 J = JBGN, JEND
               IF( IA(J).GT.ICOL ) THEN
                  NL = NL + 1
                  IF( ISYM.NE.0 ) NU = NU + 1
               ELSE
                  NU = NU + 1
               ENDIF
 10         CONTINUE
         ENDIF
 20   CONTINUE
C
      LOCIL = LOCIB
      LOCJL = LOCIL + N+1
      LOCIU = LOCJL + NL
      LOCJU = LOCIU + NU
      LOCNR = LOCJU + N+1
      LOCNC = LOCNR + N
      LOCIW = LOCNC + N
C
      LOCL = LOCRB
      LOCDIN = LOCL + NL
      LOCU = LOCDIN + N
      LOCR = LOCU + NU
      LOCZ = LOCR + N
      LOCDZ = LOCZ + N
      LOCW = LOCDZ + N
C
C         Check the workspace allocations.
      CALL DCHKW( 'DSILUR', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR )
      IF( IERR.NE.0 ) RETURN
C
      IWORK(1) = LOCIL
      IWORK(2) = LOCJL
      IWORK(3) = LOCIU
      IWORK(4) = LOCJU
      IWORK(5) = LOCL
      IWORK(6) = LOCDIN
      IWORK(7) = LOCU
      IWORK(9) = LOCIW
      IWORK(10) = LOCW
C
C         Compute the Incomplete LU decomposition.
      CALL DSILUS( N, NELT, IA, JA, A, ISYM, NL, IWORK(LOCIL),
     $     IWORK(LOCJL), RWORK(LOCL), RWORK(LOCDIN), NU, IWORK(LOCIU),
     $     IWORK(LOCJU), RWORK(LOCU), IWORK(LOCNR), IWORK(LOCNC) )
C
C         Do the Preconditioned Iterative Refinement iteration.
      CALL DIR(N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSLUI,
     $     ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK(LOCR),
     $     RWORK(LOCZ), RWORK(LOCDZ), RWORK, IWORK)
      RETURN
C------------- LAST LINE OF DSILUR FOLLOWS ----------------------------
      END
*DECK ISDIR
      FUNCTION ISDIR(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL,
     $     ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK, IWORK,
     $     BNRM, SOLNRM)
C***BEGIN PROLOGUE  ISDIR
C***REFER TO  DIR, DSJAC, DSGS
C***DATE WRITTEN   871119   (YYMMDD)
C***REVISION DATE  880320   (YYMMDD)
C***CATEGORY NO.  D2A4, D2B4
C***KEYWORDS  LIBRARY=SLATEC(SLAP),
C             TYPE=DOUBLE PRECISION(ISDIR-S),
C             Linear system, Sparse, Stop Test
C***AUTHOR  Greenbaum, Anne, Courant Institute
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-300
C             Livermore, CA 94550 (415) 423-3141
C             seager@lll-crg.llnl.gov
C***PURPOSE  Preconditioned Iterative Refinement Stop Test.
C            This routine calculates the stop test for the iterative
C            refinement iteration scheme.  It returns a nonzero if the
C            error estimate (the type of which is determined by ITOL)
C            is less than the user specified tolerance TOL.
C***DESCRIPTION
C *Usage:
C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER
C     INTEGER IERR, IUNIT, IWORK(USER DEFINED)
C     DOUBLE PRECISION B(N), X(N), A(N), TOL, ERR, R(N), Z(N), DZ(N)
C     DOUBLE PRECISION RWORK(USER DEFINED), BNRM, SOLNRM
C     EXTERNAL MSOLVE
C
C     IF( ISDIR(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL,
C    $     ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK, IWORK,
C    $     BNRM, SOLNRM) .NE. 0 ) THEN ITERATION DONE
C
C *Arguments:
C N      :IN       Integer.
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right-hand side vector.
C X      :IN       Double Precision X(N).
C         The current approximate solution vector.
C NELT   :IN       Integer.
C         Number of Non-Zeros stored in A.
C IA     :IN       Integer IA(NELT).
C JA     :IN       Integer JA(NELT).
C A      :IN       Double Precision A(NELT).
C         These arrays contain the matrix data structure for A.
C         It could take any form.  See "C *Description" in the
C         DIR routine.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all nonzero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the upper
C         or lower triangle of the matrix is stored.
C MSOLVE :EXT      External.
C         Name of a routine which solves a linear system Mz = r for
C         z given r with the preconditioning matrix M (M is supplied via
C         RWORK and IWORK arrays.  The name of the MSOLVE routine must
C         be declared external in the calling program.  The calling
C         sequence to MSOLVE is:
C             CALL MSOLVE(N, R, Z, RWORK, IWORK)
C         Where N is the number of unknowns, R is the right-hand side
C         vector, and z is the solution upon return.  RWORK is a double
C         precision array that can be used to pass necessary
C         preconditioning information and/or workspace to MSOLVE.
C         IWORK is an integer work array for the same purpose as RWORK.
C ITOL   :IN       Integer.
C         Flag to indicate type of convergence criterion.
C         If ITOL=1, iteration stops when the 2-norm of the residual
C         divided by the 2-norm of the right-hand side is less than TOL.
C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
C         residual divided by the 2-norm of M-inv times the right hand
C         side is less than TOL, where M-inv is the inverse of the
C         diagonal of A.
C         ITOL=11 is often useful for checking and comparing different
C         routines.  For this case, the user must supply the "exact"
C         solution or a very accurate approximation (one with an error
C         much less than TOL) through a common block,
C                     COMMON /SOLBLK/ SOLN(1)
C         if ITOL=11, iteration stops when the 2-norm of the difference
C         between the iterative approximation and the user-supplied
C         solution divided by the 2-norm of the user-supplied solution
C         is less than TOL.  Note that this requires the user to set up
C         the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine.
C         The routine with this declaration should be loaded before the
C         stop test so that the correct length is used by the loader.
C         This procedure is not standard Fortran and may not work
C         correctly on your system (although it has worked on every
C         system the authors have tried).  If ITOL is not 11 then this
C         common block is indeed standard Fortran.
C TOL    :IN       Double Precision.
C         Convergence criterion, as described above.
C ITER   :IN       Integer.
C         Current iteration count.
C ITMAX  :IN       Integer.
C         Maximum number of iterations.
C ERR    :OUT      Double Precision.
C         Error estimate of error in the X(N) approximate solution, as
C         defined by ITOL.
C IERR   :OUT      Integer.
C         Error flag.  IERR is set to 3 if ITOL is not on of the
C         acceptable values, see above.
C IUNIT  :IN       Integer.
C         Unit number on which to write the error at each iteration,
C         if this is desired for monitoring convergence.  If unit
C         number is 0, no writing will occur.
C R      :IN       Double Precision R(N).
C         The residual R = B-AX.
C Z      :WORK     Double Precision Z(N).
C         Workspace used to hold the pseudo-residual M z = r.
C DZ     :WORK     Double Precision DZ(N).
C         Workspace used to hold temporary vector(s).
C RWORK  :WORK     Double Precision RWORK(USER DEFINABLE).
C         Double Precision array that can be used by  MSOLVE.
C IWORK  :WORK     Integer IWORK(USER DEFINABLE).
C         Integer array that can be used by MSOLVE.
C BNRM   :INOUT    Double Precision.
C         Norm of the right hand side.  Type of norm depends on ITOL.
C         Calculated only on the first call.
C SOLNRM :INOUT    Double Precision.
C         2-Norm of the true solution, SOLN.  Only computed and used
C         if ITOL = 11.
C
C *Function Return Values:
C       0 : Error estimate (determined by ITOL) is *NOT* less than the
C           specified tolerance, TOL.  The iteration must continue.
C       1 : Error estimate (determined by ITOL) is less than the
C           specified tolerance, TOL.  The iteration can be considered
C           complete.
C
C *Precision:           Double Precision
C *See Also:
C       DIR, DSJAC, DSGS
C
C *Cautions:
C     This routine will attempt to write to the fortran logical output
C     unit IUNIT, if IUNIT .ne. 0.  Thus, the user must make sure that
C     this  logical  unit  must  be  attached  to  a  file or terminal
C     before calling this routine with a non-zero  value  for   IUNIT.
C     This routine does not check for the validity of a non-zero IUNIT
C     unit number.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  MSOLVE, DNRM2
C***COMMON BLOCKS    SOLBLK
C***END PROLOGUE  ISDIR
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER
      INTEGER IUNIT, IWORK(*)
      DOUBLE PRECISION B(N), X(N), A(NELT), R(N), Z(N), DZ(N), RWORK(*)
      EXTERNAL MSOLVE
      COMMON /SOLBLK/ SOLN(1)
C
C***FIRST EXECUTABLE STATEMENT  ISDIR
      ISDIR = 0
      IF( ITOL.EQ.1 ) THEN
C         err = ||Residual||/||RightHandSide|| (2-Norms).
         IF(ITER .EQ. 0) BNRM = DNRM2(N, B, 1)
         ERR = DNRM2(N, R, 1)/BNRM
      ELSE IF( ITOL.EQ.2 ) THEN
C                  -1              -1
C         err = ||M  Residual||/||M  RightHandSide|| (2-Norms).
         IF(ITER .EQ. 0) THEN
            CALL MSOLVE(N, B, DZ, NELT, IA, JA, A, ISYM, RWORK, IWORK)
            BNRM = DNRM2(N, DZ, 1)
         ENDIF
         ERR = DNRM2(N, Z, 1)/BNRM
      ELSE IF( ITOL.EQ.11 ) THEN
C         err = ||x-TrueSolution||/||TrueSolution|| (2-Norms).
         IF( ITER.EQ.0 ) SOLNRM = DNRM2(N, SOLN, 1)
         DO 10 I = 1, N
            DZ(I) = X(I) - SOLN(I)
 10      CONTINUE
         ERR = DNRM2(N, DZ, 1)/SOLNRM
      ELSE
C
C         If we get here ITOL is not one of the acceptable values.
         ERR = 1.0E10
         IERR = 3
      ENDIF
C
      IF( IUNIT.NE.0 ) THEN
         WRITE(IUNIT,1000) ITER,ERR
      ENDIF
C
      IF( ERR.LE.TOL ) ISDIR = 1
C
      RETURN
 1000 FORMAT(5X,'ITER = ',I4,' Error Estimate = ',E16.7)
C------------- LAST LINE OF ISDIR FOLLOWS -----------------------------
      END

max_stars_repo_name: fmyuan/clm-microbe

max_issues_repo_name: fmyuan/clm-microbe

max_forks_repo_name: email-clm/CLM-Microbe


text:
program add_real
      implicit none
      INTEGER, PARAMETER      :: np=10
      REAL, DIMENSION(np)   :: A, B
      INTEGER:: i, getpid

      DO i=1, np
         A(i)=-1
         B(i)=1
      END DO

      !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO
      DO i=1, np
         A(i)=3*B(i)
      END DO
      !$OMP END TARGET TEAMS DISTRIBUTE PARALLEL DO

      DO i=1, np
         IF (A(i) .NE. 3) THEN
           WRITE(*, *) "ERROR AT INDEX ", i, "EXPECT 2 BUT RECEIVED", A(i)
           call kill(getpid(),7)
         ENDIF
      END DO

end program add_real

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
!===================================================================
!  Uses a divergent flow to push parcels toward or away from grid
!  points so as to bring the area fraction at each grid point back
!  to unity.

!  Here tri-diagonal solves are used in the vertical direction.

!           Initialise data first using inipar.f90
!===================================================================

program diverge

 !Import FFT library:
use stafft
use deriv1d

 !Import parameters and constants:
use constants

implicit none

 !Total number of parcels:
integer:: n

 !Parcel positions (x,z) and area fractions:
double precision:: x(nm),z(nm),v(nm)

 !Rates of change (dx/dt,dz/dt):
double precision:: dxdt(nm),dzdt(nm)

 !Gridded area fraction:
double precision:: volg(0:nz,nx)

!Tri-diagonal arrays:
double precision:: ap(nx),apb(nx)
double precision:: etdv(0:nz,nx),htdv(0:nz,nx)
double precision:: etd1(nz-1),htd1(nz-1)
double precision:: etda(nz),htda(nz)

 !Wavenumbers::
double precision:: hrkx(nx),rkx(nx)

 !Quantities needed in FFTs:
double precision:: xtrig(2*nx)
integer:: xfactors(5)

 !Next horizontal grid point used in bilinear interpolation:
integer:: ixp(nx)

 !Time step (does not matter):
double precision,parameter:: dt=0.01d0

double precision:: t
integer:: istep,nsteps

!---------------------------------------------------------------------
! Initialise (read data, define fixed arrays):
call initialise

write(*,*) ' Enter the number of time steps used to relax the'
write(*,*) ' gridded area fraction back to unity:'
read(*,*) nsteps

! Evolve parcels:
do istep=1,nsteps
  call evolve(istep)
enddo

! Finalise:
call par2grid

! Write final data
write(80,rec=nsteps+1) real(nsteps+1),real(volg-one)

close(40)
close(80)

write(*,*)
write(*,*) ' The rms area fraction error is in volf.asc and the'
write(*,*) ' gridded area fraction error is in dv.r4; view with'
write(*,*)
write(*,'(a,i3,1x,i3,a)') '  dataview -ndim ',nx,nz+1,' dv.r4 -mod &'
write(*,*)

contains

!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

subroutine initialise
! Initialises this program

implicit none

double precision:: a0(nx),a0b(nx),ksq(nx)
double precision,parameter:: dzisq=dzi**2
integer,parameter:: nwx=nx/2
integer:: i,ix,kx,iz

!---------------------------------------------------------------------
! Read parcel positions and volume fractions (of grid cell volume):
open(80,file='ini_parcels.asc',status='old')
read(80,*) n
do i=1,n
  read(80,*) x(i),z(i),v(i)
enddo
close(80)

t=zero

!----------------------------------------------------------
! Next grid point indices used in bi-linear interpolation:
do ix=1,nx-1
  ixp(ix)=ix+1
enddo
ixp(nx)=1

!----------------------------------------------------------
! Set up FFTs:
call initfft(nx,xfactors,xtrig)

! Define x wavenumbers:
call init_deriv(nx,ellx,hrkx)
rkx(1)=zero
do kx=1,nwx-1
  rkx(kx+1)   =hrkx(2*kx)
  rkx(nx+1-kx)=hrkx(2*kx)
enddo
rkx(nwx+1)=hrkx(nx)

 !Squared wavenumber array (used in tridiagonal solve):
ksq=rkx**2

!-----------------------------------------------------------------------
! Fixed coefficients used in the tridiagonal problems:
a0=-two*dzisq-f56*ksq
a0b=-dzisq-f13*ksq
ap=dzisq-f112*ksq
apb=dzisq-f16*ksq

!-----------------------------------------------------------------------
! Tridiagonal arrays for inversion of Poisson's equation:
htdv(:,1)=zero
etdv(:,1)=zero
do kx=2,nx
  htdv(0,kx)=one/a0b(kx)
  etdv(0,kx)=-apb(kx)*htdv(0,kx)
  do iz=1,nz-1
    htdv(iz,kx)=one/(a0(kx)+ap(kx)*etdv(iz-1,kx))
    etdv(iz,kx)=-ap(kx)*htdv(iz,kx)
  enddo
  htdv(nz,kx)=one/(a0b(kx)+apb(kx)*etdv(nz-1,kx))
enddo

! Tridiagonal arrays for the compact difference calculation of d/dz
! for fields f for which df/dz = 0 at the boundaries:
htd1(1)=one/f23
etd1(1)=-f16*htd1(1)
do iz=2,nz-2
  htd1(iz)=one/(f23+f16*etd1(iz-1))
  etd1(iz)=-f16*htd1(iz)
enddo
htd1(nz-1)=one/(f23+f16*etd1(nz-2))

! Tridiagonal arrays used for integrating in z (see vertint):
htda(1)=one/f76
etda(1)=-f16*htda(1)
do iz=2,nz-1
  htda(iz)=one/(one+f16*etda(iz-1))
  etda(iz)=-f16*htda(iz)
enddo
htda(nz)=one/(f76+f16*etda(nz-2))

!---------------------------------------------------------------
! File to contain rms area fraction - 1 & max abs value vs time:
open(40,file='volf.asc',status='replace')

! File to contain the gridded area fraction - 1 vs time:
open(80,file='dv.r4',form='unformatted',access='direct', &
                   status='replace',recl=4*(nx*(nz+1)+1))

return
end subroutine initialise

!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

subroutine evolve(istep)
! Evolves the parcel positions (x,z) one time step (Euler method).

implicit none

double precision:: phi(0:nz,nx),ud(0:nz,nx),wd(0:nz,nx)
double precision:: wbar(0:nz)
double precision:: xx,px,pxc,zz,pz,pzc
integer:: istep,ix,iz,kx
integer:: i,ix0,ix1,iz0,iz1
!DGD
double precision,parameter:: pi=acos(-1.d0)
double precision,parameter:: k=two*pi/ellx,m=pi/ellz
double precision:: phie(0:nz,nx),ude(0:nz,nx),wde(0:nz,nx)
!DGD

!-----------------------------------------------------------------
! Compute gridded area fractions:
call par2grid

! Form divergence field * dt and store in phi temporarily:
phi=volg-one

! Write data
write(80,rec=istep) real(istep),real(phi)

!DGD
!px=one
px=one
do ix=1,nx
  xx=k*(dx*dble(ix-1)-hlx)
  do iz=0,nz
    zz=m*(zmin+dz*dble(iz))
    phie(iz,ix)=sin(xx)*sin(zz)
    ude(iz,ix)=k*cos(xx)*sin(zz)
    wde(iz,ix)=m*(px+sin(xx))*cos(zz)
    phi(iz,ix)=-(m**2*px+(k**2+m**2)*sin(xx))*sin(zz)
  enddo
enddo
!DGD

!-----------------------------------------
! Forward x FFT:
call forfft(nz+1,nx,phi,xtrig,xfactors)

! Compute the x-independent part of wd by integration:
call vertint(phi(0,1),wbar)

! Invert Laplace's operator semi-spectrally with compact differences:
call lapinv1(phi)

! Compute x derivative spectrally:
call deriv(nz+1,nx,hrkx,phi,ud)

! Reverse x FFT to define x velocity component ud:
call revfft(nz+1,nx,ud,xtrig,xfactors)

! Compute z derivative by compact differences:
call diffz1(phi,wd)

! Add on the x-independent part of wd:
wd(:,1)=wd(:,1)+wbar

! Reverse x FFT:
call revfft(nz+1,nx,wd,xtrig,xfactors)

!DGD
open(77,file='nduw.r4',form='unformatted',access='direct', &
                     status='replace',recl=4*(nx*(nz+1)+1))
write(77,rec=1) real(0.),real(ud-ude)
write(77,rec=2) real(0.),real(wd-wde)
close(77)

call revfft(nz+1,nx,phi,xtrig,xfactors)
open(77,file='ndphi.r4',form='unformatted',access='direct', &
                      status='replace',recl=4*(nx*(nz+1)+1))
write(77,rec=1) real(0.),real(phi-phie)
close(77)

open(77,file='phi.r4',form='unformatted',access='direct', &
                    status='replace',recl=4*(nx*(nz+1)+1))
write(77,rec=1) real(0.),real(phi)
close(77)

ud=(ud-ude)**2
xx=sqrt((f12*sum(ud(0,:)+ud(nz,:))+sum(ud(1:nz-1,:)))/dble(ncell))
wd=(wd-wde)**2
zz=sqrt(sum(wd(1:nz-1,:))/dble(ncell))
write(*,'(a,i4,2(a,f18.16))') &
     ' nx = ',nx,'   rms ud error = ',xx,'   rms wd error = ',zz
stop
!DGD

!------------------------------------------------------------------
! Increment parcel positions usind (ud,wd) field:
do i=1,n
  xx=dxi*(x(i)-xmin)
  ix0=min(1+int(xx),nx)
  pxc=dble(ix0)-xx
  px=one-pxc
  ix1=ixp(ix0)

  zz=dzi*(z(i)-zmin)
  iz0=min(int(zz),nz-1)
  pz=zz-dble(iz0)
  pzc=one-pz
  iz1=iz0+1

  x(i)=x(i)+pxc*(pzc*ud(iz0,ix0)+pz*ud(iz1,ix0))+ &
            +px*(pzc*ud(iz0,ix1)+pz*ud(iz1,ix1))

  z(i)=z(i)+pxc*(pzc*wd(iz0,ix0)+pz*wd(iz1,ix0))+ &
            +px*(pzc*wd(iz0,ix1)+pz*wd(iz1,ix1))
enddo

! Advance time:
t=t+dt

return
end subroutine evolve

!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

subroutine par2grid

! Reverse bi-linearly interpolates the area fraction (v)
! on all parcels to obtain the gridded area fraction (volg).

 !Declarations:
implicit none

 !Local variables:
double precision:: xx,px,pxc,zz,pz,pzc
integer:: i,ix0,ix1,iz0,iz1

!--------------------------------------------------
volg=zero

do i=1,n
  xx=dxi*(x(i)-xmin)
  ix0=min(1+int(xx),nx)
  pxc=dble(ix0)-xx
  px=one-pxc
  ix1=ixp(ix0)

  zz=dzi*(z(i)-zmin)
  iz0=min(int(zz),nz-1)
  pz=zz-dble(iz0)
  pzc=one-pz
  iz1=iz0+1

  volg(iz0,ix0)=volg(iz0,ix0)+v(i)*pxc*pzc
  volg(iz1,ix0)=volg(iz1,ix0)+v(i)*pxc*pz
  volg(iz0,ix1)=volg(iz0,ix1)+v(i)*px*pzc
  volg(iz1,ix1)=volg(iz1,ix1)+v(i)*px*pz
enddo

 !Double edge values:
do ix0=1,nx
  volg(0 ,ix0)=two*volg(0 ,ix0)
  volg(nz,ix0)=two*volg(nz,ix0)
enddo

 !Write out diagnostics:
write(40,'(f7.2,2(1x,f15.12))') t,sqrt(sum((volg-one)**2)/dble(ncell)), &
                                  maxval(abs(volg-one))

return
end subroutine par2grid

!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

subroutine lapinv1(fs)
 !Inverts Laplace's operator on fs in semi-spectral space.
 !Here dfs/dz = 0 on the z boundaries.
 !Uses 4th-order compact differencing
 !*** Overwrites fs ***

 !Declarations:
implicit none

 !Passed variables:
double precision:: fs(0:nz,nx)

 !Local variables:
double precision:: rs(0:nz,nx)
integer:: kx,iz

fs(:,1)=zero
do kx=2,nx
  rs(0,kx)=f13*fs(0,kx)+f16*fs(1,kx)
  do iz=1,nz-1
    rs(iz,kx)=f112*(fs(iz-1,kx)+fs(iz+1,kx))+f56*fs(iz,kx)
  enddo
  rs(nz,kx)=f13*fs(nz,kx)+f16*fs(nz-1,kx)

  fs(0,kx)=rs(0,kx)*htdv(0,kx)
  do iz=1,nz-1
    fs(iz,kx)=(rs(iz,kx)-ap(kx)*fs(iz-1,kx))*htdv(iz,kx)
  enddo
  fs(nz,kx)=(rs(nz,kx)-apb(kx)*fs(nz-1,kx))*htdv(nz,kx)

  do iz=nz-1,0,-1
    fs(iz,kx)=etdv(iz,kx)*fs(iz+1,kx)+fs(iz,kx)
  enddo
enddo

return
end subroutine

!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

subroutine diffz1(fs,ds)
 !Calculates df/dz for a field f which has df/dz = 0 at the boundaries
 !using 4th-order compact differencing.  Here fs = f and ds = df/dz.

 !Declarations:
implicit none

 !Passed variables:
double precision:: fs(0:nz,nx),ds(0:nz,nx)

 !Local variables:
integer:: ix,iz

do ix=1,nx
  do iz=1,nz-1
    ds(iz,ix)=(fs(iz+1,ix)-fs(iz-1,ix))*hdzi
  enddo

  ds(0,ix)=zero
  ds(1,ix)=ds(1,ix)*htd1(1)
  do iz=2,nz-1
    ds(iz,ix)=(ds(iz,ix)-f16*ds(iz-1,ix))*htd1(iz)
  enddo
  ds(nz,ix)=zero

  do iz=nz-2,1,-1
    ds(iz,ix)=etd1(iz)*ds(iz+1,ix)+ds(iz,ix)
  enddo
enddo

return
end subroutine

!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

subroutine vertint(ds,fs)
 !Finds f by integrating df/dz = d, ensuring f = 0 at the boundaries
 !using 4th-order compact differencing.  Here ds = df/dz and fs = f.

 !Declarations:
implicit none

 !Passed variables:
double precision:: ds(0:nz),fs(0:nz)

 !Local variables:
double precision:: es(nz),esum
integer:: iz

!-------------------------------------------
 !First interpolate ds to a half grid as es:
do iz=1,nz
  es(iz)=f23*(ds(iz-1)+ds(iz))
enddo

es(1)=es(1)*htda(1)
do iz=2,nz
  es(iz)=(es(iz)-f16*es(iz-1))*htda(iz)
enddo

do iz=nz-1,1,-1
  es(iz)=etda(iz)*es(iz+1)+es(iz)
enddo

!-------------------------------------------
 !Next adjust es to ensure f(nz) = 0:
esum=(f1112*(es(1)+es(nz))+sum(es(2:nz-1)))/(dble(nz)-f16)
es=es-esum

 !Integrate:
fs(0)=zero
fs(1)=dz24*(23.d0*es(1)+es(2))
do iz=2,nz-1
  fs(iz)=fs(iz-1)+dz24*(es(iz-1)+22.d0*es(iz)+es(iz+1))
enddo
fs(nz)=zero

return
end subroutine

 !End main program
end program diverge

max_stars_repo_name: matt-frey/epic

max_issues_repo_name: matt-frey/epic

max_forks_repo_name: matt-frey/epic


text:
 Image(roth1.jpg, right, thumbnail, 300)Evan Rothstein was a candidate for ASUCD ASUCD Senate Senate in the Winter 2013 ASUCD Election running on the BEST slate.

Candidate Statement

Hey Aggies! My name is EVAN ROTHSTEIN; Im a second year transfer economics major, intern to senator Maxwell Kappes, and president of Muggle Quidditch at UC Davis.
I am running for senate on the BEST slate because I believe I can make ASUCD incrementally better.

I believe every student has the right to be represented at the table. I believe clubs should be in the senates consciousness. And I believe in keeping UC Davis a leader on environmental issues.

The Associated Student of UC Davis (ASUCD) is the largest club on campus, with every single undergrad paying $123 in membership dues per year.
ASUCDs mission is to enrich the student experience. It runs dozens of services from Unitrans to the CoHo and is the largest employer of undergraduate students on campus.

But ASUCD is only as good as its leaders.

Spring is budget season for ASUCD senate. It is vitally important to have someone who is financially literate and willing to maximize the efficiency of every dollar spent – all 11 million of them. Thats me.

A few years ago, the leadership of Entertainment Council (EC), with the help of Electronic Music for Change (EMC), hosted popular raves every quarter. If elected, my first priority will be to BRING BACK THE RAVES.

My second priority will be to bring modern mobile payments to campus. The modern credit card came into existence in the 1960s, yet wasnt accepted in the CoHo until the mid 1990s. Lets not wait 30 years before adopting the next big thing in payment systems: paying with your smartphone. Its time for UCD to get ahead of the curve on this.
Not only will waving your smartphone at the register make lines move faster, but it will save countless trees by replacing paper receipts with digital ones. Efficiency and environmental win!
I am in contact with the unit director of the CoHo, Darin Schluep, about brining mobile payments there.

Please vote me #1 for senate, fellow BEST candidates #25, and Paul Min & Sergio Cano for executive!

Endorsement Like Statements

I first met Evan Rothstein back when I was both an ASUCD Senator and Manager at the Ben & Jerrys that used to be at the Davis Commons shopping center. Despite him being in high school at the time, he was very community orientated and politically aware. We had many discussions about improving the city of Davis and he even knew a surprising amount of things about what was going on on the UC Davis campus. This was quite a while ago, since I was a senator in 200506 so Evan nows possesses a good deal of institutional memory. I may be a political hasbeen but I believe that Evan will do a wonderful job representing the students of UC Davis. If I had a vote, Id rank him number one.  Rob Roy

20130205 17:18:54 nbsp I heartily endorse this candidate for any office, ASUCD or not Users/StevenDaubert

20130223 18:38:18 nbsp Real talk, if anyone deserved a seat this election, he shouldve. Users/RomanRivilis

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c###outpar.for
      SUBROUTINE OUTPAR
C--------------------------------
C
C     THIS ROUTINE OUTPUTS IONOSPHERIC PARAMETERS
C
      COMMON / CON / D2R, DCL, GAMA, PI, PI2, PIO2, R2D, RZ, VOFL
      COMMON / ES / FS(3,5), HS(5)
      COMMON / GEOG / GYZ(5), RAT(5), GMDIP(5), CLCK(5), ABIY(5),
     1 ARTIC(5), SIGPAT(5), EPSPAT(5)
      COMMON / MFAC / F2M3(5),HPF2(5),ZENANG(5),ZENMAX(5),IEDP,FSECV(3)
      COMMON / OUTLAB / LABEL(11), LAYTYP(5), IEAST, INORTH, ISOUTH,
     A IWEST, LABLI, LABLJ, LABLK
      CHARACTER IEAST*1, ISOUTH*1, INORTH*1, IWEST*1, LABEL*5, LAYTYP*2,
     A LABLI*5, LABLJ*5, LABLK*5
      COMMON / OUTPRT / LINBOT(26), LINBD(14), LINTOP(15), LINTP(14),
     A GRPTYP, JOUT, LINBYP, LINES, LINMAX, LINTYP, LPAGES, NLINE
      COMMON /RON /CLAT(5), CLONG(5), GLAT(5), RD(5), FI(3,5), YI(3,5),
     1HI(3,5), HPRIM(30,5), HTRUE(30,5), FVERT(30,5),KM,KFX, AFAC(30,5),
     2HTR(50,3), FNSQ(50,3)
      COMMON / TIME / IT, GMT, UTIME(24), GMTR, XLMT(24), ITIM, JTX
      COMMON / FILES / LUI, LUO, LU2, LU5, LU6, LU15, LU16, LU20, LU25,
     A LU26, LU35
      CHARACTER IGLAT*1, ISLAT*1, ISLONG*1
      DO 160 K=1,KM
      LINES = LINES + 1
      IF(LINES - LINMAX) 110, 110, 100
C.....CALL SUBROUTINE OUTTOP TO OUTPUT HEADER INFORMATION
  100 CALL OUTTOP
      WRITE(LUO,503)
      WRITE(LUO,500) YI(1,1), HI(1,1), HS(1)
      WRITE(LUO,501)
      LINES = LINTOP(15)
  110 GL = ABS(GLAT(K) * R2D)
      XLAT = ABS(CLAT(K) * R2D)
      XLONG = ABS(CLONG(K) * R2D)
      GYZH = GYZ(K) / 2.
      IGLAT = INORTH
      ISLAT = INORTH
      ISLONG = IEAST
      IF(GLAT(K)) 115, 120, 120
  115 IGLAT = ISOUTH
  120 IF(CLAT(K)) 125, 130, 130
  125 ISLAT = ISOUTH
  130 IF(CLONG(K)) 140, 150, 150
  140 ISLONG = IWEST
  150 WRITE(LUO,502) XLAT, ISLAT, XLONG, ISLONG, CLCK(K), GMT, FI(1,K),
     1 FI(2,K), YI(2,K), HI(2,K), GYZH, FI(3,K), YI(3,K), HI(3,K),
     2 FS(1,K), FS(2,K), FS(3,K), F2M3(K), HPF2(K), RAT(K), ZENANG(K),
     3 ZENMAX(K), GL, IGLAT
  160 CONTINUE
      RETURN
  500 FORMAT(' ',32X,'YE = ',F5.1,5X,'HE = ',F5.1,5X,'HS = ',F5.1,/)
  501 FORMAT(' ',2X,'LAT',3X,'LONG',4X,'LMT',4X,'UT',4X,'E',5X,'F1',4X,
     A 'Y1',5X,'H1',2X,'FH/2',3X,'F2Z',4X,'Y2',5X,'H2',3X,'ES',2X,'MED',
     B 3X,'HI',1X,'M3000',3X,'HPF2',3X,'RAT',4X,'ZEN',2X,'ZMAX',2X,
     C 'MAGL')
  502 FORMAT(2(1X,F5.1,A1),2F6.1,F6.2,2F6.1,F7.1,3F6.1,F7.1,
     1 3F5.1,F6.2,2(F7.1,F6.1),F6.1,A1)
  503 FORMAT(' ')
      END
C--------------------------------

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
      SUBROUTINE DPMJET_OPENINP()
*
*----------------------------------------------------------------------*
* Opens a file with a given unit number
*
*
* IOUNIT: Input unit to be assiged to the file
* FILNAM: Name of the file
*
*----------------------------------------------------------------------*
*
      IMPLICIT NONE
      LOGICAL EXISTS
      INTEGER LNROOT
      INTEGER LNBLNK
      CHARACTER*1000 FILNAM
      CHARACTER*1000 CHROOT
*----------------------------------------------------------------------*
C  input/output channels
      INCLUDE 'inc/dtflka'
*----------------------------------------------------------------------*
      CHROOT=' '
      CALL GETENV('DPMJET_ROOT',CHROOT)
      LNROOT = LNBLNK(CHROOT)
      IF(LNROOT.LE.0) THEN
         FILNAM='dpmjet.dat'
      ELSE
         FILNAM=CHROOT(1:LNROOT)//'/dpmdata/dpmjet.dat'
      ENDIF
      INQUIRE(FILE=FILNAM,EXIST=EXISTS)
      IF(.NOT.EXISTS) THEN
         PRINT*,'***********************************'
         PRINT*,'*           openfile              *'
         PRINT*,'*        ---------------          *'
         PRINT*,'*   File dpmjet.dat not found     *'
         PRINT*,'*         Program STOP            *'
         PRINT*,'*  Check DPMJET_ROOT environment  *'
         PRINT*,'*           variable              *'
         PRINT*,'***********************************'
         STOP
      ENDIF
      LINp = 10
      OPEN (LINp, FILE="dpmjet.inp", STATUS="OLD")
      LDAt = 23
      OPEN (LDAt, FILE=FILNAM, STATUS="OLD")
*
      CALL PHOJET_OPENINP()
*
      RETURN
      END
      SUBROUTINE PHOJET_OPENINP()
*
*----------------------------------------------------------------------*
* Opens a file with a given unit number
*
*
* IOUNIT: Input unit to be assiged to the file
* FILNAM: Name of the file
*
*----------------------------------------------------------------------*
*
      IMPLICIT NONE
      LOGICAL EXISTS
      INTEGER LNROOT
      INTEGER LNBLNK
      CHARACTER*1000 CHROOT
*----------------------------------------------------------------------*
C  input/output channels
      INCLUDE 'inc/poinou'
*----------------------------------------------------------------------*
      CHROOT=' '
      CALL GETENV('DPMJET_ROOT',CHROOT)
      LNROOT = LNBLNK(CHROOT)
      IF(LNROOT.LE.0) THEN
         DATDir='dpmdata/'
         LENDir=8
      ELSE
         DATDir=CHROOT(1:LNROOT)//'/dpmdata/'
         LENDir=LNROOT+9
      ENDIF
      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
c { dg-do run }
        integer(kind=8) foo, bar
        complex(kind=8) c
        data c/(4d10,0)/
        foo = 4d10
        bar = c
        if (foo .ne. bar) call abort
        end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


!-----------------------------------------------------------------------------
        subroutine error(name,msg)
!
!       Generic error handler
!
!-----------------------------------------------------------------------------
        implicit none

        integer loc
        character*80 name, msg

        loc = index(name,'$')-1
        write(*,"(/,'Error in:  ',a)") name(1:loc)
        loc = index(msg,'$')-1
        write(*,"('  ',a)") msg(1:loc)

        call exit(1)

        stop
        end


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 04.02.02.02/02                        *
C  *    TEST TITLE : Setting entries of various kinds      *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

      INTEGER    MAXPM
C define a number greater than the number of initialized polymarker entries
      PARAMETER (MAXPM=300)

C PHIGS parameter types
      INTEGER    MPLBTE, MPMBTE, MTXBTE, MINBTE, MEDBTE, MPAI, MCOLI,
     1           VWTBI
      INTEGER    PMI, MTYPE, COLI, MTYPE2, COLI2

C type of returned value
      INTEGER    PSET,     PREALI
      PARAMETER (PSET = 0, PREALI = 1)

      REAL       MSZSF, MSZSF2
C used to control tests
      INTEGER    MNE, NIE, AII, LISTEL, LISTOT, DEFI(MAXPM), IDUM
      INTEGER    LUI, DPI, CNDE, NNDE

C parameters for <inquire workstation connection and type>
      INTEGER    SPECWT, SPECON

      DATA DEFI / MAXPM*0 /

      CALL INITGL ('04.02.02.02/02')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)
C open workstation
      CALL POPWK  (WKID, CONID, WTYPE)
      CALL PQWKC  (WKID, ERRIND, SPECON, SPECWT)
      CALL CHKINQ ('pqkwc', ERRIND)

      CALL PQWKSL (SPECWT, ERRIND, MPLBTE, MPMBTE, MTXBTE, MINBTE,
     1             MEDBTE, MPAI, MCOLI, VWTBI)
      CALL CHKINQ ('pqwksl', ERRIND)
C maximum number of polymarker bundle table entries
      MNE = MPMBTE

C *** *** *** *** ***    initialized entries   *** *** *** *** ***

C There should be at least 5 initialized polymarker bundles - this
C inquiry gets the first, so it should work.  This statement
C is NOT meant to test the requirement of 5 initialized bundles;
C that is done elsewhere.

      CALL SETMSG ('12', '<Inquire list of polymarker indices> ' //
     1             'should return initialized indices.')
C   nie = number of initialized polymarker entries
C   aii = an initialized index (the first)
      CALL PQEPMI (WKID, 1, ERRIND, NIE, AII)
      IF (ERRIND .NE. 0) THEN
         CALL FAIL
         GOTO 200
      ELSE
         CALL PASS
      ENDIF

C get initialized attributes
      CALL PQPMR (WKID, AII, PSET, ERRIND, MTYPE, MSZSF, COLI)
      CALL CHKINQ ('pqpmr', ERRIND)

C set different values for attributes
C make sure marker type is valid
      IF (MTYPE .LT. 5) THEN
         MTYPE2 = MTYPE + 1
      ELSE
         MTYPE2 = 2
      ENDIF
      MSZSF2 = MSZSF * 0.15
      COLI2  = COLI  + 1
      CALL PSPMR (WKID, AII, MTYPE2, MSZSF2, COLI2)

C get current number of defined entries
      CALL PQEPMI (WKID, 1, ERRIND, CNDE, IDUM)
      CALL CHKINQ ('pqepmi', ERRIND)

      CALL SETMSG ('10 12', 'Re-defining a system-initialized entry ' //
     1             'should not change the number of defined entries.')
      CALL IFPF (CNDE .EQ. NIE)

C check that aii is still among defined entries
      CALL SETMSG ('10 12', 'Initialized bundle entry should be ' //
     1             'reported as defined after being set.')
      DO 150 LISTEL = 1, NIE
         CALL PQEPMI (WKID, LISTEL, ERRIND, IDUM, PMI)
         CALL CHKINQ ('pqepmi', ERRIND)
         IF (PMI .EQ. AII) THEN
            CALL PASS
            GOTO 160
         ENDIF
150   CONTINUE
      CALL FAIL
      GOTO 200
160   CONTINUE

      CALL PQPMR (WKID, AII, PSET, ERRIND, MTYPE, MSZSF, COLI)
      CALL CHKINQ ('pqpmr', ERRIND)
      CALL SETMSG ('11 14', 'New values should be set for ' //
     1             'initialized bundle entry.')
      CALL IFPF (MTYPE .EQ. MTYPE2 .AND.
     1           MSZSF .EQ. MSZSF2 .AND.
     2           COLI  .EQ. COLI2)

C *** *** *** *** ***    undefined entries   *** *** *** *** ***
200   CONTINUE

C get lowest undefined index
C get current number of defined entries
      CALL PQEPMI (WKID, 1, ERRIND, CNDE, IDUM)
      CALL CHKINQ ('pqepmi', ERRIND)
C defi is originally an array of zeros
      DO 250 LISTEL = 1, CNDE
         CALL PQEPMI (WKID, LISTEL, ERRIND, IDUM, PMI)
         CALL CHKINQ ('pqepmi', ERRIND)
C mark as defined any index .le. cnde+1
C This will leave at least one UNmarked.
         IF (PMI .LE. CNDE+1) DEFI (PMI) = 1
250   CONTINUE

C first remaining zero indicates undefined entry
      DO 260 PMI = 1, CNDE+1
         IF (DEFI (PMI) .EQ. 0) THEN
C determine lowest undefined index (not in defined list) = lui
            LUI = PMI
            GOTO 270
         ENDIF
260   CONTINUE
      CALL UNMSG
     1     ('Cannot find undefined index - error in test algorithm.')

C established lui as lowest undefined index
270   CONTINUE
      CALL SETMSG ('12 15 16', '<Inquire polymarker representation> ' //
     1             'should return attributes for bundle #1 when '     //
     2             'inquiring an undefined entry as REALIZED.')

C inquire undefined bundle as REALIZED
      CALL PQPMR  (WKID, LUI, PREALI, ERRIND, MTYPE,  MSZSF,  COLI)
      CALL CHKINQ ('pqpmr', ERRIND)
C inquire bundle #1 as REALIZED
      CALL PQPMR  (WKID, 1,   PREALI, ERRIND, MTYPE2, MSZSF2, COLI2)
      CALL CHKINQ ('pqpmr', ERRIND)
C should have same attributes
      CALL IFPF (MTYPE .EQ. MTYPE2 .AND.
     1           MSZSF .EQ. MSZSF2 .AND.
     2           COLI  .EQ. COLI2)

C  If number of initialized entries >= implementation max then skip test
C  of setting undefined entry - goto redefined_test.
      IF (NIE .GE. MNE) THEN
         CALL INMSG ('Skipping undefined-test because all ' //
     1               'polymarker entries are initialized.')
         GOTO 300
      ENDIF

C Check that program allows enough room
      IF (NIE .GE. MAXPM) THEN
         CALL INMSG ('Skipping undefined-test because number ' //
     1               'of initialized entries exceed program limits.')
         GOTO 300
      ENDIF

C set values for attributes
      MTYPE2 = 2
      MSZSF2 = 0.25
      COLI2  = 4
      CALL PSPMR (WKID, LUI, MTYPE2, MSZSF2, COLI2)

      CALL PQEPMI (WKID, 1, ERRIND, LISTOT, IDUM)
      CALL CHKINQ ('pqepmi', ERRIND)

      CALL SETMSG ('10 12', 'Defining an undefined entry should ' //
     1             'increment the number of defined entries.')
      CALL IFPF (LISTOT .EQ. CNDE + 1)

C check that lui is now among defined entries
      CALL SETMSG ('10 12', 'Undefined bundle entry should be ' //
     1             'reported as defined after set.')
      DO 280 LISTEL = 1, LISTOT
         CALL PQEPMI (WKID, LISTEL, ERRIND, IDUM, PMI)
         CALL CHKINQ ('pqepmi', ERRIND)
         IF (PMI .EQ. LUI) THEN
            CALL PASS
            GOTO 290
         ENDIF
280   CONTINUE
      CALL FAIL
      GOTO 300
290   CONTINUE

      CALL PQPMR (WKID, LUI, PSET, ERRIND, MTYPE, MSZSF, COLI)
      CALL CHKINQ ('pqpmr', ERRIND)
      CALL SETMSG ('11 14', 'New values should be set for ' //
     1             'undefined bundle entry.')
      CALL IFPF (MTYPE .EQ. MTYPE2 .AND.
     1           MSZSF .EQ. MSZSF2 .AND.
     2           COLI  .EQ. COLI2)

C *** *** *** *** ***    redefined entries   *** *** *** *** ***
300   CONTINUE

C get current number defined entries
      CALL PQEPMI (WKID, 1, ERRIND, CNDE, IDUM)
      CALL CHKINQ ('pqepmi', ERRIND)

C let dpi = a defined polymarker index from middle of list
      CALL PQEPMI (WKID, (CNDE+1)/2, ERRIND, IDUM, DPI)
      CALL CHKINQ ('pqepmi', ERRIND)

C set to first set of values
      MTYPE = 4
      MSZSF = 0.35
      COLI  = 3
      CALL PSPMR (WKID, DPI, MTYPE, MSZSF, COLI)
C re-set to different values
      MTYPE2 = 3
      MSZSF2 = 0.40
      COLI2  = 2
      CALL PSPMR (WKID, DPI, MTYPE2, MSZSF2, COLI2)

C get list size - nnde = new number of defined entries
      CALL PQEPMI (WKID, 1, ERRIND, NNDE, IDUM)
      CALL CHKINQ ('pqepmi', ERRIND)

      CALL SETMSG ('10 12', 'Re-defining a user-defined entry ' //
     1             'should not change number of defined entries.')
      CALL IFPF (NNDE .EQ. CNDE)

C check that dpi is among defined entries
      CALL SETMSG ('10 12', 'Redefined bundle entry should be ' //
     1             'reported as defined after set.')
      DO 320 LISTEL = 1, NNDE
         CALL PQEPMI (WKID, LISTEL, ERRIND, IDUM, PMI)
         CALL CHKINQ ('pqepmi', ERRIND)
         IF (PMI .EQ. DPI) THEN
            CALL PASS
            GOTO 330
         ENDIF
320   CONTINUE
      CALL FAIL
      GOTO 399
330   CONTINUE

      CALL PQPMR (WKID, DPI, PSET, ERRIND, MTYPE, MSZSF, COLI)
      CALL CHKINQ ('pqpmr', ERRIND)
      CALL SETMSG ('11 14', 'New values should be set for ' //
     1             'redefined bundle entry.')
      CALL IFPF (MTYPE .EQ. MTYPE2 .AND.
     1           MSZSF .EQ. MSZSF2 .AND.
     2           COLI  .EQ. COLI2)

399   CONTINUE

C *** *** *** *** ***  access all defined entries  *** *** *** *** ***

      CALL SETMSG ('12 13', '<Inquire list of polymarker indices> '  //
     1             'and <inquire polymarker representation> should ' //
     2             'detect all defined entries.')
C get listel-th element of defined polymarker indices
      DO 410 LISTEL = 1, NNDE
         CALL PQEPMI (WKID, LISTEL, ERRIND, IDUM, PMI)
         IF (ERRIND .NE. 0) THEN
            CALL FAIL
            GOTO 415
         ENDIF
         CALL PQPMR (WKID, PMI, PSET, ERRIND, MTYPE, MSZSF, COLI)
         IF (ERRIND .NE. 0) THEN
            CALL FAIL
            GOTO 415
         ENDIF
410   CONTINUE
      CALL PASS
415   CONTINUE

666   CONTINUE
      CALL ENDIT

      END


text:
program quad_test
use genquad, only: get_quadrature_rule
implicit none
  logical:: diagnostics_on = .true.
  integer:: N_scheme = 4
  integer:: N_integrand = 4
  integer:: N_res = 59
  real:: vmax = 4.0
  real:: pi

  real,dimension(:,:,:), allocatable:: errors, errors_gs2
  real,dimension(:,:), allocatable:: wgt_table_test, analytic_result
  real,dimension(:), allocatable:: v_table_test, abscissae, int_wgts, integrand
  real,dimension(:), allocatable:: gs2_wgts, gs2_vgrid, gs2_integrand, uniform_vgrid
  real,dimension(:), allocatable:: simpson_wgts, vgrid_overresolved, wgts_overresolved
  real:: delta_v , numerical_result, zero = 0.0, gs2_result, maxerr_v, maxerr_w
  integer::  iv,ischeme,ires,iint, Ntot
  integer, dimension(:),allocatable:: Nv_res
  real,external:: w1,w2,w3,w4
  real,dimension(10):: hermitenodes,hermiteweights

  allocate(Nv_res(N_res))

  pi = 3.14159265358979323846264338327950288419716939938

  do ires = 1,N_res
    Nv_res(ires) =ires+1
  end do

  open(unit=22,file="results.out", status="replace")
  open(unit=23,file="wgt1.plt", status="replace")
  open(unit=24,file="wgt2.plt", status="replace")
  open(unit=25,file="wgt3.plt", status="replace")
  open(unit=26,file="wgt4.plt", status="replace")

  allocate(analytic_result(N_scheme,N_integrand))
  allocate(abscissae(1:maxval(Nv_res(:))))
  allocate(int_wgts(1:maxval(Nv_res(:))))
  allocate(integrand(1:maxval(Nv_res(:))))
  allocate(gs2_integrand(1:maxval(Nv_res(:))))
  allocate(gs2_wgts(1:maxval(Nv_res(:))))
  allocate(gs2_vgrid(1:maxval(Nv_res(:))))
  allocate(errors(N_scheme,N_integrand,N_res))
  allocate(errors_gs2(N_scheme,N_integrand,N_res))

  ! Define analytic integrals
  ! first index is weight scheme
  ! second index is integrand
  ! Calculated from mathematica
  analytic_result(1,1) = 23.48501352449811
  analytic_result(1,2) = 0.03246362468013194
  analytic_result(1,3) = 0.801369616941368
  analytic_result(1,4) = 27.61348974269745

  analytic_result(2,1) = 619.0
  analytic_result(2,2) = 0.05882352941176441
  analytic_result(2,3) = 0.949590264857363
  analytic_result(2,4) = 12.64928118705395

  analytic_result(3,1) = 59.22503999731359
  analytic_result(3,2) = 3.39405260298183
  analytic_result(3,3) = 3.604857301284351
  analytic_result(3,4) = 58.20370297469304

  analytic_result(4,1) = 302.4775729802816
  analytic_result(4,2) = 0.02621488484812071
  analytic_result(4,3) = 1.267867002167052
  analytic_result(4,4) = 16.38784893315581

  do ischeme = 1,4
    write(22,*) "*****************************"
    write(22,'(A27,I2)') " Quadrature weight Scheme ", ischeme
    write(22,*) "*****************************"
    do ires = 1,N_res
      write(22,'(A5,I2)') "Nv = ", Nv_res(ires)
!      write(*,'(A5,I2)') "Nv = ", Nv_res(ires)

      if (ischeme .EQ. 1 ) then
        call get_quadrature_rule(w1,Nv_res(ires),0.0,1.0,abscissae,int_wgts,.true.,.true.)
      else if (ischeme .EQ. 2) then
        call get_quadrature_rule(w2,Nv_res(ires),0.0,1.0,abscissae,int_wgts,.false.,.true.)
      else if (ischeme .EQ. 3) then
        call get_quadrature_rule(w3,Nv_res(ires),-1.0,4.0,abscissae,int_wgts,.false.,.false.)
      else if (ischeme .EQ. 4) then
        call get_quadrature_rule(w4,Nv_res(ires),0.0,4.0,abscissae,int_wgts,.false.,.false.)
      end if

      write(22,*) "Calculated quadrature abscissae and weights:"
      do iv = 1,Nv_res(ires)
         write(22,*) abscissae(iv), int_wgts(iv)
      end do
      do iint = 1,N_integrand
        call set_integrand(iint,Nv_res(ires),abscissae(1:Nv_res(ires)),integrand(1:Nv_res(ires)))
        call perform_integral(Nv_res(ires),integrand(1:Nv_res(ires)),int_wgts(1:Nv_res(ires)),numerical_result)

        write(22,*) "Analytic result:  ", analytic_result(ischeme,iint)
        write(22,*) "Numerical result: ", numerical_result

        write(22,*) "---------------------"
        write(22,*) "Error = ", abs( numerical_result - analytic_result(ischeme,iint))/abs(numerical_result)
        write(22,*) "---------------------"
        write(22,*) " "
        errors(ischeme,iint,ires) = abs(numerical_result - analytic_result(ischeme,iint))/abs(analytic_result(ischeme,iint))
      end do
    end do
  end do

  close(22)

     do ires = 1,N_res
        write(23,*) Nv_res(ires), (errors(1,iint,ires),iint=1,N_integrand),(errors_gs2(1,iint,ires),iint=1,N_integrand)
        write(24,*) Nv_res(ires), (errors(2,iint,ires),iint=1,N_integrand),(errors_gs2(2,iint,ires),iint=1,N_integrand)
        write(25,*) Nv_res(ires), (errors(3,iint,ires),iint=1,N_integrand),(errors_gs2(3,iint,ires),iint=1,N_integrand)
        write(26,*) Nv_res(ires), (errors(4,iint,ires),iint=1,N_integrand),(errors_gs2(4,iint,ires),iint=1,N_integrand)
     end do

  close(23)
  close(24)
  close(25)
  close(26)

  open(unit=27,file="hermitecomp.dat", status="replace")
  call get_quadrature_rule(w1,20,0.0,1.0,abscissae,int_wgts,.true.,.true.)

  ! Actual Hermite quadrature rule for N=20, from Abramowitz and Stegun
  hermitenodes = (/ 0.2453407083009, 0.7374737285454, 1.2340762153953, 1.7385377121166, 2.2549740020893, 2.7888060584281, &
       3.3478545673832, 3.9447640401156, 4.6036824495507, 5.3874808900112  /)

  hermiteweights = (/4.622436696006e-1, 2.866755053628e-1, 1.090172060200e-1, 2.481052088746e-2, 3.243773342238e-3, &
      2.283386360163e-4, 7.802556478532e-6, 1.086069370769e-7, 4.399340992273e-10, 2.229393645534e-13 /)

  maxerr_v = 0.0
  maxerr_w = 0.0
  do iv=1,10
    write(27,*) iv, " & ", abscissae(iv), " & ", -hermitenodes(10-iv+1), " & ", int_wgts(iv), " & ", hermiteweights(10-iv+1) &
      , " \\ "
    maxerr_v = max(maxerr_v, abs( (abscissae(iv)+hermitenodes(10-iv+1))/hermitenodes(10-iv+1)))
    maxerr_w = max(maxerr_w, abs( (int_wgts(iv)-hermiteweights(10-iv+1))/hermiteweights(10-iv+1)))
  end do
  do iv=11,20
    write(27,*) iv, " & ", abscissae(iv), " & ", hermitenodes(iv-10), " & ", int_wgts(iv), " & ", hermiteweights(iv-10)," \\ "
    maxerr_v = max(maxerr_v, abs( (abscissae(iv)-hermitenodes(iv-10))/hermitenodes(iv-10)))
    maxerr_w = max(maxerr_w, abs( (int_wgts(iv)-hermiteweights(iv-10))/hermiteweights(iv-10)))
  end do
  write(*,*) "Max error in Hermite abscissae = ", maxerr_v
  write(*,*) "Max error in Hermite weights = ", maxerr_w

  close(27)

contains

subroutine set_integrand(opt,N,v,f)
implicit none
  integer,intent(in):: opt, N
  real,dimension(1:), intent(in):: v
  real,dimension(1:), intent(inout):: f
  integer:: iv
!  real,external:: g1,g2,g3,g4

  do iv =1,N
    select case (opt)
      case(1)
        f(iv) = g1(v(iv))
      case(2)
        f(iv) = g2(v(iv))
      case(3)
        f(iv) = g3(v(iv))
      case(4)
        f(iv) = g4(v(iv))
    end select
  end do

end subroutine set_integrand

  ! Integrand 1:
  ! Polynomial
real function g1(v)
implicit none
  real,intent(in):: v
  g1 = 9.0 + v*(-10.0 + v*(7.0 + v*(-3.0 + v*(1.0 + v*5.0))))
  return
end function

  ! Integrand 2:
  ! Triognometric
real function g2(v)
implicit none
  real,intent(in):: v
  g2 = cos(4.0*v)
  return
end function

  ! Integrand 2:
  ! Hyperbolic
real function g3(v)
implicit none
  real,intent(in):: v
  g3 = tanh(2.0*v+1.0)
  return
end function
  ! Integrand 4:
  ! Exponential
real function g4(v)
implicit none
  real,intent(in):: v
  g4 = 10.0*exp(-v**4+2.0*v**2)
  return
end function

subroutine perform_integral(N,f,w,answer)
implicit none
   real,dimension(1:), intent(in):: f,w
   real, intent(out):: answer
   integer, intent(in):: N
   real:: sum
   integer:: iv

   sum = 0.0
   do iv= 1,N
     sum = sum + w(iv)*f(iv)
   end do

   answer = sum

end subroutine perform_integral

end program

  ! Scheme 1: w(v) =  exp(-v)
  ! Exponential
real function w1(v)
implicit none
  real,intent(in):: v
  w1 = exp(-v**2)
  return
end function

real function w2(v)
implicit none
  real,intent(in):: v
  w2 = exp(-v)
  return
end function

real function w3(v)
implicit none
  real,intent(in):: v
  w3 = 1.0/(0.1+10.0*abs(v)**3)
!  w3 = 1.0
  return
end function

real function w4(v)
implicit none
  real,intent(in):: v
  w4 = 1.0/(1.0+v**2)
  return
end function

max_stars_repo_name: gjwilkie/genquad

max_issues_repo_name: gjwilkie/genquad

max_forks_repo_name: gjwilkie/genquad