dslugm.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	dslugm (N, B, X, NELT, IA, JA, A, ISYM, NSAVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW)
subroutine	dslui (N, B, X, NELT, IA, JA, A, ISYM, RWORK, IWORK)
subroutine	dslui2 (N, B, X, IL, JL, L, DINV, IU, JU, U)
subroutine	dsmv (N, X, Y, NELT, IA, JA, A, ISYM)
subroutine	dchkw (NAME, LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR)
subroutine	ds2y (N, NELT, IA, JA, A, ISYM)
subroutine	dsilus (N, NELT, IA, JA, A, ISYM, NL, IL, JL, L, DINV, NU, IU, JU, U, NROW, NCOL)
subroutine	qs2i1d (IA, JA, A, N, KFLAG)

Function/Subroutine Documentation

dchkw()

subroutine dchkw	(	character, dimension(*)	NAME,
		integer	LOCIW,
		integer	LENIW,
		integer	LOCW,
		integer	LENW,
		integer	IERR,
		integer	ITER,
		double precision	ERR
	)

C***BEGIN PROLOGUE  DCHKW
C***SUBSIDIARY
C***PURPOSE  SLAP WORK/IWORK Array Bounds Checker.
C            This routine checks the work array lengths and interfaces
C            to the SLATEC error handler if a problem is found.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  R2
C***TYPE      DOUBLE PRECISION (SCHKW-S, DCHKW-D)
C***KEYWORDS  ERROR CHECKING, SLAP, WORKSPACE CHECKING
C***AUTHOR  Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     CHARACTER*(*) NAME
C     INTEGER LOCIW, LENIW, LOCW, LENW, IERR, ITER
C     DOUBLE PRECISION ERR
C
C     CALL DCHKW( NAME, LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR )
C
C *Arguments:
C NAME   :IN       Character*(*).
C         Name of the calling routine.  This is used in the output
C         message, if an error is detected.
C LOCIW  :IN       Integer.
C         Location of the first free element in the integer workspace
C         array.
C LENIW  :IN       Integer.
C         Length of the integer workspace array.
C LOCW   :IN       Integer.
C         Location of the first free element in the double precision
C         workspace array.
C LENRW  :IN       Integer.
C         Length of the double precision workspace array.
C IERR   :OUT      Integer.
C         Return error flag.
C               IERR = 0 => All went well.
C               IERR = 1 => Insufficient storage allocated for
C                           WORK or IWORK.
C ITER   :OUT      Integer.
C         Set to zero on return.
C ERR    :OUT      Double Precision.
C         Set to the smallest positive magnitude if all went well.
C         Set to a very large number if an error is detected.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   880225  DATE WRITTEN
C   881213  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   900805  Changed XERRWV calls to calls to XERMSG.  (RWC)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   910502  Corrected XERMSG calls to satisfy Section 6.2.2 of ANSI
C           X3.9-1978.  (FNF)
C   910506  Made subsidiary.  (FNF)
C   920511  Added complete declaration section.  (WRB)
C   921015  Added code to initialize ITER and ERR when IERR=0.  (FNF)
C***END PROLOGUE  DCHKW
C     .. Scalar Arguments ..
      DOUBLE PRECISION err
      INTEGER ierr, iter, leniw, lenw, lociw, locw
      CHARACTER name*(*)
C     .. Local Scalars ..
      CHARACTER xern1*8, xern2*8, xernam*8
C     .. External Functions ..
      DOUBLE PRECISION d1mach
      EXTERNAL d1mach
C     .. External Subroutines ..
      EXTERNAL xermsg
C***FIRST EXECUTABLE STATEMENT  DCHKW
C
C         Check the Integer workspace situation.
C
      ierr = 0
      iter = 0
      err = d1mach(1)
      IF( lociw.GT.leniw ) THEN
         ierr = 1
         err = d1mach(2)
         xernam = name
         WRITE (xern1, '(I8)') lociw
         WRITE (xern2, '(I8)') leniw
         CALL xermsg ('SLATEC', 'DCHKW',
     $      'In ' // xernam // ', INTEGER work array too short.  ' //
     $      'IWORK needs ' // xern1 // '; have allocated ' // xern2,
     $      1, 1)
      ENDIF
C
C         Check the Double Precision workspace situation.
      IF( locw.GT.lenw ) THEN
         ierr = 1
         err = d1mach(2)
         xernam = name
         WRITE (xern1, '(I8)') locw
         WRITE (xern2, '(I8)') lenw
         CALL xermsg ('SLATEC', 'DCHKW',
     $      'In ' // xernam // ', DOUBLE PRECISION work array too ' //
     $      'short.  RWORK needs ' // xern1 // '; have allocated ' //
     $      xern2, 1, 1)
      ENDIF
      RETURN
C------------- LAST LINE OF DCHKW FOLLOWS ----------------------------

ds2y()

subroutine ds2y	(	integer	N,
		integer	NELT,
		integer, dimension(nelt)	IA,
		integer, dimension(nelt)	JA,
		double precision, dimension(nelt)	A,
		integer	ISYM
	)

C***BEGIN PROLOGUE  DS2Y
C***PURPOSE  SLAP Triad to SLAP Column Format Converter.
C            Routine to convert from the SLAP Triad to SLAP Column
C            format.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D1B9
C***TYPE      DOUBLE PRECISION (SS2Y-S, DS2Y-D)
C***KEYWORDS  LINEAR SYSTEM, SLAP SPARSE
C***AUTHOR  Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM
C     DOUBLE PRECISION A(NELT)
C
C     CALL DS2Y( N, NELT, IA, JA, A, ISYM )
C
C *Arguments:
C N      :IN       Integer
C         Order of the Matrix.
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    Double Precision 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 used, this format is
C         translated to the SLAP Column format by this routine.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all non-zero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the lower
C         triangle of the matrix is stored.
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.  If the SLAP Triad format is give
C       as input then this routine transforms it into SLAP Column
C       format.  The way this routine tells which format is given as
C       input is to look at JA(N+1).  If JA(N+1) = NELT+1 then we
C       have the SLAP Column format.  If that equality does not hold
C       then it is assumed that the IA, JA, A arrays contain the
C       SLAP Triad format.
C
C       =================== S L A P Triad format ===================
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
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***REFERENCES  (NONE)
C***ROUTINES CALLED  QS2I1D
C***REVISION HISTORY  (YYMMDD)
C   871119  DATE WRITTEN
C   881213  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   910502  Corrected C***FIRST EXECUTABLE STATEMENT line.  (FNF)
C   920511  Added complete declaration section.  (WRB)
C   930701  Updated CATEGORY section.  (FNF, WRB)
C***END PROLOGUE  DS2Y
C     .. Scalar Arguments ..
      INTEGER isym, n, nelt
C     .. Array Arguments ..
      DOUBLE PRECISION a(nelt)
      INTEGER ia(nelt), ja(nelt)
C     .. Local Scalars ..
      DOUBLE PRECISION temp
      INTEGER i, ibgn, icol, iend, itemp, j
C     .. External Subroutines ..
      EXTERNAL qs2i1d
C***FIRST EXECUTABLE STATEMENT  DS2Y
C
C         Check to see if the (IA,JA,A) arrays are in SLAP Column
C         format.  If it's not then transform from SLAP Triad.
C
      IF( ja(n+1).EQ.nelt+1 ) RETURN
C
C         Sort into ascending order by COLUMN (on the ja array).
C         This will line up the columns.
C
      CALL qs2i1d( ja, ia, a, nelt, 1 )
