      SUBROUTINE NISLVI(N,NFMAX,IOPT,IFAIL,LIWK,IWK,LPIWK,LTIWK,
     $                  LRWK,RWK,LPRWK,LTRWK)
C*    Begin Prologue SOLINI
      INTEGER N,NFMAX
      INTEGER IOPT(50)
      INTEGER IFAIL
      INTEGER LIWK
      INTEGER IWK(LIWK)
      INTEGER LPIWK,LTIWK
      INTEGER LRWK
      DOUBLE PRECISION RWK(LRWK)
      INTEGER LPRWK,LTRWK
C     ------------------------------------------------------------
C
C*    Summary :
C
C     S O L I N I : Initialize linear algebra subprograms for 
C                   factorization of a sparse (N,N)-matrix
C
C*    Input parameters (* marks inout parameters)
C     ===========================================
C
C     N             Int    Order of the linear system
C     NFMAX         Int    Dimension of the matrix array A -
C                          the maximum number of nonzero elements
C                          of the sparse matrix, NFMAX <= 32767
C     IOPT(50)      Int    Option vector passed from NLEQ1S
C
C*    Output parameters
C     =================
C
C     IFAIL         Int    Error indicator returned by this routine:
C                          = 0 matrix decomposition successfull
C                          =10 supplied workspace too small
C
C*    Workspace parameters
C     ====================
C
C     LIWK          Int    Length of integer workspace (In)
C     IWK(LIWK)     Int    Integer Workspace supplied for linear solver
C     LPIWK         Int    Length of integer Workspace used by 
C                          linear solver across successive calls - i.e.
C                          'permanent' (out)       
C     LTIWK         Int    Length of integer Workspace used by 
C                          linear solver - but not across successive
C                          calls - i.e. 'temporary' (out)       
C     LRWK          Int    Length of real workspace (in)      
C     RWK(LRWK)     Dble   Real Workspace supplied for linear solver
C     LPRWK         Int    Length of real Workspace used by 
C                          linear solver across successive calls - i.e.
C                          'permanent' (out)
C     LTRWK         Int    Length of real Workspace used by 
C                          linear solver - but not across successive
C                          calls - i.e. 'temporary' (out)
C
C*    Subroutines called:  none
C
C     ------------------------------------------------------------
C*    End Prologue
      INTEGER IRM1,IRM2
      LOGICAL QSUCC
C*    Begin
      IFAIL = 0
      LPIWK = 8*NFMAX+11*N+12
      LTIWK = 0
      LPRWK = 5*NFMAX+2*N+8
      LTRWK = 0
      NN    = 5*NFMAX
      NN1   = 3*NFMAX
      IRM2  = LIWK-(LPIWK+LTIWK)
      IRM1  = MIN0(IRM2, LRWK-(LPRWK+LTRWK))
      IF (IRM1.LT.0) THEN
        IFAIL=10
        IWK(1)=-IRM1
      ELSE
        NN  = NN + IRM1
        LPRWK = LPRWK + IRM1
        LPIWK = LPIWK + IRM1
        NN1 = NN1 + (IRM2-IRM1)
        LPIWK = LPIWK + (IRM2-IRM1)
      ENDIF
      IWK(LIWK-1) = NN
      IWK(LIWK)   = NN1
C     ( Real workspace starts at RWK(NRWKFR) )
      IF (LIWK.LT.LPIWK+LTIWK.OR.LRWK.LT.LPRWK+LTRWK
     $   .OR.IFAIL.NE.0) IFAIL=10
      RETURN
      END
C
      SUBROUTINE NIFACT(N,NFILL,A,IROW,ICOL,IOPT,QSTRUC,IFAIL,LIWK,IWK,
     $LRWK,RWK)
C*    Begin Prologue FACT
      INTEGER N,NFILL
      DOUBLE PRECISION A(NFILL)
      INTEGER IROW(NFILL),ICOL(NFILL)
      INTEGER IOPT(50)
      LOGICAL QSTRUC
      INTEGER IFAIL
      INTEGER LIWK
      INTEGER IWK(LIWK)
      INTEGER LRWK
      DOUBLE PRECISION RWK(LRWK)
