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 LVALUE, LINDEX
      LOGICAL QSUCC
C*    Begin
      MPRERR = IOPT(11)
      LUERR = IOPT(12)
      IFAIL = 0
C     Minimum required storage length
      LVALUE = 2*NFMAX
      LINDEX = 3*NFMAX+2*N+1
C
      LPIWK = 82+LINDEX
      LTIWK = 0
      LPRWK = 30+LVALUE
      LTRWK = 5*N
C     
      IF (LIWK.LT.LPIWK+LTIWK.OR.LRWK.LT.LPRWK+LTRWK
     $   .OR.IFAIL.NE.0) IFAIL=10
C     Increase matrix storage to maximum available length
      LVALUE = LVALUE + LRWK - LPRWK - LTRWK
      LINDEX = LINDEX + LIWK - LPIWK - LTIWK
      LPIWK = 82+LINDEX
      LPRWK = 30+LVALUE
C     Initialize MA38
      CALL MA38ID(IWK(1), RWK(1), IWK(21))
      IWK(21) = LUERR
      IWK(22) = LUERR
      IWK(41) = LVALUE
      IWK(42) = LINDEX
      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
      INTEGER I
C*    Begin
      MPRERR = IOPT(11)
      LUERR = IOPT(12)
      IFAIL = 0
      DO 10 I=1,NFILL
         IWK(82+I) = IROW(I)
         IWK(82+NFILL+I) = ICOL(I)
         RWK(30+I) = A(I)
10    CONTINUE
      IF (QSTRUC) THEN
         CALL MA38AD(N,NFILL,1,.FALSE.,IWK(41),IWK(42),RWK(31),IWK(83), 
     $               IWK(1),RWK(1),IWK(21),IWK(43),RWK(11))
      ELSE
         CALL MA38BD(N,NFILL,1,.FALSE.,IWK(41),IWK(42),RWK(31),IWK(83), 
     $               IWK(1),RWK(1),IWK(21),IWK(43),RWK(11)) 
      ENDIF
      IFAIL = IWK(43)
      IF (IFAIL.EQ.1) IFAIL=-9
      IF (IFAIL.EQ.4) IFAIL=1
      IF (IFAIL.NE.0) WRITE(LUERR,10001) IWK(41),IWK(42)
      RETURN
10001 FORMAT(1X,'Current sizes are:'/,1X,'LVALUE = ',I10,'  LINDEX = ',
     $       I10)
      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
      INTEGER I,IB
C*    Begin
      CALL MA38CD(N, 0,.FALSE.,IWK(41),IWK(42),RWK(31),IWK(83),IWK(1), 
     $            B,RWK(31+IWK(41)),RWK(31+IWK(41)+N),RWK(1),IWK(21),
     $            IWK(43),RWK(11)) 
      IB = 30+IWK(41)
      DO 10 I=1,N
        B(I)=RWK(IB+I)
10    CONTINUE
      RETURN
      END
C