c----------------------------------------------------------------------c
c                          S P A R S K I T                             c
c----------------------------------------------------------------------c
c                   ITERATIVE SOLVERS MODULE                           c
c----------------------------------------------------------------------c
c contents: (May 7, 1990)  -- revised ilut Feb, 11, 1992               c
c----------                                                            c
c                                                                      c
c pgmres  : preconditioned GMRES solver                                c
c ilut    : incomplete LU factorization with dual truncation strategy  c
c ilutp   : ILUT with partial pivoting                                 c
c ilu0    : simple ILU(0) preconditioning                              c
c milu0   : MILU(0) preconditioning                                    c
c lusol   : forward followed by backward triangular solve (Precond.)   c
c lutsol  : solving v = (LU)^{-T} u (used for preconditioning)         c
c qsplit  : quick split routine used by ilut to sort out the k largest c
c           elements in absolute value                                 c
c                                                                      c
c Note: all preconditioners are preprocessors to pgmres.               c
c usage: call preconditioner then call pgmres.                         c
c                                                                      c
c----------------------------------------------------------------------c
       subroutine pgmres (n, im, rhs, sol, vv, eps, maxits, iout,
     *                    aa, ja, ia, alu, jlu, ju, ierr)
c-----------------------------------------------------------------------
       implicit real*8 (a-h,o-z)
       integer n, im, maxits, iout, ierr, ja(*), ia(n+1), jlu(*), ju(n)
       real*8 vv(n,*), rhs(n), sol(n), aa(*), alu(*), eps
c----------------------------------------------------------------------*
c                                                                      *
c                 *** ILUT - Preconditioned GMRES ***                  *
c                                                                      *
c----------------------------------------------------------------------*
c This is a simple version of the ILUT preconditioned GMRES algorithm. *
c The ILUT preconditioner uses a dual strategy for dropping elements   *
c instead  of the usual level of-fill-in approach. See details in ILUT *
c subroutine documentation. PGMRES uses the L and U matrices generated *
c from the subroutine ILUT to precondition the GMRES algorithm.        *
c The preconditioning is applied to the right. The stopping criterion  *
c utilized is based simply on reducing the residual norm by epsilon.   *
c This preconditioning is more reliable than ilu0 but requires more    *
c storage. It seems to be much less prone to difficulties related to   *
c strong nonsymmetries in the matrix. We recommend using a nonzero tol *
c (tol=.005 or .001 usually give good results) in ILUT. Use a large    *
c lfil whenever possible (e.g. lfil = 5 to 10). The higher lfil the    *
c more reliable the code is. Efficiency may also be much improved.     *
c Note that lfil=n and tol=0.0 in ILUT  will yield the same factors as *
c Gaussian elimination without pivoting.                               *
c                                                                      *
c ILU(0) and MILU(0) are also provided for comparison purposes         *
c USAGE: first call ILUT or ILU0 or MILU0 to set up preconditioner and *
c then call pgmres.                                                    *
c----------------------------------------------------------------------*
c Coded by Y. Saad - This version dated May, 7, 1990.                  *
c----------------------------------------------------------------------*
c parameters                                                           *
c-----------                                                           *
c on entry:                                                            *
c==========                                                            *
c                                                                      *
c n     == integer. The dimension of the matrix.                       *
c im    == size of krylov subspace:  should not exceed 50 in this      *
c          version (can be reset by changing parameter command for     *
c          kmax below)                                                 *
c rhs   == real vector of length n containing the right hand side.     *
c          Destroyed on return.                                        *
c sol   == real vector of length n containing an initial guess to the  *
c          solution on input. approximate solution on output           *
c eps   == tolerance for stopping criterion. process is stopped        *
c          as soon as ( ||.|| is the euclidean norm):                  *
c          || current residual||/||initial residual|| <= eps           *
c maxits== maximum number of iterations allowed                        *
c iout  == output unit number number for printing intermediate results *
c          if (iout .le. 0) nothing is printed out.                    *
c                                                                      *
c aa, ja,                                                              *
c ia    == the input matrix in compressed sparse row format:           *
c          aa(1:nnz)  = nonzero elements of A stored row-wise in order *
c          ja(1:nnz) = corresponding column indices.                   *
c          ia(1:n+1) = pointer to beginning of each row in aa and ja.  *
c          here nnz = number of nonzero elements in A = ia(n+1)-ia(1)  *
c                                                                      *
c alu,jlu== A matrix stored in Modified Sparse Row format containing   *
c           the L and U factors, as computed by subroutine ilut.       *
c                                                                      *
c ju     == integer array of length n containing the pointers to       *
c           the beginning of each row of U in alu, jlu as computed     *
c           by subroutine ILUT.                                        *
c                                                                      *
c on return:                                                           *
c==========                                                            *
c sol   == contains an approximate solution (upon successful return).  *
c ierr  == integer. Error message with the following meaning.          *
c          ierr = 0 --> successful return.                             *
c          ierr = 1 --> convergence not achieved in itmax iterations.  *
c          ierr =-1 --> the initial guess seems to be the exact        *
c                       solution (initial residual computed was zero)  *
c                                                                      *
c----------------------------------------------------------------------*
c                                                                      *
c work arrays:                                                         *
c=============                                                         *
c vv    == work array of length  n x (im+1) (used to store the Arnoli  *
c          basis)                                                      *
c----------------------------------------------------------------------*
c subroutines called :                                                 *
c amux   : SPARSKIT routine to do the matrix by vector multiplication  *
c          delivers y=Ax, given x  -- see SPARSKIT/BLASSM/amux         *
c lusol : combined forward and backward solves (Preconditioning ope.) *
c BLAS1  routines.                                                     *
c----------------------------------------------------------------------*
       parameter (kmax=50)
       real*8 hh(kmax+1,kmax), c(kmax), s(kmax), rs(kmax+1),t