C     ------------------------------------------------------------
C
C*    Summary :
C
C     F A C T : Call linear algebra subprogram for factorization of
C               a (N,N)-matrix
C
C*    Input parameters
C     ================
C
C     N             Int    Order of the linear system
C     NFILL         Int    Number of nonzero elements in the sparse
C                          matrix, NFILL <= NFMAX
C     A(NFMAX)      Dble   Matrix real evalues storage. 
C                          See routine NIINT.
C     IROW(NFMAX)   ShInt  Matrix row indices storage.
C     ICOL(NFMAX)   ShInt  Matrix column indices storage.
C     IOPT(50)      Int    Option vector passed from NLEQ1S
C     QSTRUC        Logic  (Input) Tells, if the sparse structure
C                          of the current matrix is different from 
C                          that of the matrix supplied in the previous
C                          call:
C                          = .TRUE.  : Structure is different
C                          = .FALSE. : Structure is the same as before
C
C*    Output parameters
C     =================
C
C     IFAIL         Int    Error indicator returned by this routine:
C                          = 0 matrix decomposition successfull
C                          = 1 singular matrix
C                          other : See IFLAG values of Y12MxF routines.
C
C*    Workspace parameters
C     ====================
C
C     LIWK          Int    Length of integer workspace (In)
C     IWK(LIWK)     Int    Integer Workspace supplied for linear solver
C     LRWK          Int    Length of real workspace (in)      
C     RWK(LRWK)     Dble   Real Workspace supplied for linear solver
C
C*    Subroutines called:  Y12MBF, Y12MCF
C
C     References:
C
CÊ    /1/ Z. Zlatev:
C         On some pivotal strategies in Gaussian elimination by 
C         sparse technique.
C         SIAM Journal on Numerical Analysis, 17 (1980), 18-30.
C
CÊ    /2/ Z. Zlatev:
C         Computational Methods for General Sparse Matrices.
C         Kluwer Academic Publishers, Dordrecht, Boston, London, 1991.
C
C     /3/ Z. Zlatev and I. Dimov:
C         Computational and Environmental Challenges in Environmental
C         Modelling.
C         Elsevier, Amsterdam, Boston, Heidelberg, London, New York,
C         Oxford, Paris, San Diego, San Francisco, Singapore, Sidney,
C         Tokyo, 2006 (Chapter 6).
C
C     /4/ Z. Zlatev, J. Wasniewsky and K. Shaumburg:
C         Y12M Ð Solution of Large and Sparse Systems of 
C         Linear Algebraic Equations.
C         Spinger-Verlag, Berlin, Heidelberg, New York, 1981.
C
C     ------------------------------------------------------------
C*    End Prologue
C*    Begin
      MPRERR = IOPT(11)
      LUERR = IOPT(12)
      IFAIL = 0
      NN  = IWK(LIWK-1)
      NN1 = IWK(LIWK)
      IHA = N*11
      DO I=1,NFILL
         RWK(I)    = A(I)
         IWK(NN+I) = IROW(I)
         IWK(I)    = ICOL(I)
      ENDDO
      IF (QSTRUC) THEN
         IWK(NN+NN1+IHA+4) = 1 ! IFLAG(4)
      ELSE
         IWK(NN+NN1+IHA+4) = 2 ! IFLAG(4)
      ENDIF
      DO I=1,8
        RWK(NN+I) = 0.0D0
      ENDDO
      CALL Y12MBF(N, NFILL, RWK(1), IWK(1), NN, IWK(NN+1), NN1,
     1            IWK(NN+NN1+1), N, RWK(NN+1), IWK(NN+NN1+IHA+1),
     2            IFAIL)
      IF (IFAIL.NE.0) THEN
         IF (MPRERR.GT.0) WRITE(LUERR,*) "Y12MBF FAILED"
         RETURN
      ENDIF
      RWK(NN+1) = 9.0      ! AFLAG(1) = STABILITY FACTOR
      RWK(NN+2) = 1.0D-12  ! AFLAG(2) = DROP TOLERANCE
      RWK(NN+3) = 1.0D16   ! AFLAG(3) = STOP GROWTH FACTOR
      RWK(NN+4) = 1.0D-12
      DO I=1,N
        RWK(NN+8+N+I)=0.0D0
      ENDDO
      IWK(NN+NN1+IHA+2) = 3 ! IFLAG(2)
      IWK(NN+NN1+IHA+3) = 1 ! IFLAG(3)  general pivot search
      IWK(NN+NN1+IHA+5) = 2 ! IFLAG(5)
      CALL Y12MCF(N, NFILL, RWK(1), IWK(1), NN, IWK(NN+1), NN1,
     1            RWK(NN+9), RWK(NN+9+N), IWK(NN+NN1+1), N, RWK(NN+1),
     2            IWK(NN+NN1+IHA+1), IFAIL)
      IF (IFAIL.NE.0 .AND. MPRERR.GT.0) WRITE(LUERR,*) "Y12MCF FAILED"
      RETURN
      END
C
      SUBROUTINE NISOLV(N,NFILL,A,IROW,ICOL,B,IOPT,IFAIL,LIWK,IWK,
     $LRWK,RWK)
C*    Begin Prologue SOLVE
      INTEGER N,NFILL
      DOUBLE PRECISION A(NFILL)
      INTEGER IROW(NFILL),ICOL(NFILL)
      DOUBLE PRECISION B(N)
      INTEGER IOPT(50)
      INTEGER IFAIL
      INTEGER LIWK
      INTEGER IWK(LIWK)
      INTEGER LRWK
      DOUBLE PRECISION RWK(LRWK)
C     ------------------------------------------------------------
C
C*    Summary :
C
C     S O L V E : Call linear algebra subprogram for solution of
C                 the linear system A*Z = B
C
C*    Parameters
C     ==========
C
C     N,NFILL,A,IROW,ICOL,IOPT,IFAIL,LIWK,IWK,LRWK,RWK :
C                        See description for subroutine NIFACT.          
C     B(N)       Dble    In:  Right hand side of the linear system
C                        Out: Solution of the linear system
C
C     Subroutines called: Y12MDF
C
C     ------------------------------------------------------------
C*    End Prologue
C*    Begin
      NN  = IWK(LIWK-1)
      NN1 = IWK(LIWK)
      IHA = N*11
      IWK(NN+NN1+IHA+5)=3 
      CALL Y12MDF(N, RWK(1), NN, B, RWK(NN+9), IWK(1),
     1            IWK(NN+NN1+1), N, IWK(NN+NN1+IHA+1), IFAIL)      
      RETURN
      END
C
