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