c-------------------------------------------------------------
c arnoldi size should not exceed kmax=50 in this version..
c to reset modify paramter kmax accordingly.
c-------------------------------------------------------------
       data epsmac/1.d-16/
       n1 = n + 1
       its = 0
c-------------------------------------------------------------
c outer loop starts here..
c-------------- compute initial residual vector --------------
       call amux (n, sol, vv, aa, ja, ia)
       do 21 j=1,n
          vv(j,1) = rhs(j) - vv(j,1)
 21    continue
c-------------------------------------------------------------
 20    ro = dnrm2(n, vv, 1)
c      if (iout .gt. 0 .and. its .eq. 0)
c    *      write(iout, 199) its, ro
       if (ro .eq. 0.0d0) goto 999
       t = 1.0d0/ ro
       do 210 j=1, n
          vv(j,1) = vv(j,1)*t
 210   continue
       if (its .eq. 0) eps1=eps*ro
c     ** initialize 1-st term  of rhs of hessenberg system..
       rs(1) = ro
       i = 0
 4     i=i+1
       its = its + 1
       i1 = i + 1
       call lusol (n, vv(1,i), rhs, alu, jlu, ju)
       call amux (n, rhs, vv(1,i1), aa, ja, ia)
c-----------------------------------------
c     modified gram - schmidt...
c-----------------------------------------
       do 55 j=1, i
          t = ddot(n, vv(1,j),1,vv(1,i1),1)
          hh(j,i) = t
          call daxpy(n, -t, vv(1,j), 1, vv(1,i1), 1)
 55    continue
       t = dnrm2(n, vv(1,i1), 1)
       hh(i1,i) = t
       if ( t .eq. 0.0d0) goto 58
       t = 1.0d0/t
       do 57  k=1,n
          vv(k,i1) = vv(k,i1)*t
 57    continue
c
c     done with modified gram schimd and arnoldi step..
c     now  update factorization of hh
c
 58    if (i .eq. 1) goto 121
c--------perfrom previous transformations  on i-th column of h
       do 66 k=2,i
          k1 = k-1
          t = hh(k1,i)
          hh(k1,i) = c(k1)*t + s(k1)*hh(k,i)
          hh(k,i) = -s(k1)*t + c(k1)*hh(k,i)
 66    continue
 121   gam = sqrt(hh(i,i)**2 + hh(i1,i)**2)
c
c     if gamma is zero then any small value will do...
c     will affect only residual estimate
c
       if (gam .eq. 0.0d0) gam = epsmac
c
c     get  next plane rotation
c
       c(i) = hh(i,i)/gam
       s(i) = hh(i1,i)/gam
       rs(i1) = -s(i)*rs(i)
       rs(i) =  c(i)*rs(i)
c
c     detrermine residual norm and test for convergence-
c
       hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i)
       ro = abs(rs(i1))
 131   format(1h ,2e14.4)
c      if (iout .gt. 0)
c    *      write(iout, 199) its, ro
       if (i .lt. im .and. (ro .gt. eps1))  goto 4
c
c     now compute solution. first solve upper triangular system.
c
       rs(i) = rs(i)/hh(i,i)
       do 30 ii=2,i
          k=i-ii+1
          k1 = k+1
          t=rs(k)
          do 40 j=k1,i
             t = t-hh(k,j)*rs(j)
 40       continue
          rs(k) = t/hh(k,k)
 30    continue
c
c     form linear combination of v(*,i)'s to get solution
c
       t = rs(1)
       do 15 k=1, n
          rhs(k) = vv(k,1)*t
 15    continue
       do 16 j=2, i
          t = rs(j)
          do 161 k=1, n
             rhs(k) = rhs(k)+t*vv(k,j)
 161      continue
 16    continue
c
c     call preconditioner.
c
       call lusol (n, rhs, rhs, alu, jlu, ju)
       do 17 k=1, n
          sol(k) = sol(k) + rhs(k)
 17    continue
c
c     restart outer loop  when necessary
c
       if (ro .le. eps1) goto 990
       if (its .ge. maxits) goto 991
c
c     else compute residual vector and continue..
c
       do 24 j=1,i
          jj = i1-j+1
          rs(jj-1) = -s(jj-1)*rs(jj)
          rs(jj) = c(jj-1)*rs(jj)
 24    continue
       do 25  j=1,i1
          t = rs(j)
          if (j .eq. 1)  t = t-1.0d0
          call daxpy (n, t, vv(1,j), 1,  vv, 1)
 25    continue
 199   format('   its =', i4, ' res. norm =', d20.6)
c     restart outer loop.
       goto 20
 990   ierr = 0
       return
 991   ierr = 1
       return
 999   continue
       ierr = -1
       return
c-----------------end of pgmres ---------------------------------------
c-----------------------------------------------------------------------
       end
c-----------------------------------------------------------------------
       subroutine ilut (n,a,ja,ia,lfil,tol,alu,jlu,ju,iwk,
     *                  wu,wl,jr,jwl,jwu,ierr)
c-----------------------------------------------------------------------
       implicit real*8 (a-h,o-z)
       real*8 a(*), alu(*), wu(n+1), wl(n), tol
       integer ja(*),ia(n+1),jlu(*),ju(n),jr(n), jwu(n),
     *      jwl(n), n, lfil, iwk, ierr
