subroutine spars_solve (nnn,nnz,ia,ja,sam,srhs,sp)
c
c=======================================================================
c This subroutine is a front end to the SPARSKIT package.            ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <iounits.h>
#include <convinfo.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,ierr,lfil,nnn,nnz,nrow,nwk,nzmin
#ifdef fixprecond
     &        ,fxpcnd,ncall
#endif
      integer ia(nnz),ipar(16),iw(nmax*2),ja(nnz),jau(nzmax),ju(nmax)
      logical bad
      FLOAT
     &      sam(nnz),srhs(nnn),sp(nnn)
#ifdef iforttime
      real dum
      real dtime
#endif
      double precision tol
      double precision am(nzmax),au(nzmax),fpar(16),p(nmax),rhs(nmax),
     &                 sol(nmax),wk(nmax*40),xran(nmax)
c
      external cg, runrc, ilut
c
#ifdef fixprecond
      save fxpcnd,ncall,au,jau,ju
c
      data fxpcnd,ncall /5,0/
c
#endif
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Check storage space.
c-----------------------------------------------------------------------
c
      bad   = .false.
      nzmin = nnz
      if (cgstat.ne.0) nzmin=max( nzmin, (nmax*2*(fillin+1)) )
c
      if (nnn.gt.nmax) then
         write (stdout,*) 'Stopping in spars_solve : set NMAX to at ',
     &                    'least ',nnn,' in param.h'
         bad = .true.
         if (cgstat.ne.0) nzmin=max( nzmin, (nnn*2*(fillin+1)) )
      endif
      if (nzmax.lt.nzmin) then
         write (stdout,*) 'Stopping in spars_solve : set NZMAX to at ',
     &                    'least ',nzmin,' in param.h'
         bad = .true.
      endif
c
      if (bad) call exitus ('SPARS_SOLVE')
c
c-----------------------------------------------------------------------
c  pass single to double precision
c-----------------------------------------------------------------------
c
      do 10 i=1,nnz
         am(i)=dble(sam(i))
  10  continue
      do 20 i=1,nnn
         rhs(i)=dble(srhs(i))
         p(i)=dble(sp(i))
  20  continue
c 
c-----------------------------------------------------------------------
c  set the parameters for the iterative solvers
c-----------------------------------------------------------------------
c
      ipar(2) = cgstat
      ipar(3) = 1
      ipar(4) = lwk
      ipar(5) = 16
      ipar(6) = maxits
      fpar(1) = tolrel
      fpar(2) = tolabs
c
      nrow = nnn
c
c-----------------------------------------------------------------------
c  set-up the preconditioner ILUT(15, 1E-4)  ! new definition of lfil
c-----------------------------------------------------------------------
#ifdef timesolver
c
# ifndef iforttime
      call dtime(dttpcg)
# else
      dum = dtime(dttpcg)
# endif
#endif
c
#ifndef fixprecond
      if(ipar(2).ne.0) then
#else
      if ((ipar(2).ne.0).and.(ncall.lt.fxpcnd)) then
         ncall = ncall + 1
#endif
         lfil = fillin
         tol = tolpcg
         nwk = nzmax
         call ilut (nrow,am,ja,ia,lfil,tol,au,jau,ju,nwk,
     *     wk,iw,ierr)
         if(ierr.ne.0) then
            write (stdout,900) ierr
            if (ierr.gt.0) then
               write (stdout,910) ierr
            elseif (ierr.eq.-1) then
               write (stdout,920) nrow
            elseif (ierr.eq.-2) then
               write (stdout,930)
            elseif (ierr.eq.-3) then
               write (stdout,940)
            elseif (ierr.eq.-4) then
               write (stdout,950) lfil
            elseif (ierr.eq.-5) then
               write (stdout,960)
            else
               write (stdout,970)
            endif
            call exitus('SPARS_SOLVE')
         endif
      endif
c
c-----------------------------------------------------------------------
c  set initial guess
c-----------------------------------------------------------------------
c
      do 30 i = 1, nrow
         xran(i) = p(i)
         sol(i)=0.d0
  30  continue
c
c-----------------------------------------------------------------------
c  Solve sparse matrix.
c-----------------------------------------------------------------------
c
      call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,am,ja,ia,au,jau,ju,
     +     mican,cg,ierr,stdout)
c
#ifdef timesolver
# ifndef iforttime
      call dtime(dttpcg)
# else
      dum = dtime(dttpcg)
# endif
c
#endif
c-----------------------------------------------------------------------
c  Pass to single precision and check for errors.
c-----------------------------------------------------------------------
c
      do 40 i = 1, nrow
         sp(i) = real(sol(i))
  40  continue
c
      resini=real(fpar(3))
      restarget=real(fpar(4))
      residu=real(fpar(5))
      if(ierr.ne.0) then
        call exitus('SPARS_SOLVE')
      endif
c
      return
c
 900  format (/'***Error:  SPARS_SOLVE - error flag returned from ILUT;'
     &        ' IERR = ',i10/11x,'Probable meaning:')
 910  format (11x,'zero pivot encountered at step number ',i10)
 920  format (11x,'Input matrix may be wrong.  (The elimination process'
     &        ' has generated a'/11x,'row in L or U whose length is > ',
     &        'N=',i10,'.)')
 930  format (11x,'The matrix L overflows the array AL.')
 940  format (11x,'The matrix U overflows the array ALU.')
 950  format (11x,'Illegal value for LFIL =',i10)
 960  format (11x,'zero row encountered.')
 970  format (11x,'unknown error.')
c
      end