C
C         Loop over each column to see where the column indices change
C         in the column index array ja.  This marks the beginning of the
C         next column.
C
CVD$R NOVECTOR
      ja(1) = 1
      DO 20 icol = 1, n-1
         DO 10 j = ja(icol)+1, nelt
            IF( ja(j).NE.icol ) THEN
               ja(icol+1) = j
               GOTO 20
            ENDIF
 10      CONTINUE
 20   CONTINUE
      ja(n+1) = nelt+1
C
C         Mark the n+2 element so that future calls to a SLAP routine
C         utilizing the YSMP-Column storage format will be able to tell.
C
      ja(n+2) = 0
C
C         Now loop through the IA array making sure that the diagonal
C         matrix element appears first in the column.  Then sort the
C         rest of the column in ascending order.
C
      DO 70 icol = 1, n
         ibgn = ja(icol)
         iend = ja(icol+1)-1
         DO 30 i = ibgn, iend
            IF( ia(i).EQ.icol ) THEN
C
C              Swap the diagonal element with the first element in the
C              column.
C
               itemp = ia(i)
               ia(i) = ia(ibgn)
               ia(ibgn) = itemp
               temp = a(i)
               a(i) = a(ibgn)
               a(ibgn) = temp
               GOTO 40
            ENDIF
 30      CONTINUE
 40      ibgn = ibgn + 1
         IF( ibgn.LT.iend ) THEN
            DO 60 i = ibgn, iend
               DO 50 j = i+1, iend
                  IF( ia(i).GT.ia(j) ) THEN
                     itemp = ia(i)
                     ia(i) = ia(j)
                     ia(j) = itemp
                     temp = a(i)
                     a(i) = a(j)
                     a(j) = temp
                  ENDIF
 50            CONTINUE
 60         CONTINUE
         ENDIF
 70   CONTINUE
      RETURN
C------------- LAST LINE OF DS2Y FOLLOWS ----------------------------

dsilus()

subroutine dsilus	(	integer	N,
		integer	NELT,
		integer, dimension(nelt)	IA,
		integer, dimension(nelt)	JA,
		double precision, dimension(nelt)	A,
		integer	ISYM,
		integer	NL,
		integer, dimension(nl)	IL,
		integer, dimension(nl)	JL,
		double precision, dimension(nl)	L,
		double precision, dimension(n)	DINV,
		integer	NU,
		integer, dimension(nu)	IU,
		integer, dimension(nu)	JU,
		double precision, dimension(nu)	U,
		integer, dimension(n)	NROW,
		integer, dimension(n)	NCOL
	)