c----------------------------------------------------------------------*
c                      *** ILUT preconditioner ***                     *
c                      ---------------------------                     *
c      incomplete LU factorization with dual truncation mechanism      *
c      VERSION 2 : sorting  done for both L and U.                     *
c                                                                      *
c----------------------------------------------------------------------*
c---- coded by Youcef Saad May, 5, 1990. ------------------------------*
c---- Dual drop-off strategy works as follows.                         *
c                                                                      *
c     1) Theresholding in L and U as set by tol. Any element whose size*
c        is less than some tolerance (relative to the norm of current  *
c        row in u) is dropped.                                         *
c                                                                      *
c     2) Keeping only the largest lfil+il(i) elements in the i-th row  *
c        of L and the largest lfil+iu(i) elements in the i-th row of   *
c        U where il(i), iu(i) are the original number of nonzero       *
c        elements of the L-part and the U-part of the i-th row of A    *
c                                                                      *
c Flexibility: one can use tol=0 to get a strategy based on keeping the*
c largest elements in each row of L and U. Taking tol .ne. 0 but lfil=n*
c will give the usual threshold strategy (however, fill-in is then     *
c impredictible).                                                      *
c                                                                      *
c----------------------------------------------------------------------*
c PARAMETERS
c-----------
c
c on entry:
c==========
c n       = integer. The dimension of the matrix A.
c
c a,ja,ia = matrix stored in Compressed Sparse Row format.
c
c lfil    = integer. The fill-in parameter. Each row of L and
c           each row of U will have a maximum of lfil elements
c           in addition to their original number of nonzero elements.
c           Thus storage can be determined beforehand.
c           lfil must be .ge. 0.
c
c iwk     = integer. The minimum length of arrays alu and jlu
c
c On return:
c===========
c
c alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
c           the L and U factors together. The diagonal (stored in
c           alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
c           contains the i-th row of L (excluding the diagonal entry=1)
c           followed by the i-th row of U.
c
c ju      = integer array of length n containing the pointers to
c           the beginning of each row of U in the matrix alu,jlu.
c
c ierr    = integer. Error message with the following meaning.
c           ierr  = 0    --> successful return.
c           ierr .gt. 0  --> zero pivot encountered at step number ierr.
c           ierr  = -1   --> Error. input matrix may be wrong.
c                            (The elimination process has generated a
c                            row in L or U whose length is .gt.  n.)
c           ierr  = -2   --> The matrix L overflows the array al.
c           ierr  = -3   --> The matrix U overflows the array alu.
c           ierr  = -4   --> Illegal value for lfil.
c           ierr  = -5   --> zero row encountered.
c
c work arrays:
c=============
c jr,jwu,jwl 	  = integer work arrays of length n.
c wu, wl          = real work arrays of length n+1, and n resp.
c
c Notes:
c ------
c A must have all nonzero diagonal elements.
c-----------------------------------------------------------------------
        if (lfil .lt. 0) goto 998
c-------------------------------
c initialize ju0 (points to next element to be added to alu,jlu)
c and pointer.
c
	ju0 = n+2
	jlu(1) = ju0
c
c  integer double pointer array.
c
	do 1 j=1, n
		jr(j)  = 0
 1           continue
c-----------------------------------------------------------------------
c  beginning of main loop.
c-----------------------------------------------------------------------
	do 500 ii = 1, n
           j1 = ia(ii)
           j2 = ia(ii+1) - 1
           tnorm = 0.0d0
           do 501 k=j1,j2
              tnorm = tnorm+abs(a(k))
 501          continue
              if (tnorm .eq. 0.0) goto 999
              tnorm = tnorm/real(j2-j1+1)
c
c--- unpack L-part and U-part of row of A in arrays wl, wu --
c
       lenu = 1
       lenl = 0
       jwu(1) = ii
       wu(1) = 0.0
       jr(ii) = 1
c
       do 170  j = j1, j2
           k = ja(j)
           t = a(j)
           if (abs(t) .lt. tol*tnorm .and. k .ne. ii) goto 170
           if (k .lt. ii) then
              lenl = lenl+1
              jwl(lenl) = k
              wl(lenl) = t
              jr(k) = lenl
           else if (k .eq. ii) then
              wu(1) = t
           else
              lenu = lenu+1
              jwu(lenu) = k
              wu(lenu) = t
              jr(k) = lenu
           endif
 170	continue
        tnorm = tnorm/real(j2-j1+1)
        lenl0 = lenl
        lenu0 = lenu
        jj = 0
        nl = 0
c-------------------------------------------------------------------
c---------------------- eliminate previous rows --------------------
c-------------------------------------------------------------------
 150    jj = jj+1
        if (jj .gt. lenl) goto 160
c-------------------------------------------------------------------
c in order to do the elimination in the correct order we need to
c exchange the current row number with the one that has
c smallest column number, among jj,jj+1,...,lenl.
c-------------------------------------------------------------------
        jrow = jwl(jj)
        k = jj
c
c determine smallest column index
c
        do 151 j=jj+1,lenl
           if (jwl(j) .lt. jrow) then
              jrow = jwl(j)
              k = j
           endif
 151    continue
c
c exchange in jwl
c
       if (k .ne. jj) then
           j = jwl(jj)
           jwl(jj) = jwl(k)
           jwl(k) = j
c
c exchange in jr
c
           jr(jrow) = jj
           jr(j) = k
c
c exchange in wl
c
           s = wl(jj)
           wl(jj) = wl(k)
           wl(k) = s
        endif
c
        if (jrow .ge. ii) goto 160
c---------get the multiplier for row to be eliminated: jrow
        fact = wl(jj)*alu(jrow)
c zero out element in row by setting jr(jrow) = 0
        jr(jrow) = 0
        if (abs(fact)*wu(n+2-jrow) .le. tol*tnorm) goto 150
c-------------------------------------------------------------------
c------------ combine current row and row jrow ---------------------
c-------------------------------------------------------------------
        do 203 k = ju(jrow), jlu(jrow+1)-1
           s = fact*alu(k)
           j = jlu(k)
           jpos = jr(j)
c
c if fill-in element is small then disregard:
c
           if (abs(s) .lt. tol*tnorm .and. jpos .eq. 0) goto 203
           if (j .ge. ii) then
c
c     dealing with upper part.
c
              if (jpos .eq. 0) then
c     this is a fill-in element
                 lenu = lenu+1
                 if (lenu .gt. n) goto 995
                 jwu(lenu) = j
                 jr(j) = lenu
                 wu(lenu) = - s
              else
c     no fill-in element --
                 wu(jpos) = wu(jpos) - s
              endif
           else
c
c     dealing with lower part.
c
              if (jpos .eq. 0) then
c     this is a fill-in element
                 lenl = lenl+1
                 if (lenl .gt. n) goto 995
                 jwl(lenl) = j
                 jr(j) = lenl
                 wl(lenl) = - s
              else
c     no fill-in element --
                 wl(jpos) = wl(jpos) - s
              endif
           endif
 203	continue
        nl = nl+1
        wl(nl) = fact
        jwl(nl)  = jrow
	goto 150
c----------------------------------------------------------
c------------ update l-matrix -----------------------------
c----------------------------------------------------------
 160    len = min0(nl,lenl0+lfil)
c 160    len = min0(nl,lfil)

  	call qsplit (wl,jwl,nl,len)
c
        do 204 k=1, len
           if (ju0 .gt. iwk) goto 996
           alu(ju0) =  wl(k)
           jlu(ju0) =  jwl(k)
           ju0 = ju0+1
 204    continue
c
c  save pointer to beginning of row ii of U
c
        ju(ii) = ju0
c
c  reset double-pointer jr to zero (L-part - except first
c  jj-1 elements which have already been reset)
	do 306 k= jj, lenl
              jr(jwl(k)) = 0
 306	continue
        len = min0(lenu,lenu0+lfil+1) ! +1 added 11/23/94
c        len = min0(lenu,lfil)
	call qsplit (wu(2), jwu(2), lenu-1,len)
c----------------------------------------------------------
c------------ update u-matrix -----------------------------
c----------------------------------------------------------
        t = abs(wu(1))
        if (len + ju0 .gt. iwk) goto 997
        do 302 k=2, len
           jlu(ju0) = jwu(k)
           alu(ju0) = wu(k)
           t = t + abs(wu(k) )
           ju0 = ju0+1
 302	continue
c
c     save norm in wu (backwards). Norm is in fact average abs value
c
        wu(n+2-ii) = t / real(len+1)
c
c     store inverse of diagonal element of u
c
        if (wu(1) .eq. 0.0) wu(1) = (0.0001 + tol)*tnorm
c
        alu(ii) = 1.0d0/ wu(1)
c
c     update pointer to beginning of next row of U.
c
	jlu(ii+1) = ju0
c
c     reset double-pointer jr to zero (U-part)
c
	do 308 k=1, lenu
           jr(jwu(k)) = 0
 308	continue
c-----------------------------------------------------------------------
c     end main loop
c-----------------------------------------------------------------------
 500  continue
        ierr = 0
        return
c
c     zero pivot :
c
c 900    ierr = ii
c        return
c
c     incomprehensible error. Matrix must be wrong.
c
 995    ierr = -1
        return
c
c     insufficient storage in L.
c
 996    ierr = -2
        return
c
c     insufficient storage in U.
c
 997    ierr = -3
        return
c
c     illegal lfil entered.
c
 998    ierr = -4
        return
c
c     zero row encountered
c
 999    ierr = -5
        return
c---------------- end of ilut  -----------------------------------------
c-----------------------------------------------------------------------
        end
c----------------------------------------------------------------------
      subroutine ilutp (n,a,ja,ia,lfil,droptol,permtol,mbloc,alu,jlu,
     *    ju,iwk,wu,wl,jr,jwl,jwu,iperm,ierr)
c-----------------------------------------------------------------------
c     SPARSKIT ROUTINE   --  ILUT with PIVOTING --
c-----------------------------------------------------------------------
c      implicit none
      integer n,ja(*),ia(n+1),lfil,jlu(*),ju(n),jr(n),jwu(n),
     *     jwl(n),iwk,iperm(2*n),ierr
      real*8 a(*), alu(*), wu(n+1), wl(n), droptol
c----------------------------------------------------------------------*
c                      *** ILUT preconditioner ***                     *
c                      ---------------------------                     *
c      incomplete LU factorization with dual truncation mechanism      *
c      VERSION 2 : sorting  done for both L and U.                     *
c                                                                      *
c----------------------------------------------------------------------*
c---- coded by Youcef Saad Sep 8, 1993 . ------------------------------*
c---- Dual drop-off strategy works as follows.                         *
c                                                                      *
c 1) Theresholding in L and U as set by droptol. Any element whose size*
c    is less than some tolerance (relative to the norm of current      *
c    row in u) is dropped.                                             *
c                                                                      *
c 2) Keeping only the largest lfil+il(i) elements in the i-th row      *
c    of L and the largest lfil+iu(i) elements in the i-th row of       *
c    U where il(i), iu(i) are the original number of nonzero           *
c    elements of the L-part and the U-part of the i-th row of A        *
c                                                                      *
c column permuting is used --                                          *
c  see also comments in ilut                                           *
c----------------------------------------------------------------------*
c PARAMETERS
c-----------
c
c on entry:
c==========
c n       = integer. The dimension of the matrix A.
c
c a,ja,ia = matrix stored in Compressed Sparse Row format.
c           ONE RETURN THE COLUMNS OF A ARE PERMUTED.
c
c lfil    = integer. The fill-in parameter. Each row of L and
c           each row of U will have a maximum of lfil elements
c           in addition to their original number of nonzero elements.
c           Thus storage can be determined beforehand.
c           lfil must be .ge. 0.
c
c droptol = tolerance used for dropping elements in L and U.
c           elements are dropped if they are .lt. norm(row) x droptol
c           row = row being eliminated
c
c permtol = tolerance ratio used for determning whether to permute
c           two columns. We will permute two columns only if
c           a(i,j)*permtol .gt. a(i,i) [good values 0.1 to 0.01]
c
c mbloc   = if desired, permuting can be done only within the diagonal
c           blocks of size mbloc. Useful for PDE problems with several
c           degrees of freedom.. If feature not wanted take mbloc=n.
c
c iwk     = integer. The minimum length of arrays alu and jlu
c           to work properly, the code requires that iwk be
c
c                      .ge. nnz + 2*lfil*n + 2
c
c           where nnz = original number of nonzero elements in A.
c           if iwk is not large enough the code will stop prematurely
c           with ierr = -2 or ierr = -3 (see below).
c
c On return:
c===========
c
c alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
c           the L and U factors together. The diagonal (stored in
c           alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
c           contains the i-th row of L (excluding the diagonal entry=1)
c           followed by the i-th row of U.
c
c ju      = integer array of length n containing the pointers to
c           the beginning of each row of U in the matrix alu,jlu.
c iperm   = contains the permutation arrays ..
c           iperm(1:n) = old numbers of unknowns
c           iperm(n+1:2*n) = reverse permutation = new unknowns.
c
c ierr    = integer. Error message with the following meaning.
c           ierr  = 0    --> successful return.
c           ierr .gt. 0  --> zero pivot encountered at step number ierr.
c           ierr  = -1   --> Error. input matrix may be wrong.
c                            (The elimination process has generated a
c                            row in L or U whose length is .gt.  n.)
c           ierr  = -2   --> The matrix L overflows the array al.
c           ierr  = -3   --> The matrix U overflows the array alu.
c           ierr  = -4   --> Illegal value for lfil.
c           ierr  = -5   --> zero row encountered.
c
c work arrays:
c=============
c jr,jwu,jwl 	  = integer work arrays of length n.
c wu, wl          = real work arrays of length n+1, and n resp.
c
c Notes:
c ------
c A must have all nonzero diagonal elements.
c U -- MATRIX   STORED IN UNPERMUTED FORMAT TO AVOID PERMUTATION ARRAYS
c code working -- CODED BY Y. SAAD, SEPT 9, 1993.
c-----------------------------------------------------------------------
c     local variables
c
      integer k,i,j,jrow,ju0,ii,j1,j2,jpos,len,imax,lenu,lenl,
     *     nl,jj,lenl0, lenu0,mbloc,icut
      real*8 s, tmp, tnorm,xmax,xmax0, fact, abs, t, permtol
c
      if (lfil .lt. 0) goto 998
c-------------------------------
c initialize ju0 (points to next element to be added to alu,jlu)
c and pointer.
c-----------------------------------------------------------------------
      ju0 = n+2
      jlu(1) = ju0
c
c  integer double pointer array.
c
      do 1 j=1, n
         jr(j)  = 0
         iperm(j) = j
         iperm(n+j) = j
 1    continue
c-----------------------------------------------------------------------
c  beginning of main loop.
c-----------------------------------------------------------------------
      do 500 ii = 1, n
         j1 = ia(ii)
         j2 = ia(ii+1) - 1
         tnorm = 0.0d0
         do 501 k=j1,j2
            tnorm = tnorm+abs(a(k))
 501     continue
         if (tnorm .eq. 0.0) goto 999
         tnorm = tnorm/(j2-j1+1)
c
c--- unpack L-part and U-part of row of A in arrays wl, wu --
c
         lenu = 1
         lenl = 0
         jwu(1) = ii
         wu(1) = 0.0
         jr(ii) = 1
c-----------------------------------------------------------------------
         do 170  j = j1, j2
            k = iperm(n+ja(j))
            t = a(j)
            if (abs(t) .lt. droptol*tnorm .and. k .ne. ii) goto 170
            if (k .lt. ii) then
               lenl = lenl+1
               jwl(lenl) = k
               wl(lenl) = t
               jr(k) = lenl
            else if (k .eq. ii) then
               wu(1) = t
            else
               lenu = lenu+1
               jwu(lenu) = k
               wu(lenu) = t
               jr(k) = lenu
            endif
 170     continue
c-----------------------------------------------------------------------
         tnorm = tnorm/(j2-j1+1)
         lenl0 = lenl
         lenu0 = lenu
         jj = 0
         nl = 0
c-------------------------------------------------------------------
c---------------------- eliminate previous rows --------------------
c-------------------------------------------------------------------
 150     jj = jj+1
         if (jj .gt. lenl) goto 160
c-------------------------------------------------------------------
c in order to do the elimination in the correct order we need to
c exchange the current row number with the one that has
c smallest column number, among jj,jj+1,...,lenl.
c-------------------------------------------------------------------
         jrow = jwl(jj)
         k = jj
c
c determine smallest column index
c
         do 151 j=jj+1,lenl
            if (jwl(j) .lt. jrow) then
               jrow = jwl(j)
               k = j
            endif
 151     continue
c
c     exchange in jwl
c
         if (k .ne. jj) then
            j = jwl(jj)
            jwl(jj) = jwl(k)
            jwl(k) = j
c
c     exchange in jr
c
            jr(jrow) = jj
            jr(j) = k
c
c     exchange in wl
c
            s = wl(jj)
            wl(jj) = wl(k)
            wl(k) = s
         endif
c-----------------------------------------------------------------------
         if (jrow .ge. ii) goto 160
c
c     get the multiplier for row to be eliminated: jrow
c
         fact = wl(jj)*alu(jrow)
c     zero out element in row by setting jr(jrow) = 0
c
         jr(jrow) = 0
c
         if (abs(fact)*wu(n+2-jrow) .le. droptol*tnorm) goto 150
c-------------------------------------------------------------------
c------------combine current row and row jrow ---------------------
c-------------------------------------------------------------------
         do 203 k = ju(jrow), jlu(jrow+1)-1
            s = fact*alu(k)
c     new column number
            j = iperm(n+jlu(k))
            jpos = jr(j)
c
c     if fill-in element is small then disregard:
c
            if (abs(s) .lt. droptol*tnorm .and. jpos .eq. 0) goto 203
            if (j .ge. ii) then
c
c     dealing with upper part.
c
               if (jpos .eq. 0) then
c     this is a fill-in element
                  lenu = lenu+1
                  if (lenu .gt. n) goto 995
                  jwu(lenu) = j
                  jr(j) = lenu
                  wu(lenu) = - s
               else
c     no fill-in element --
                  wu(jpos) = wu(jpos) - s
               endif
            else
c
c     dealing with lower part.
c
               if (jpos .eq. 0) then
c     this is a fill-in element
                 lenl = lenl+1
                 if (lenl .gt. n) goto 995
                 jwl(lenl) = j
                 jr(j) = lenl
                 wl(lenl) = - s
              else
c     no fill-in element --
                 wl(jpos) = wl(jpos) - s
              endif
           endif
 203	continue
        nl = nl+1
        wl(nl) = fact
        jwl(nl)  = jrow
	goto 150
c----------------------------------------------------------
c------------ update l-matrix -----------------------------
c----------------------------------------------------------
 160    len = min(nl,lenl0+lfil)
c     160    len = min0(nl,lfil)

  	call qsplit (wl,jwl,nl,len)
c
c     store L-part -- in original coordinates ..
c
        do 204 k=1, len
           if (ju0 .gt. iwk) goto 996
           alu(ju0) =  wl(k)
           jlu(ju0) = iperm(jwl(k))
c     jlu(ju0) = jwl(k)
           ju0 = ju0+1
 204    continue
c
c     save pointer to beginning of row ii of U
c
        ju(ii) = ju0
c
c     reset double-pointer jr to zero (L-part - except first
c     jj-1 elements which have already been reset)
c
	do 306 k= jj, lenl
           jr(jwl(k)) = 0
 306	continue
c----------------------------------------------------------
c------------update u-matrix -----------------------------
c----------------------------------------------------------
        len = min(lenu,lenu0+lfil+1) ! +1 added 11/23/94
c     len = min0(lenu,lfil)
	call qsplit (wu(2), jwu(2), lenu-1,len)
c
        imax = 1
        xmax = abs(wu(imax))
        xmax0 = xmax
c
        icut = ii - 1 + mbloc - mod(ii-1,mbloc)
c
        do k=2,len
           t = abs(wu(k))
           if (t .gt. xmax .and. t*permtol .gt. xmax0 .and.
     *          jwu(k) .le. icut) then
              imax = k
              xmax = t
           endif
        enddo
c
c     exchange wu's
c
        tmp = wu(1)
        wu(1) = wu(imax)
        wu(imax) = tmp
c
c     update iperm and reverse iperm
c
        j = jwu(imax)
        i = iperm(ii)
        iperm(ii) = iperm(j)
        iperm(j) = i
c     reverse iperm
        iperm(n+iperm(ii)) = ii
        iperm(n+iperm(j)) = j
c
        t = abs(wu(k))
        if (len + ju0 .gt. iwk) goto 997
c
c     store U-part in original coordinates
c
        do 302 k=2, len
           jlu(ju0) = iperm(jwu(k))
           alu(ju0) = wu(k)
           t = t + abs(wu(k) )
           ju0 = ju0+1
 302	continue
c
c     save norm in wu (backwards). Norm is in fact average abs value
c
        wu(n+2-ii) = t / (len+1)
c
c     store inverse of diagonal element of u
c
        if (wu(1) .eq. 0.0) wu(1) = (1.0D-4 + droptol)*tnorm
c
        alu(ii) = 1.0d0/ wu(1)
c
c     update pointer to beginning of next row of U.
c
	jlu(ii+1) = ju0
c
c     reset double-pointer jr to zero (U-part)
c
	do 308 k=1, lenu
           jr(jwu(k)) = 0
 308	continue
c-----------------------------------------------------------------------
c     end main loop
c-----------------------------------------------------------------------
 500  continue
c
c     permute all column indices of LU ...
c
c
c     call dvperm(n,alu,iperm(n+1))
c
c     do ii =1, n
c     do k = ju(ii), jlu(ii+1)-1
c     jlu(k) = iperm(jlu(k))
c     enddo
c     enddo
c-----------------------------------------------------------------------
      do k = jlu(1),jlu(n+1)-1
         jlu(k) = iperm(n+jlu(k))
      enddo
c
c     ...and A
c
      do k=1, ia(n+1)-1
         ja(k) = iperm(n+ja(k))
      enddo
c
      ierr = 0
      return
c
c     zero pivot :
c
c     900    ierr = ii
c     return
c
c     incomprehensible error. Matrix must be wrong.
c
 995  ierr = -1
      return
c
c     insufficient storage in L.
c
 996  ierr = -2
      return
c
c     insufficient storage in U.
c
 997  ierr = -3
      return
c
c     illegal lfil entered.
c
 998  ierr = -4
      return
c
c     zero row encountered
c
 999  ierr = -5
      return
c----------------end-of-ilutp-------------------------------------------
c-----------------------------------------------------------------------
      end
c----------------------------------------------------------------------
	subroutine ilu0 (n, a, ja, ia, alu, jlu, ju, iw, ierr)
	implicit real*8 (a-h,o-z)
	real*8 a(*), alu(*)
        integer ja(*), ia(*), ju(*), jlu(*), iw(*)
c------------------ right preconditioner ------------------------------*
c                    ***   ilu(0) preconditioner.   ***                *
c----------------------------------------------------------------------*
c Note that this has been coded in such a way that it can be used
c with pgmres. Normally, since the data structure of the L+U matrix is
c the same as that the A matrix, savings can be made. In fact with
c some definitions (not correct for general sparse matrices) all we
c need in addition to a, ja, ia is an additional diagonal.
c ILU0 is not recommended for serious problems. It is only provided
c here for comparison purposes.
c-----------------------------------------------------------------------
c
c on entry:
c---------
c n       = dimension of matrix
c a, ja,
c ia      = original matrix in compressed sparse row storage.
c
c on return:
c-----------
c alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
c           the L and U factors together. The diagonal (stored in
c           alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
c           contains the i-th row of L (excluding the diagonal entry=1)
c           followed by the i-th row of U.
c
c ju	  = pointer to the diagonal elements in alu, jlu.
c
c ierr	  = integer indicating error code on return
c	     ierr = 0 --> normal return
c	     ierr = k --> code encountered a zero pivot at step k.
c work arrays:
c-------------
c iw	    = integer work array of length n.
c------------
c IMPORTANT
c-----------
c it is assumed that the the elements in the input matrix are stored
c    in such a way that in each row the lower part comes first and
c    then the upper part. To get the correct ILU factorization, it is
c    also necessary to have the elements of L sorted by increasing
c    column number. It may therefore be necessary to sort the
c    elements of a, ja, ia prior to calling ilu0. This can be
c    achieved by transposing the matrix twice using csrcsc.
c
c-----------------------------------------------------------------------
        ju0 = n+2
        jlu(1) = ju0
c
c initialize work vector to zero's
c
	do 31 i=1, n
           iw(i) = 0
 31     continue
c
c main loop
c
	do 500 ii = 1, n
           js = ju0
c
c generating row number ii of L and U.
c
           do 100 j=ia(ii),ia(ii+1)-1
c
c     copy row ii of a, ja, ia into row ii of alu, jlu (L/U) matrix.
c
              jcol = ja(j)
              if (jcol .eq. ii) then
                 alu(ii) = a(j)
                 iw(jcol) = ii
                 ju(ii)  = ju0
              else
                 alu(ju0) = a(j)
                 jlu(ju0) = ja(j)
                 iw(jcol) = ju0
                 ju0 = ju0+1
              endif
 100       continue
           jlu(ii+1) = ju0
           jf = ju0-1
           jm = ju(ii)-1
c
c     exit if diagonal element is reached.
c
           do 150 j=js, jm
              jrow = jlu(j)
              tl = alu(j)*alu(jrow)
              alu(j) = tl
c
c     perform  linear combination
c
              do 140 jj = ju(jrow), jlu(jrow+1)-1
                 jw = iw(jlu(jj))
                 if (jw .ne. 0) alu(jw) = alu(jw) - tl*alu(jj)
 140          continue
 150       continue
c
c     invert  and store diagonal element.
c
           if (alu(ii) .eq. 0.0d0) goto 600
           alu(ii) = 1.0d0/alu(ii)
c
c     reset pointer iw to zero
c
           iw(ii) = 0
           do 201 i = js, jf
 201          iw(jlu(i)) = 0
 500       continue
           ierr = 0
           return
c
c     zero pivot :
c
 600       ierr = ii
c
           return
c------- end of ilu0 ---------------------------------------------------
c-----------------------------------------------------------------------
           end
c----------------------------------------------------------------------
	subroutine milu0 (n, a, ja, ia, alu, jlu, ju, iw, ierr)
	implicit real*8 (a-h,o-z)
	real*8 a(*), alu(*)
	integer ja(*), ia(*), ju(*), jlu(*), iw(*)
c----------------------------------------------------------------------*
c                *** simple milu(0) preconditioner. ***                *
c----------------------------------------------------------------------*
c Note that this has been coded in such a way that it can be used
c with pgmres. Normally, since the data structure of a, ja, ia is
c the same as that of a, ja, ia, savings can be made. In fact with
c some definitions (not correct for general sparse matrices) all we
c need in addition to a, ja, ia is an additional diagonal.
c Ilu0 is not recommended for serious problems. It is only provided
c here for comparison purposes.
c-----------------------------------------------------------------------
c
c on entry:
c----------
c n       = dimension of matrix
c a, ja,
c ia      = original matrix in compressed sparse row storage.
c
c on return:
c----------
c alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
c           the L and U factors together. The diagonal (stored in
c           alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
c           contains the i-th row of L (excluding the diagonal entry=1)
c           followed by the i-th row of U.
c
c ju	  = pointer to the diagonal elements in alu, jlu.
c
c ierr	  = integer indicating error code on return
c	     ierr = 0 --> normal return
c	     ierr = k --> code encountered a zero pivot at step k.
c work arrays:
c-------------
c iw	    = integer work array of length n.
c------------
c Note (IMPORTANT):
c-----------
C it is assumed that the the elements in the input matrix are ordered
c    in such a way that in each row the lower part comes first and
c    then the upper part. To get the correct ILU factorization, it is
c    also necessary to have the elements of L ordered by increasing
c    column number. It may therefore be necessary to sort the
c    elements of a, ja, ia prior to calling milu0. This can be
c    achieved by transposing the matrix twice using csrcsc.
c-----------------------------------------------------------
          ju0 = n+2
          jlu(1) = ju0
c initialize work vector to zero's
	do 31 i=1, n
 31           iw(i) = 0
c
c-------------- MAIN LOOP ----------------------------------
c
	do 500 ii = 1, n
           js = ju0
c
c generating row number ii or L and U.
c
           do 100 j=ia(ii),ia(ii+1)-1
c
c copy row ii of a, ja, ia into row ii of alu, jlu (L/U) matrix.
c
              jcol = ja(j)
              if (jcol .eq. ii) then
                 alu(ii) = a(j)
                 iw(jcol) = ii
                 ju(ii)  = ju0
              else
                 alu(ju0) = a(j)
                 jlu(ju0) = ja(j)
                 iw(jcol) = ju0
                 ju0 = ju0+1
              endif
 100       continue
           jlu(ii+1) = ju0
           jf = ju0-1
           jm = ju(ii)-1
c s accumulates fill-in values
           s = 0.0d0
           do 150 j=js, jm
              jrow = jlu(j)
              tl = alu(j)*alu(jrow)
              alu(j) = tl
c-----------------------perform linear combination --------
              do 140 jj = ju(jrow), jlu(jrow+1)-1
                 jw = iw(jlu(jj))
                 if (jw .ne. 0) then
                       alu(jw) = alu(jw) - tl*alu(jj)
                    else
                       s = s + tl*alu(jj)
                    endif
 140          continue
 150       continue
c----------------------- invert and store diagonal element.
           alu(ii) = alu(ii)-s
           if (alu(ii) .eq. 0.0d0) goto 600
           alu(ii) = 1.0d0/alu(ii)
c----------------------- reset pointer iw to zero
           iw(ii) = 0
           do 201 i = js, jf
 201          iw(jlu(i)) = 0
 500       continue
           ierr = 0
           return
c     zero pivot :
 600       ierr = ii
           return
c------- end of milu0 --------------------------------------------------
c-----------------------------------------------------------------------
           end
c-----------------------------------------------------------------------
	subroutine lusol (n, y, x, alu, jlu, ju)
        real*8 x(n), y(n), alu(*)
	integer n, jlu(*), ju(*)
c-----------------------------------------------------------------------
c
c performs a forward followed by a backward solve
c for LU matrix as produced by  ILUT
c
c-----------------------------------------------------------------------
c local variables
c
        integer i,k
c
c forward solve
c
        do 40 i = 1, n
           x(i) = y(i)
           do 41 k=jlu(i),ju(i)-1
              x(i) = x(i) - alu(k)* x(jlu(k))
 41        continue
 40     continue
c
c     backward solve.
c
	do 90 i = n, 1, -1
	   do 91 k=ju(i),jlu(i+1)-1
              x(i) = x(i) - alu(k)*x(jlu(k))
 91	   continue
           x(i) = alu(i)*x(i)
 90     continue
c
  	return
c----------------end of lusol ------------------------------------------
c-----------------------------------------------------------------------
	end
c-----------------------------------------------------------------------
	subroutine lutsol (n, y, x, alu, jlu, ju) 
        real*8 x(n), y(n), alu(*)
	integer n, jlu(*), ju(*)
c-----------------------------------------------------------------------
c
c Given a LU decomposition of a matrix stored in (alu, jlu, ju),
c this routine performs x = (LU)^{-T} y, i.e. LU transposed solve.
c 
c-----------------------------------------------------------------------
c local variables
c
        integer i,k
c
        do 10 i = 1, n
           x(i) = y(i)
 10     continue
c
c forward solve (with U^T)
c
        do 20 i = 1, n
           x(i) = x(i) * alu(i)
           do 30 k=ju(i),jlu(i+1)-1
              x(jlu(k)) = x(jlu(k)) - alu(k)* x(i)
 30        continue
 20     continue
c     
c     backward solve (with L^T)
c     
	do 40 i = n, 1, -1 
	   do 50 k=jlu(i),ju(i)-1
              x(jlu(k)) = x(jlu(k)) - alu(k)*x(i)
 50        continue
 40     continue
c
  	return
c----------------end of lutsol -----------------------------------------
c-----------------------------------------------------------------------
	end
c----------------------------------------------------------------------- 
        subroutine qsplit  (a, ind, n, ncut)
        real*8 a(n)
        integer ind(n), n, ncut
c-----------------------------------------------------------------------
c     does a quick-sort split of a real array.
c     on input a(1:n). is a real array
c     on output a(1:n) is permuted such that its elements satisfy:
c
c     abs(a(i)) .ge. abs(a(ncut)) for i .lt. ncut and
c     abs(a(i)) .le. abs(a(ncut)) for i .gt. ncut
c
c     ind(1:n) is an integer array which permuted in the same way as a(*).
c-----------------------------------------------------------------------
        real*8 tmp, abskey
        integer itmp, first, last
c-----
        first = 1
        last = n
        if (ncut .lt. first .or. ncut .gt. last) return
c
c     outer loop -- while mid .ne. ncut do
c
 1      mid = first
        abskey = abs(a(mid))
        do 2 j=first+1, last
           if (abs(a(j)) .gt. abskey) then
              mid = mid+1
c     interchange
              tmp = a(mid)
              itmp = ind(mid)
              a(mid) = a(j)
              ind(mid) = ind(j)
              a(j)  = tmp
              ind(j) = itmp
           endif
 2      continue
c
c     interchange
c
        tmp = a(mid)
        a(mid) = a(first)
        a(first)  = tmp
c
        itmp = ind(mid)
        ind(mid) = ind(first)
        ind(first) = itmp
c
c     test for while loop
c
        if (mid .eq. ncut) return
        if (mid .gt. ncut) then
           last = mid-1
        else
           first = mid+1
        endif
        goto 1
c----------------end-of-qsplit------------------------------------------
c-----------------------------------------------------------------------
        end