C***BEGIN PROLOGUE  DSILUS
C***PURPOSE  Incomplete LU Decomposition Preconditioner SLAP Set Up.
C            Routine to generate the incomplete LDU decomposition of a
C            matrix.  The unit lower triangular factor L is stored by
C            rows and the unit upper triangular factor U is stored by
C            columns.  The inverse of the diagonal matrix D is stored.
C            No fill in is allowed.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D2E
C***TYPE      DOUBLE PRECISION (SSILUS-S, DSILUS-D)
C***KEYWORDS  INCOMPLETE LU FACTORIZATION, ITERATIVE PRECONDITION,
C             NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
C***AUTHOR  Greenbaum, Anne, (Courant Institute)
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM
C     INTEGER NL, IL(NL), JL(NL), NU, IU(NU), JU(NU)
C     INTEGER NROW(N), NCOL(N)
C     DOUBLE PRECISION A(NELT), L(NL), DINV(N), U(NU)
C
C     CALL DSILUS( N, NELT, IA, JA, A, ISYM, NL, IL, JL, L,
C    $    DINV, NU, IU, JU, U, NROW, NCOL )
C
C *Arguments:
C N      :IN       Integer
C         Order of the Matrix.
C NELT   :IN       Integer.
C         Number of elements in arrays IA, JA, and A.
C IA     :IN       Integer IA(NELT).
C JA     :IN       Integer JA(NELT).
C A      :IN       Double Precision A(NELT).
C         These arrays should hold the matrix A in the SLAP Column
C         format.  See "Description", below.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all non-zero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the lower
C         triangle of the matrix is stored.
C NL     :OUT      Integer.
C         Number of non-zeros in the L array.
C IL     :OUT      Integer IL(NL).
C JL     :OUT      Integer JL(NL).
C L      :OUT      Double Precision L(NL).
C         IL, JL, L  contain the unit lower triangular factor of  the
C         incomplete decomposition  of some  matrix stored  in   SLAP
C         Row format.     The   Diagonal  of ones  *IS*  stored.  See
C         "DESCRIPTION", below for more details about the SLAP format.
C NU     :OUT      Integer.
C         Number of non-zeros in the U array.
C IU     :OUT      Integer IU(NU).
C JU     :OUT      Integer JU(NU).
C U      :OUT      Double Precision     U(NU).
C         IU, JU, U contain   the unit upper triangular factor of the
C         incomplete  decomposition    of some matrix  stored in SLAP
C         Column  format.   The Diagonal of ones   *IS*  stored.  See
C         "Description", below  for  more  details  about  the   SLAP
C         format.
C NROW   :WORK     Integer NROW(N).
C         NROW(I) is the number of non-zero elements in the I-th row
C         of L.
C NCOL   :WORK     Integer NCOL(N).
C         NCOL(I) is the number of non-zero elements in the I-th
C         column of U.
C
C *Description
C       IL, JL, L should contain the unit  lower triangular factor of
C       the incomplete decomposition of the A matrix  stored in SLAP
C       Row format.  IU, JU, U should contain  the unit upper factor
C       of the  incomplete decomposition of  the A matrix  stored in
C       SLAP Column format This ILU factorization can be computed by
C       the DSILUS routine. The diagonals (which are all one's) are
C       stored.
C
C       =================== S L A P Column format ==================
C
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       ==================== S L A P Row format ====================
C
C       This routine requires  that the matrix A  be  stored  in the
C       SLAP  Row format.   In this format  the non-zeros are stored
C       counting across  rows (except for the diagonal  entry, which
C       must  appear first  in each  "row")  and  are stored  in the
C       double precision  array A.  In other words, for each row  in
C       the matrix  put the diagonal  entry in A.   Then put in  the
C       other  non-zero elements  going across  the row  (except the
C       diagonal) in order.  The JA array holds the column index for
C       each non-zero.  The IA array holds the offsets  into the JA,
C       A  arrays  for  the   beginning  of  each  row.    That  is,
C       JA(IA(IROW)),A(IA(IROW)) are the first elements of the IROW-
C       th row in  JA and A,  and  JA(IA(IROW+1)-1), A(IA(IROW+1)-1)
C       are  the last elements  of the  IROW-th row.   Note  that we
C       always have  IA(N+1) = NELT+1, where N is the number of rows
C       in the matrix  and  NELT is the  number of non-zeros  in the
C       matrix.
C
C       Here is an example of the SLAP Row storage format for a  5x5
C       Matrix (in the A and JA arrays '|' denotes the end of a row):
C
C           5x5 Matrix         SLAP Row 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 12 15 | 22 21 | 33 35 | 44 | 55 51 53
C       |21 22  0  0  0|  JA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  IA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C***SEE ALSO  SILUR
C***REFERENCES  1. Gene Golub and Charles Van Loan, Matrix Computations,
C                  Johns Hopkins University Press, Baltimore, Maryland,
C                  1983.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   890404  DATE WRITTEN
C   890404  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   920511  Added complete declaration section.  (WRB)
C   920929  Corrected format of reference.  (FNF)
C   930701  Updated CATEGORY section.  (FNF, WRB)
C***END PROLOGUE  DSILUS
C     .. Scalar Arguments ..
      INTEGER isym, n, nelt, nl, nu
C     .. Array Arguments ..
      DOUBLE PRECISION a(nelt), dinv(n), l(nl), u(nu)
      INTEGER ia(nelt), il(nl), iu(nu), ja(nelt), jl(nl), ju(nu),
     +        ncol(n), nrow(n)
C     .. Local Scalars ..
      DOUBLE PRECISION temp
      INTEGER i, ibgn, icol, iend, indx, indx1, indx2, indxc1, indxc2,
     +        indxr1, indxr2, irow, itemp, j, jbgn, jend, jtemp, k, kc,
     +        kr
C***FIRST EXECUTABLE STATEMENT  DSILUS
C
C         Count number of elements in each row of the lower triangle.
C
      DO 10 i=1,n
         nrow(i) = 0
         ncol(i) = 0
 10   CONTINUE
CVD$R NOCONCUR
CVD$R NOVECTOR
      DO 30 icol = 1, n
         jbgn = ja(icol)+1
         jend = ja(icol+1)-1
         IF( jbgn.LE.jend ) THEN
            DO 20 j = jbgn, jend
               IF( ia(j).LT.icol ) THEN
                  ncol(icol) = ncol(icol) + 1
               ELSE
                  nrow(ia(j)) = nrow(ia(j)) + 1
                  IF( isym.NE.0 ) ncol(ia(j)) = ncol(ia(j)) + 1
               ENDIF
 20         CONTINUE
         ENDIF
 30   CONTINUE
      ju(1) = 1
      il(1) = 1
      DO 40 icol = 1, n
         il(icol+1) = il(icol) + nrow(icol)
         ju(icol+1) = ju(icol) + ncol(icol)
         nrow(icol) = il(icol)
         ncol(icol) = ju(icol)
 40   CONTINUE
C
C         Copy the matrix A into the L and U structures.
      DO 60 icol = 1, n
         dinv(icol) = a(ja(icol))
         jbgn = ja(icol)+1
         jend = ja(icol+1)-1
         IF( jbgn.LE.jend ) THEN
            DO 50 j = jbgn, jend
               irow = ia(j)
               IF( irow.LT.icol ) THEN
C         Part of the upper triangle.
                  iu(ncol(icol)) = irow
                  u(ncol(icol)) = a(j)
                  ncol(icol) = ncol(icol) + 1
               ELSE
C         Part of the lower triangle (stored by row).
                  jl(nrow(irow)) = icol
                  l(nrow(irow)) = a(j)
                  nrow(irow) = nrow(irow) + 1
                  IF( isym.NE.0 ) THEN
C         Symmetric...Copy lower triangle into upper triangle as well.
                     iu(ncol(irow)) = icol
                     u(ncol(irow)) = a(j)
                     ncol(irow) = ncol(irow) + 1
                  ENDIF
               ENDIF
 50         CONTINUE
         ENDIF
 60   CONTINUE
C
C         Sort the rows of L and the columns of U.
      DO 110 k = 2, n
         jbgn = ju(k)
         jend = ju(k+1)-1
         IF( jbgn.LT.jend ) THEN
            DO 80 j = jbgn, jend-1
               DO 70 i = j+1, jend
                  IF( iu(j).GT.iu(i) ) THEN
                     itemp = iu(j)
                     iu(j) = iu(i)
                     iu(i) = itemp
                     temp = u(j)
                     u(j) = u(i)
                     u(i) = temp
                  ENDIF
 70            CONTINUE
 80         CONTINUE
         ENDIF
         ibgn = il(k)
         iend = il(k+1)-1
         IF( ibgn.LT.iend ) THEN
            DO 100 i = ibgn, iend-1
               DO 90 j = i+1, iend
                  IF( jl(i).GT.jl(j) ) THEN
                     jtemp = ju(i)
                     ju(i) = ju(j)
                     ju(j) = jtemp
                     temp = l(i)
                     l(i) = l(j)
                     l(j) = temp
                  ENDIF
 90            CONTINUE
 100        CONTINUE
         ENDIF
 110  CONTINUE
C
C         Perform the incomplete LDU decomposition.
      DO 300 i=2,n
C
C           I-th row of L
         indx1 = il(i)
         indx2 = il(i+1) - 1
         IF(indx1 .GT. indx2) GO TO 200
         DO 190 indx=indx1,indx2
            IF(indx .EQ. indx1) GO TO 180
            indxr1 = indx1
            indxr2 = indx - 1
            indxc1 = ju(jl(indx))
            indxc2 = ju(jl(indx)+1) - 1
            IF(indxc1 .GT. indxc2) GO TO 180
 160        kr = jl(indxr1)
 170        kc = iu(indxc1)
            IF(kr .GT. kc) THEN
               indxc1 = indxc1 + 1
               IF(indxc1 .LE. indxc2) GO TO 170
            ELSEIF(kr .LT. kc) THEN
               indxr1 = indxr1 + 1
               IF(indxr1 .LE. indxr2) GO TO 160
            ELSEIF(kr .EQ. kc) THEN
               l(indx) = l(indx) - l(indxr1)*dinv(kc)*u(indxc1)
               indxr1 = indxr1 + 1
               indxc1 = indxc1 + 1
               IF(indxr1 .LE. indxr2 .AND. indxc1 .LE. indxc2) GO TO 160
            ENDIF
 180        l(indx) = l(indx)/dinv(jl(indx))
 190     CONTINUE
C
C         I-th column of U
 200     indx1 = ju(i)
         indx2 = ju(i+1) - 1
         IF(indx1 .GT. indx2) GO TO 260
         DO 250 indx=indx1,indx2
            IF(indx .EQ. indx1) GO TO 240
            indxc1 = indx1
            indxc2 = indx - 1
            indxr1 = il(iu(indx))
            indxr2 = il(iu(indx)+1) - 1
            IF(indxr1 .GT. indxr2) GO TO 240
 210        kr = jl(indxr1)
 220        kc = iu(indxc1)
            IF(kr .GT. kc) THEN
               indxc1 = indxc1 + 1
               IF(indxc1 .LE. indxc2) GO TO 220
            ELSEIF(kr .LT. kc) THEN
               indxr1 = indxr1 + 1
               IF(indxr1 .LE. indxr2) GO TO 210
            ELSEIF(kr .EQ. kc) THEN
               u(indx) = u(indx) - l(indxr1)*dinv(kc)*u(indxc1)
               indxr1 = indxr1 + 1
               indxc1 = indxc1 + 1
               IF(indxr1 .LE. indxr2 .AND. indxc1 .LE. indxc2) GO TO 210
            ENDIF
 240        u(indx) = u(indx)/dinv(iu(indx))
 250     CONTINUE
C
C         I-th diagonal element
 260     indxr1 = il(i)
         indxr2 = il(i+1) - 1
         IF(indxr1 .GT. indxr2) GO TO 300
         indxc1 = ju(i)
         indxc2 = ju(i+1) - 1
         IF(indxc1 .GT. indxc2) GO TO 300
 270     kr = jl(indxr1)
 280     kc = iu(indxc1)
         IF(kr .GT. kc) THEN
            indxc1 = indxc1 + 1
            IF(indxc1 .LE. indxc2) GO TO 280
         ELSEIF(kr .LT. kc) THEN
            indxr1 = indxr1 + 1
            IF(indxr1 .LE. indxr2) GO TO 270
         ELSEIF(kr .EQ. kc) THEN
            dinv(i) = dinv(i) - l(indxr1)*dinv(kc)*u(indxc1)
            indxr1 = indxr1 + 1
            indxc1 = indxc1 + 1
            IF(indxr1 .LE. indxr2 .AND. indxc1 .LE. indxc2) GO TO 270
         ENDIF
C
 300  CONTINUE
C
C         Replace diagonal elements by their inverses.
CVD$ VECTOR
      DO 430 i=1,n
         dinv(i) = 1.0d0/dinv(i)
 430  CONTINUE
C
      RETURN
C------------- LAST LINE OF DSILUS FOLLOWS ----------------------------

dslugm()

subroutine dslugm	(	integer	N,
		double precision, dimension(n)	B,
		double precision, dimension(n)	X,
		integer	NELT,
		integer, dimension(nelt)	IA,
		integer, dimension(nelt)	JA,
		double precision, dimension(nelt)	A,
		integer	ISYM,
		integer	NSAVE,
		integer	ITOL,
		double precision	TOL,
		integer	ITMAX,
		integer	ITER,
		double precision	ERR,
		integer	IERR,
		integer	IUNIT,
		double precision, dimension(lenw)	RWORK,
		integer	LENW,
		integer, dimension(leniw)	IWORK,
		integer	LENIW
	)

C***BEGIN PROLOGUE  DSLUGM
C***PURPOSE  Incomplete LU GMRES iterative sparse Ax=b solver.
C            This routine uses the generalized minimum residual
C            (GMRES) method with incomplete LU factorization for
C            preconditioning to solve possibly non-symmetric linear
C            systems of the form: Ax = b.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D2A4, D2B4
C***TYPE      DOUBLE PRECISION (SSLUGM-S, DSLUGM-D)
C***KEYWORDS  GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
C             NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
C***AUTHOR  Brown, Peter, (LLNL), pnbrown@llnl.gov
C           Hindmarsh, Alan, (LLNL), alanh@llnl.gov
C           Seager, Mark K., (LLNL), seager@llnl.gov
C             Lawrence Livermore National Laboratory
C             PO Box 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C***DESCRIPTION
C
C *Usage:
C      INTEGER   N, NELT, IA(NELT), JA(NELT), ISYM, NSAVE, ITOL
C      INTEGER   ITMAX, ITER, IERR, IUNIT, LENW, IWORK(LENIW), LENIW
C      DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(LENW)
C
C      CALL DSLUGM(N, B, X, NELT, IA, JA, A, ISYM, NSAVE,
C     $     ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT,
C     $     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     :IN       Integer IA(NELT).
C JA     :IN       Integer JA(NELT).
C A      :IN       Double Precision 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 non-zero 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 NSAVE  :IN       Integer.
C         Number of direction vectors to save and orthogonalize against.
C         Must be greater than 1.
C ITOL   :IN       Integer.
C         Flag to indicate the type of convergence criterion used.
C         ITOL=0  Means the  iteration stops when the test described
C                 below on  the  residual RL  is satisfied.  This is
C                 the  "Natural Stopping Criteria" for this routine.
C                 Other values  of   ITOL  cause  extra,   otherwise
C                 unnecessary, computation per iteration and     are
C                 therefore  much less  efficient.  See  ISDGMR (the
C                 stop test routine) for more information.
C         ITOL=1  Means   the  iteration stops   when the first test
C                 described below on  the residual RL  is satisfied,
C                 and there  is either right  or  no preconditioning
C                 being used.
C         ITOL=2  Implies     that   the  user    is   using    left
C                 preconditioning, and the second stopping criterion
C                 below is used.
C         ITOL=3  Means the  iteration stops   when  the  third test
C                 described below on Minv*Residual is satisfied, and
C                 there is either left  or no  preconditioning begin
C                 used.
C         ITOL=11 is    often  useful  for   checking  and comparing
C                 different routines.  For this case, the  user must
C                 supply  the  "exact" solution or  a  very accurate
C                 approximation (one with  an  error much less  than
C                 TOL) through a common block,
C                     COMMON /DSLBLK/ SOLN( )
C                 If ITOL=11, iteration stops when the 2-norm of the
C                 difference between the iterative approximation and
C                 the user-supplied solution  divided by the  2-norm
C                 of the  user-supplied solution  is  less than TOL.
C                 Note that this requires  the  user to  set up  the
C                 "COMMON     /DSLBLK/ SOLN(LENGTH)"  in the calling
C                 routine.  The routine with this declaration should
C                 be loaded before the stop test so that the correct
C                 length is used by  the loader.  This procedure  is
C                 not standard Fortran and may not work correctly on
C                 your   system (although  it  has  worked  on every
C                 system the authors have tried).  If ITOL is not 11
C                 then this common block is indeed standard Fortran.
C TOL    :INOUT    Double Precision.
C         Convergence criterion, as described below.  If TOL is set
C         to zero on input, then a default value of 500*(the smallest
C         positive magnitude, machine epsilon) is used.
C ITMAX  :IN       Integer.
C         Maximum number of iterations.  This routine uses the default
C         of NRMAX = ITMAX/NSAVE to determine the when each restart
C         should occur.  See the description of NRMAX and MAXL in
C         DGMRES for a full and frightfully interesting discussion of
C         this topic.
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.  Letting norm() denote the Euclidean
C         norm, ERR is defined as follows...
C         If ITOL=0, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
C                               for right or no preconditioning, and
C                         ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
C                                norm(SB*(M-inverse)*B),
C                               for left preconditioning.
C         If ITOL=1, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
C                               since right or no preconditioning
C                               being used.
C         If ITOL=2, then ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
C                                norm(SB*(M-inverse)*B),
C                               since left preconditioning is being
C                               used.
C         If ITOL=3, then ERR =  Max  |(Minv*(B-A*X(L)))(i)/x(i)|
C                               i=1,n
C         If ITOL=11, then ERR = norm(SB*(X(L)-SOLN))/norm(SB*SOLN).
C IERR   :OUT      Integer.
C         Return error flag.
C               IERR = 0 => All went well.
C               IERR = 1 => Insufficient storage allocated for
C                           RGWK or IGWK.
C               IERR = 2 => Routine DPIGMR failed to reduce the norm
C                           of the current residual on its last call,
C                           and so the iteration has stalled.  In
C                           this case, X equals the last computed
C                           approximation.  The user must either
C                           increase MAXL, or choose a different
C                           initial guess.
C               IERR =-1 => Insufficient length for RGWK array.
C                           IGWK(6) contains the required minimum
C                           length of the RGWK array.
C               IERR =-2 => Inconsistent ITOL and JPRE values.
C         For IERR <= 2, RGWK(1) = RHOL, which is the norm on the
C         left-hand-side of the relevant stopping test defined
C         below associated with the residual for the current
C         approximation X(L).
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 of size LENW.
C LENW   :IN       Integer.
C         Length of the double precision workspace, RWORK.
C         LENW >= 1 + N*(NSAVE+7) +  NSAVE*(NSAVE+3)+NL+NU.
C         Here NL is the number of non-zeros in the lower triangle of
C         the matrix (including the diagonal) and NU is the number of
C         non-zeros in the upper triangle of the matrix (including the
C         diagonal).
C         For the recommended values,  RWORK  has size at least
C         131 + 17*N + NL + NU.
C IWORK  :INOUT    Integer IWORK(LENIW).
C         Used to hold pointers into the RWORK array.
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.
C         LENIW >= NL+NU+4*N+32.
C
C *Description:
C       DSLUGM solves a linear system A*X = B rewritten in the form:
C
C        (SB*A*(M-inverse)*(SX-inverse))*(SX*M*X) = SB*B,
C
C       with right preconditioning, or
C
C        (SB*(M-inverse)*A*(SX-inverse))*(SX*X) = SB*(M-inverse)*B,
C
C       with left preconditioning, where A is an n-by-n double precision
C       matrix, X and B are N-vectors, SB and SX are  diagonal scaling
C       matrices, and M is the Incomplete LU factorization of A.  It
C       uses preconditioned  Krylov subpace   methods  based on  the
C       generalized minimum residual  method (GMRES).   This routine
C       is a  driver  routine  which  assumes a SLAP   matrix   data
C       structure   and  sets  up  the  necessary  information to do
C       diagonal  preconditioning  and calls the main GMRES  routine
C       DGMRES for the   solution   of the linear   system.   DGMRES
C       optionally   performs  either  the full    orthogonalization
C       version of the  GMRES algorithm or  an incomplete variant of
C       it.  Both versions use restarting of the linear iteration by
C       default, although the user can disable this feature.
C
C       The GMRES  algorithm generates a sequence  of approximations
C       X(L) to the  true solution of the above  linear system.  The
C       convergence criteria for stopping the  iteration is based on
C       the size  of the  scaled norm of  the residual  R(L)  =  B -
C       A*X(L).  The actual stopping test is either:
C
C               norm(SB*(B-A*X(L))) .le. TOL*norm(SB*B),
C
C       for right preconditioning, or
C
C               norm(SB*(M-inverse)*(B-A*X(L))) .le.
C                       TOL*norm(SB*(M-inverse)*B),
C
C       for left preconditioning, where norm() denotes the Euclidean
C       norm, and TOL is  a positive scalar less  than one  input by
C       the user.  If TOL equals zero  when DSLUGM is called, then a
C       default  value  of 500*(the   smallest  positive  magnitude,
C       machine epsilon) is used.  If the  scaling arrays SB  and SX
C       are used, then  ideally they  should be chosen  so  that the
C       vectors SX*X(or SX*M*X) and  SB*B have all their  components
C       approximately equal  to  one in  magnitude.  If one wants to
C       use the same scaling in X  and B, then  SB and SX can be the
C       same array in the calling program.
C
C       The following is a list of the other routines and their
C       functions used by GMRES:
C       DGMRES  Contains the matrix structure independent driver
C               routine for GMRES.
C       DPIGMR  Contains the main iteration loop for GMRES.
C       DORTH   Orthogonalizes a new vector against older basis vectors.
C       DHEQR   Computes a QR decomposition of a Hessenberg matrix.
C       DHELS   Solves a Hessenberg least-squares system, using QR
C               factors.
C       RLCALC  Computes the scaled residual RL.
C       XLCALC  Computes the solution XL.
C       ISDGMR  User-replaceable stopping routine.
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 on
C       is being used and act accordingly.
C
C       =================== S L A P Triad format ===================
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
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 *Side Effects:
C       The SLAP Triad format (IA, JA, A) is modified internally to be
C       the SLAP Column format.  See above.
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 is attached to a file or terminal before calling
C     this routine with a non-zero value for IUNIT.  This routine does
C     not check for the validity of a non-zero IUNIT unit number.
C
C***REFERENCES  1. Peter N. Brown and A. C. Hindmarsh, Reduced Storage
C                  Matrix Methods in Stiff ODE Systems, Lawrence Liver-
C                  more National Laboratory Report UCRL-95088, Rev. 1,
C                  Livermore, California, June 1987.
C***ROUTINES CALLED  DCHKW, DGMRES, DS2Y, DSILUS, DSLUI, DSMV
C***REVISION HISTORY  (YYMMDD)
C   890404  DATE WRITTEN
C   890404  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   920407  COMMON BLOCK renamed DSLBLK.  (WRB)
C   920511  Added complete declaration section.  (WRB)
C   920929  Corrected format of references.  (FNF)
C   921019  Corrected NEL to NL.  (FNF)
C***END PROLOGUE  DSLUGM
C         The following is for optimized compilation on LLNL/LTSS Crays.
CLLL. OPTIMIZE
C     .. Parameters ..
      INTEGER locrb, locib
      parameter(locrb=1, locib=11)
C     .. Scalar Arguments ..
      DOUBLE PRECISION err, tol
      INTEGER ierr, isym, iter, itmax, itol, iunit, leniw, lenw, n,
     +        nelt, nsave
C     .. Array Arguments ..
      DOUBLE PRECISION a(nelt), b(n), rwork(lenw), x(n)
      INTEGER ia(nelt), iwork(leniw), ja(nelt)
C     .. Local Scalars ..
      INTEGER icol, j, jbgn, jend, locdin, locigw, locil, lociu, lociw,
     +        locjl, locju, locl, locnc, locnr, locrgw, locu, locw,
     +        myitol, nl, nu
C     .. External Subroutines ..
      EXTERNAL dchkw, dgmres, ds2y, dsilus, dslui, dsmv
C***FIRST EXECUTABLE STATEMENT  DSLUGM
C
      ierr = 0
      err  = 0
      IF( nsave.LE.1 ) THEN
         ierr = 3
         RETURN
      ENDIF
C
C         Change the SLAP input matrix IA, JA, A to SLAP-Column format.
      CALL ds2y( n, nelt, ia, ja, a, isym )
C
C         Count number of Non-Zero elements preconditioner ILU matrix.
C         Then set up the work arrays.  We assume MAXL=KMP=NSAVE.
      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
      locigw = locib
      locil = locigw + 20
      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
      locrgw = locu + nu
      locw = locrgw + 1+n*(nsave+6)+nsave*(nsave+3)
C
C         Check the workspace allocations.
      CALL dchkw( 'DSLUGM', 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         Perform the Incomplete LU Preconditioned Generalized Minimum
C         Residual iteration algorithm.  The following DGMRES
C         defaults are used MAXL = KMP = NSAVE, JSCAL = 0,
C         JPRE = -1, NRMAX = ITMAX/NSAVE
      iwork(locigw  ) = nsave
      iwork(locigw+1) = nsave
      iwork(locigw+2) = 0
      iwork(locigw+3) = -1
      iwork(locigw+4) = itmax/nsave
      myitol = 0
C
      CALL dgmres( n, b, x, nelt, ia, ja, a, isym, dsmv, dslui,
     $     myitol, tol, itmax, iter, err, ierr, iunit, rwork, rwork,
     $     rwork(locrgw), lenw-locrgw, iwork(locigw), 20,
     $     rwork, iwork )
C
      IF( iter.GT.itmax ) ierr = 2
      RETURN
C------------- LAST LINE OF DSLUGM FOLLOWS ----------------------------

dslui()

subroutine dslui	(	integer	N,
		double precision, dimension(n)	B,
		double precision, dimension(n)	X,
		integer	NELT,
		integer, dimension(nelt)	IA,
		integer, dimension(nelt)	JA,
		double precision, dimension(nelt)	A,
		integer	ISYM,
		double precision, dimension(*)	RWORK,
		integer, dimension(10)	IWORK
	)

C***BEGIN PROLOGUE  DSLUI
C***PURPOSE  SLAP MSOLVE for LDU Factorization.
C            This routine acts as an interface between the SLAP generic
C            MSOLVE calling convention and the routine that actually
C                           -1
C            computes  (LDU)  B = X.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D2E
C***TYPE      DOUBLE PRECISION (SSLUI-S, DSLUI-D)
C***KEYWORDS  ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM SOLVE,
C             SLAP, SPARSE
C***AUTHOR  Greenbaum, Anne, (Courant Institute)
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C       It is assumed that RWORK and IWORK have initialized with
C       the information required for DSLUI2:
C          IWORK(1) = Starting location of IL in IWORK.
C          IWORK(2) = Starting location of JL in IWORK.
C          IWORK(3) = Starting location of IU in IWORK.
C          IWORK(4) = Starting location of JU in IWORK.
C          IWORK(5) = Starting location of L in RWORK.
C          IWORK(6) = Starting location of DINV in RWORK.
C          IWORK(7) = Starting location of U in RWORK.
C       See the DESCRIPTION of DSLUI2 for details.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DSLUI2
C***REVISION HISTORY  (YYMMDD)
C   871119  DATE WRITTEN
C   881213  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   920511  Added complete declaration section.  (WRB)
C   921113  Corrected C***CATEGORY line.  (FNF)
C   930701  Updated CATEGORY section.  (FNF, WRB)
C***END PROLOGUE  DSLUI
C     .. Scalar Arguments ..
      INTEGER isym, n, nelt
C     .. Array Arguments ..
      DOUBLE PRECISION a(nelt), b(n), rwork(*), x(n)
      INTEGER ia(nelt), iwork(10), ja(nelt)
C     .. Local Scalars ..
      INTEGER locdin, locil, lociu, locjl, locju, locl, locu
C     .. External Subroutines ..
      EXTERNAL dslui2
C***FIRST EXECUTABLE STATEMENT  DSLUI
C
C         Pull out the locations of the arrays holding the ILU
C         factorization.
C
      locil = iwork(1)
      locjl = iwork(2)
      lociu = iwork(3)
      locju = iwork(4)
      locl = iwork(5)
      locdin = iwork(6)
      locu = iwork(7)
C
C         Solve the system LUx = b
      CALL dslui2(n, b, x, iwork(locil), iwork(locjl), rwork(locl),
     $     rwork(locdin), iwork(lociu), iwork(locju), rwork(locu) )
C
      RETURN
C------------- LAST LINE OF DSLUI FOLLOWS ----------------------------

dslui2()

subroutine dslui2	(	integer	N,
		double precision, dimension(n)	B,
		double precision, dimension(n)	X,
		integer, dimension(*)	IL,
		integer, dimension(*)	JL,
		double precision, dimension(*)	L,
		double precision, dimension(n)	DINV,
		integer, dimension(*)	IU,
		integer, dimension(*)	JU,
		double precision, dimension(*)	U
	)

      use omp_lib
C***BEGIN PROLOGUE  DSLUI2
C***PURPOSE  SLAP Backsolve for LDU Factorization.
C            Routine to solve a system of the form  L*D*U X = B,
C            where L is a unit lower triangular matrix, D is a diagonal
C            matrix, and U is a unit upper triangular matrix.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D2E
C***TYPE      DOUBLE PRECISION (SSLUI2-S, DSLUI2-D)
C***KEYWORDS  ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM SOLVE,
C             SLAP, SPARSE
C***AUTHOR  Greenbaum, Anne, (Courant Institute)
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     INTEGER N, IL(NL), JL(NL), IU(NU), JU(NU)
C     DOUBLE PRECISION B(N), X(N), L(NL), DINV(N), U(NU)
C
C     CALL DSLUI2( N, B, X, IL, JL, L, DINV, IU, JU, U )
C
C *Arguments:
C N      :IN       Integer
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right hand side.
C X      :OUT      Double Precision X(N).
C         Solution of L*D*U x = b.
C IL     :IN       Integer IL(NL).
C JL     :IN       Integer JL(NL).
C L      :IN       Double Precision L(NL).
C         IL, JL, L contain the unit  lower triangular factor of the
C         incomplete decomposition of some matrix stored in SLAP Row
C         format.  The diagonal of ones *IS* stored.  This structure
C         can   be   set  up  by   the  DSILUS  routine.   See   the
C         "Description", below  for more   details about   the  SLAP
C         format.  (NL is the number of non-zeros in the L array.)
C DINV   :IN       Double Precision DINV(N).
C         Inverse of the diagonal matrix D.
C IU     :IN       Integer IU(NU).
C JU     :IN       Integer JU(NU).
C U      :IN       Double Precision U(NU).
C         IU, JU, U contain the unit upper triangular factor  of the
C         incomplete decomposition  of  some  matrix stored in  SLAP
C         Column format.   The diagonal of ones  *IS* stored.   This
C         structure can be set up  by the DSILUS routine.  See   the
C         "Description", below   for  more   details about  the SLAP
C         format.  (NU is the number of non-zeros in the U array.)
C
C *Description:
C       This routine is supplied with  the SLAP package as a routine
C       to  perform  the  MSOLVE operation  in   the  SIR and   SBCG
C       iteration routines for  the  drivers DSILUR and DSLUBC.   It
C       must  be called  via   the  SLAP  MSOLVE  calling   sequence
C       convention interface routine DSLUI.
C         **** THIS ROUTINE ITSELF DOES NOT CONFORM TO THE ****
C               **** SLAP MSOLVE CALLING CONVENTION ****
C
C       IL, JL, L should contain the unit lower triangular factor of
C       the incomplete decomposition of the A matrix  stored in SLAP
C       Row format.  IU, JU, U should contain  the unit upper factor
C       of the  incomplete decomposition of  the A matrix  stored in
C       SLAP Column format This ILU factorization can be computed by
C       the DSILUS routine. The diagonals (which are all one's) are
C       stored.
C
C       =================== S L A P Column format ==================
C
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       ==================== S L A P Row format ====================
C
C       This routine requires  that the matrix A  be  stored  in the
C       SLAP  Row format.   In this format  the non-zeros are stored
C       counting across  rows (except for the diagonal  entry, which
C       must  appear first  in each  "row")  and  are stored  in the
C       double precision  array A.  In other words, for each row  in
C       the matrix  put the diagonal  entry in A.   Then put in  the
C       other  non-zero elements  going across  the row  (except the
C       diagonal) in order.  The JA array holds the column index for
C       each non-zero.  The IA array holds the offsets  into the JA,
C       A  arrays  for  the   beginning  of  each  row.    That  is,
C       JA(IA(IROW)),A(IA(IROW)) are the first elements of the IROW-
C       th row in  JA and A,  and  JA(IA(IROW+1)-1), A(IA(IROW+1)-1)
C       are  the last elements  of the  IROW-th row.   Note  that we
C       always have  IA(N+1) = NELT+1, where N is the number of rows
C       in the matrix  and  NELT is the  number of non-zeros  in the
C       matrix.
C
C       Here is an example of the SLAP Row storage format for a  5x5
C       Matrix (in the A and JA arrays '|' denotes the end of a row):
C
C           5x5 Matrix         SLAP Row 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 12 15 | 22 21 | 33 35 | 44 | 55 51 53
C       |21 22  0  0  0|  JA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  IA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C       With  the SLAP  format  the "inner  loops" of  this  routine
C       should vectorize   on machines with   hardware  support  for
C       vector gather/scatter operations.  Your compiler may require
C       a  compiler directive  to  convince   it that there  are  no
C       implicit vector  dependencies.  Compiler directives  for the
C       Alliant FX/Fortran and CRI CFT/CFT77 compilers  are supplied
C       with the standard SLAP distribution.
C
C***SEE ALSO  DSILUS
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   871119  DATE WRITTEN
C   881213  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   920511  Added complete declaration section.  (WRB)
C   921113  Corrected C***CATEGORY line.  (FNF)
C   930701  Updated CATEGORY section.  (FNF, WRB)
C***END PROLOGUE  DSLUI2
C     .. Scalar Arguments ..
      INTEGER n
C     .. Array Arguments ..
      DOUBLE PRECISION b(n), dinv(n), l(*), u(*), x(n)
      INTEGER il(*), iu(*), jl(*), ju(*)
C     .. Local Scalars ..
      INTEGER i, icol, irow, j, jbgn, jend
C***FIRST EXECUTABLE STATEMENT  DSLUI2
C
C         Solve  L*Y = B,  storing result in X, L stored by rows.
C
c$omp parallel default(none)
c$omp& shared(x,b,n)
c$omp& private(i)
c$omp do
      DO 10 i = 1, n
         x(i) = b(i)
 10   CONTINUE
c$omp end do
c$omp end parallel
c
      DO 30 irow = 2, n
         jbgn = il(irow)
         jend = il(irow+1)-1
         IF( jbgn.LE.jend ) THEN
CLLL. OPTION ASSERT (NOHAZARD)
CDIR$ IVDEP
CVD$ ASSOC
CVD$ NODEPCHK
            DO 20 j = jbgn, jend
               x(irow) = x(irow) - l(j)*x(jl(j))
 20         CONTINUE
         ENDIF
 30   CONTINUE
c
C
C         Solve  D*Z = Y,  storing result in X.
c$omp parallel default(none)
c$omp& shared(n,x,dinv)
c$omp& private(i)
c$omp do
      DO 40 i=1,n
         x(i) = x(i)*dinv(i)
 40   CONTINUE
c$omp end do
c$omp end parallel
C
C         Solve  U*X = Z, U stored by columns.
      DO 60 icol = n, 2, -1
         jbgn = ju(icol)
         jend = ju(icol+1)-1
         IF( jbgn.LE.jend ) THEN
CLLL. OPTION ASSERT (NOHAZARD)
CDIR$ IVDEP
CVD$ NODEPCHK
            DO 50 j = jbgn, jend
               x(iu(j)) = x(iu(j)) - u(j)*x(icol)
 50         CONTINUE
         ENDIF
 60   CONTINUE
C
      RETURN
C------------- LAST LINE OF DSLUI2 FOLLOWS ----------------------------

dsmv()

subroutine dsmv	(	integer	N,
		double precision, dimension(n)	X,
		double precision, dimension(n)	Y,
		integer	NELT,
		integer, dimension(nelt)	IA,
		integer, dimension(nelt)	JA,
		double precision, dimension(nelt)	A,
		integer	ISYM
	)

      use omp_lib
C***BEGIN PROLOGUE  DSMV
C***PURPOSE  SLAP Column Format Sparse Matrix Vector Product.
C            Routine to calculate the sparse matrix vector product:
C            Y = A*X.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D1B4
C***TYPE      DOUBLE PRECISION (SSMV-S, DSMV-D)
C***KEYWORDS  MATRIX VECTOR MULTIPLY, SLAP, SPARSE
C***AUTHOR  Greenbaum, Anne, (Courant Institute)
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     INTEGER  N, NELT, IA(NELT), JA(NELT), ISYM
C     DOUBLE PRECISION X(N), Y(N), A(NELT)
C
C     CALL DSMV(N, X, Y, NELT, IA, JA, A, ISYM )
C
C *Arguments:
C N      :IN       Integer.
C         Order of the Matrix.
C X      :IN       Double Precision X(N).
C         The vector that should be multiplied by the matrix.
C Y      :OUT      Double Precision Y(N).
C         The product of the matrix and the 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 should hold the matrix A in the SLAP Column
C         format.  See "Description", below.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all non-zero 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
C *Description
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       With  the SLAP  format  the "inner  loops" of  this  routine
C       should vectorize   on machines with   hardware  support  for
C       vector gather/scatter operations.  Your compiler may require
C       a  compiler directive  to  convince   it that there  are  no
C       implicit vector  dependencies.  Compiler directives  for the
C       Alliant FX/Fortran and CRI CFT/CFT77 compilers  are supplied
C       with the standard SLAP distribution.
C
C *Cautions:
C     This   routine   assumes  that  the matrix A is stored in SLAP
C     Column format.  It does not check  for  this (for  speed)  and
C     evil, ugly, ornery and nasty things  will happen if the matrix
C     data  structure  is,  in fact, not SLAP Column.  Beware of the
C     wrong data structure!!!
C
C***SEE ALSO  DSMTV
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   871119  DATE WRITTEN
C   881213  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   920511  Added complete declaration section.  (WRB)
C   930701  Updated CATEGORY section.  (FNF, WRB)
C***END PROLOGUE  DSMV
C     .. Scalar Arguments ..
      INTEGER isym, n, nelt
C     .. Array Arguments ..
      DOUBLE PRECISION a(nelt), x(n), y(n)
      INTEGER ia(nelt), ja(nelt)
C     .. Local Scalars ..
      INTEGER i, ibgn, icol, iend, irow, j, jbgn, jend
C***FIRST EXECUTABLE STATEMENT  DSMV
      if(isym.ne.1) then
C
C         Zero out the result vector.
C
c$omp parallel default(none)
c$omp& shared(n,y,ja,ia,a,x)
c$omp& private(i,icol,ibgn,iend)
c$omp do
         DO 10 i = 1, n
            y(i) = 0
 10      CONTINUE
c$omp end do
C
C         Multiply by A.
C
c$omp do
         DO 30 icol = 1, n
            ibgn = ja(icol)
            iend = ja(icol+1)-1
            DO 20 i = ibgn, iend
c$omp atomic
               y(ia(i)) = y(ia(i)) + a(i)*x(icol)
 20         CONTINUE
 30      CONTINUE
c$omp end do
c$omp end parallel
C
      else
C
C         The matrix is symmetric.  Need to get the other half in...
C         This loops assumes that the diagonal is the first entry in
C         each column.
C
C
C         Zero out the result vector.
C
c$omp parallel default(none)
c$omp& shared(n,y,ja,ia,a,x)
c$omp& private(i,j,icol,irow,ibgn,iend,jbgn,jend)
c$omp do
         DO 100 i = 1, n
            y(i) = 0
 100     CONTINUE
c$omp end do
C
C         Multiply by A.
C
c$omp do
         DO 300 icol = 1, n
            ibgn = ja(icol)
            iend = ja(icol+1)-1
            DO 200 i = ibgn, iend
c$omp atomic
               y(ia(i)) = y(ia(i)) + a(i)*x(icol)
 200        CONTINUE
 300     CONTINUE
c$omp end do
C
c$omp do
         DO 50 irow = 1, n
            jbgn = ja(irow)+1
            jend = ja(irow+1)-1
            IF( jbgn.GT.jend ) GOTO 50
            DO 40 j = jbgn, jend
               y(irow) = y(irow) + a(j)*x(ia(j))
 40         CONTINUE
 50      CONTINUE
c$omp end do
c$omp end parallel
      ENDIF
      RETURN
C------------- LAST LINE OF DSMV FOLLOWS ----------------------------

qs2i1d()

subroutine qs2i1d	(	integer, dimension(n)	IA,
		integer, dimension(n)	JA,
		double precision, dimension(n)	A,
		integer	N,
		integer	KFLAG
	)

C***BEGIN PROLOGUE  QS2I1D
C***SUBSIDIARY
C***PURPOSE  Sort an integer array, moving an integer and DP array.
C            This routine sorts the integer array IA and makes the same
C            interchanges in the integer array JA and the double pre-
C            cision array A.  The array IA may be sorted in increasing
C            order or decreasing order.  A slightly modified QUICKSORT
C            algorithm is used.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  N6A2A
C***TYPE      DOUBLE PRECISION (QS2I1R-S, QS2I1D-D)
C***KEYWORDS  SINGLETON QUICKSORT, SLAP, SORT, SORTING
C***AUTHOR  Jones, R. E., (SNLA)
C           Kahaner, D. K., (NBS)
C           Seager, M. K., (LLNL) seager@llnl.gov
C           Wisniewski, J. A., (SNLA)
C***DESCRIPTION
C     Written by Rondall E Jones
C     Modified by John A. Wisniewski to use the Singleton QUICKSORT
C     algorithm. date 18 November 1976.
C
C     Further modified by David K. Kahaner
C     National Bureau of Standards
C     August, 1981
C
C     Even further modification made to bring the code up to the
C     Fortran 77 level and make it more readable and to carry
C     along one integer array and one double precision array during
C     the sort by
C     Mark K. Seager
C     Lawrence Livermore National Laboratory
C     November, 1987
C     This routine was adapted from the ISORT routine.
C
C     ABSTRACT
C         This routine sorts an integer array IA and makes the same
C         interchanges in the integer array JA and the double precision
C         array A.
C         The array IA may be sorted in increasing order or decreasing
C         order.  A slightly modified quicksort algorithm is used.
C
C     DESCRIPTION OF PARAMETERS
C        IA - Integer array of values to be sorted.
C        JA - Integer array to be carried along.
C         A - Double Precision array to be carried along.
C         N - Number of values in integer array IA to be sorted.
C     KFLAG - Control parameter
C           = 1 means sort IA in INCREASING order.
C           =-1 means sort IA in DECREASING order.
C
C***SEE ALSO  DS2Y
C***REFERENCES  R. C. Singleton, Algorithm 347, An Efficient Algorithm
C                 for Sorting With Minimal Storage, Communications ACM
C                 12:3 (1969), pp.185-7.
C***ROUTINES CALLED  XERMSG
C***REVISION HISTORY  (YYMMDD)
C   761118  DATE WRITTEN
C   890125  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   900805  Changed XERROR calls to calls to XERMSG.  (RWC)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   910506  Made subsidiary to DS2Y and corrected reference.  (FNF)
C   920511  Added complete declaration section.  (WRB)
C   920929  Corrected format of reference.  (FNF)
C   921012  Corrected all f.p. constants to double precision.  (FNF)
C***END PROLOGUE  QS2I1D
CVD$R NOVECTOR
CVD$R NOCONCUR
C     .. Scalar Arguments ..
      INTEGER kflag, n
C     .. Array Arguments ..
      DOUBLE PRECISION a(n)
      INTEGER ia(n), ja(n)
C     .. Local Scalars ..
      DOUBLE PRECISION r, ta, tta
      INTEGER i, iit, ij, it, j, jjt, jt, k, kk, l, m, nn
C     .. Local Arrays ..
      INTEGER il(21), iu(21)
C     .. External Subroutines ..
      EXTERNAL xermsg
C     .. Intrinsic Functions ..
      INTRINSIC abs, int
C***FIRST EXECUTABLE STATEMENT  QS2I1D
      nn = n
      IF (nn.LT.1) THEN
         CALL xermsg ('SLATEC', 'QS2I1D',
     $      'The number of values to be sorted was not positive.', 1, 1)
         RETURN
      ENDIF
      IF( n.EQ.1 ) RETURN
      kk = abs(kflag)
      IF ( kk.NE.1 ) THEN
         CALL xermsg ('SLATEC', 'QS2I1D',
     $      'The sort control parameter, K, was not 1 or -1.', 2, 1)
         RETURN
      ENDIF
C
C     Alter array IA to get decreasing order if needed.
C
      IF( kflag.LT.1 ) THEN
         DO 20 i=1,nn
            ia(i) = -ia(i)
 20      CONTINUE
      ENDIF
C
C     Sort IA and carry JA and A along.
C     And now...Just a little black magic...
      m = 1
      i = 1
      j = nn
      r = .375d0
 210  IF( r.LE.0.5898437d0 ) THEN
         r = r + 3.90625d-2
      ELSE
         r = r-.21875d0
      ENDIF
 225  k = i
C
C     Select a central element of the array and save it in location
C     it, jt, at.
C
      ij = i + int((j-i)*r)
      it = ia(ij)
      jt = ja(ij)
      ta = a(ij)
C
C     If first element of array is greater than it, interchange with it.
C
      IF( ia(i).GT.it ) THEN
         ia(ij) = ia(i)
         ia(i)  = it
         it     = ia(ij)
         ja(ij) = ja(i)
         ja(i)  = jt
         jt     = ja(ij)
         a(ij)  = a(i)
         a(i)   = ta
         ta     = a(ij)
      ENDIF
      l=j
C
C     If last element of array is less than it, swap with it.
C
      IF( ia(j).LT.it ) THEN
         ia(ij) = ia(j)
         ia(j)  = it
         it     = ia(ij)
         ja(ij) = ja(j)
         ja(j)  = jt
         jt     = ja(ij)
         a(ij)  = a(j)
         a(j)   = ta
         ta     = a(ij)
C
C     If first element of array is greater than it, swap with it.
C
         IF ( ia(i).GT.it ) THEN
            ia(ij) = ia(i)
            ia(i)  = it
            it     = ia(ij)
            ja(ij) = ja(i)
            ja(i)  = jt
            jt     = ja(ij)
            a(ij)  = a(i)
            a(i)   = ta
            ta     = a(ij)
         ENDIF
      ENDIF
C
C     Find an element in the second half of the array which is
C     smaller than it.
C
  240 l=l-1
      IF( ia(l).GT.it ) GO TO 240
C
C     Find an element in the first half of the array which is
C     greater than it.
C
  245 k=k+1
      IF( ia(k).LT.it ) GO TO 245
C
C     Interchange these elements.
C
      IF( k.LE.l ) THEN
         iit   = ia(l)
         ia(l) = ia(k)
         ia(k) = iit
         jjt   = ja(l)
         ja(l) = ja(k)
         ja(k) = jjt
         tta   = a(l)
         a(l)  = a(k)
         a(k)  = tta
         GOTO 240
      ENDIF
C
C     Save upper and lower subscripts of the array yet to be sorted.
C
      IF( l-i.GT.j-k ) THEN
         il(m) = i
         iu(m) = l
         i = k
         m = m+1
      ELSE
         il(m) = k
         iu(m) = j
         j = l
         m = m+1
      ENDIF
      GO TO 260
C
C     Begin again on another portion of the unsorted array.
C
  255 m = m-1
      IF( m.EQ.0 ) GO TO 300
      i = il(m)
      j = iu(m)
  260 IF( j-i.GE.1 ) GO TO 225
      IF( i.EQ.j ) GO TO 255
      IF( i.EQ.1 ) GO TO 210
      i = i-1
  265 i = i+1
      IF( i.EQ.j ) GO TO 255
      it = ia(i+1)
      jt = ja(i+1)
      ta =  a(i+1)
      IF( ia(i).LE.it ) GO TO 265
      k=i
  270 ia(k+1) = ia(k)
      ja(k+1) = ja(k)
      a(k+1)  =  a(k)
      k = k-1
      IF( it.LT.ia(k) ) GO TO 270
      ia(k+1) = it
      ja(k+1) = jt
      a(k+1)  = ta
      GO TO 265
C
C     Clean up, if necessary.
C
  300 IF( kflag.LT.1 ) THEN
         DO 310 i=1,nn
            ia(i) = -ia(i)
 310     CONTINUE
      ENDIF
      RETURN
C------------- LAST LINE OF QS2I1D FOLLOWS ----------------------------