sd_solve_1.f

Go to the documentation of this file.

!                   *********************
                    SUBROUTINE sd_solve_1
!                   *********************
!
     &(npoin,nsegb,gloseg,maxseg,da,xa,xinc,rhs,infogr,typext)
!
!***********************************************************************
! BIEF   V7P0                                    21/07/2011
!***********************************************************************
!
!brief    DIRECT RESOLUTION OF A SYMMETRICAL LINEAR SYSTEM WITH
!+                MINIMUM DEGREE PERMUTATION AND LDLT DECOMPOSITION.
!+
!+            FROM SEGMENT STORAGE TO COMPACT STORAGE (MORSE).
!code
!+ IMPORTANT NOTE: INSPIRED FROM PACKAGE CMLIB3 - YALE UNIVERSITE-YSMP
!+
!+
!+               YALE SPARSE MATRIX PACKAGE - NONSYMMETRIC CODES
!+                    SOLVING THE SYSTEM OF EQUATIONS MX = B
!+                        (UNCOMPRESSED POINTER STORAGE)
!+
!+    I.   CALLING SEQUENCES
!+         THE COEFFICIENT MATRIX CAN BE PROCESSED BY AN ORDERING ROUTINE
!+    (E.G., TO REDUCE FILLIN OR ENSURE NUMERICAL STABILITY) BEFORE USING
!+    THE REMAINING SUBROUTINES.  IF NO REORDERING IS DONE, THEN SET
!+    R(I) = C(I) = IC(I) = I  FOR I=1,...,N.  THE CALLING SEQUENCE IS --
!+        (      (MATRIX ORDERING))
!+         NSF   (SYMBOLIC FACTORIZATION TO DETERMINE WHERE FILLIN WILL
!+                 OCCUR DURING NUMERIC FACTORIZATION)
!+         NNF   (NUMERIC FACTORIZATION INTO PRODUCT LDU OF UNIT LOWER
!+                 TRIANGULAR MATRIX L, DIAGONAL MATRIX D, AND UNIT UPPER
!+                 TRIANGULAR MATRIX U, AND SOLUTION OF LINEAR SYSTEM)
!+         NNS   (SOLUTION OF LINEAR SYSTEM FOR ADDITIONAL RIGHT-HAND
!+     OR          SIDE USING LDU FACTORIZATION FROM NNF)
!+         NNT   (SOLUTION OF TRANSPOSED LINEAR SYSTEM FOR ADDITIONAL
!+                 RIGHT-HAND SIDE USING LDU FACTORIZATION FROM NNF)
!+
!+    II.  STORAGE OF SPARSE MATRICES
!+         THE NONZERO ENTRIES OF THE COEFFICIENT MATRIX M ARE STORED
!+    ROW-BY-ROW IN_SS1 THE ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO
!+    ENTRIES IN_SS1 EACH ROW, WE NEED TO KNOW IN_SS1 WHICH COLUMN EACH ENTRY
!+    LIES.  THE COLUMN INDICES WHICH CORRESPOND TO THE NONZERO ENTRIES
!+    OF M ARE STORED IN_SS1 THE ARRAY JA;  I.E., IF  A(K) = M(I,J),  THEN
!+    JA(K) = J.  IN_SS1 ADDITION, WE NEED TO KNOW WHERE EACH ROW STARTS AND
!+    HOW LONG IT IS.  THE INDEX POSITIONS IN_SS1 JA AND A WHERE THE ROWS OF
!+    M BEGIN ARE STORED IN_SS1 THE ARRAY IA;  I.E., IF M(I,J) IS THE FIRST
!+    NONZERO ENTRY (STORED) IN_SS1 THE I-TH ROW AND A(K) = M(I,J),  THEN
!+    IA(I) = K.  MOREOVER, THE INDEX IN_SS1 JA AND A OF THE FIRST LOCATION
!+    FOLLOWING THE LAST ELEMENT IN_SS1 THE LAST ROW IS STORED IN_SS1 IA(N+1).
!+    THUS, THE NUMBER OF ENTRIES IN_SS1 THE I-TH ROW IS GIVEN BY
!+    IA(I+1) - IA(I),  THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED
!+    CONSECUTIVELY IN_SS1
!+            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),
!+    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN_SS1
!+            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).
!+    FOR EXAMPLE, THE 5 BY 5 MATRIX
!+                ( 1. 0. 2. 0. 0.)
!+                ( 0. 3. 0. 0. 0.)
!+            M = ( 0. 4. 5. 6. 0.)
!+                ( 0. 0. 0. 7. 0.)
!+                ( 0. 0. 0. 8. 9.)
!+    WOULD BE STORED AS
!+               \ 1  2  3  4  5  6  7  8  9
!+            ---+--------------------------
!+            IA \ 1  3  4  7  8 10
!+            JA \ 1  3  2  2  3  4  4  4  5
!+             A \ 1. 2. 3. 4. 5. 6. 7. 8. 9.         .
!+
!+         THE STRICT TRIANGULAR PORTIONS OF THE MATRICES L AND U ARE
!+    STORED IN_SS1 THE SAME FASHION USING THE ARRAYS  IL, JL, L  AND
!+    IU, JU, U  RESPECTIVELY.  THE DIAGONAL ENTRIES OF L AND U ARE
!+    ASSUMED TO BE EQUAL TO ONE AND ARE NOT STORED.  THE ARRAY D
!+    CONTAINS THE RECIPROCALS OF THE DIAGONAL ENTRIES OF THE MATRIX D.
!+
!+    III. ADDITIONAL STORAGE SAVINGS
!+         IN_SS1 NSF, R AND IC CAN BE THE SAME ARRAY IN_SS1 THE CALLING
!+    SEQUENCE IF NO REORDERING OF THE COEFFICIENT MATRIX HAS BEEN DONE.
!+         IN_SS1 NNF, R, C AND IC CAN ALL BE THE SAME ARRAY IF NO REORDERING
!+    HAS BEEN DONE.  IF ONLY THE ROWS HAVE BEEN REORDERED, THEN C AND IC
!+    CAN BE THE SAME ARRAY.  IF THE ROW AND COLUMN ORDERINGS ARE THE
!+    SAME, THEN R AND C CAN BE THE SAME ARRAY.  Z AND ROW CAN BE THE
!+    SAME ARRAY.
!+         IN_SS1 NNS OR NNT, R AND C CAN BE THE SAME ARRAY IF NO REORDERING
!+    HAS BEEN DONE OR IF THE ROW AND COLUMN ORDERINGS ARE THE SAME.  Z
!+    AND B CAN BE THE SAME ARRAY;  HOWEVER, THEN B WILL BE DESTROYED.
!+
!+    IV.  PARAMETERS
!+         FOLLOWING IS A LIST OF PARAMETERS TO THE PROGRAMS.  NAMES ARE
!+    UNIFORM AMONG THE VARIOUS SUBROUTINES.  CLASS ABBREVIATIONS ARE --
!+       N - INTEGER VARIABLE
!+       F - REAL VARIABLE
!+       V - SUPPLIES A VALUE TO A SUBROUTINE
!+       R - RETURNS A RESULT FROM A SUBROUTINE
!+       I - USED INTERNALLY BY A SUBROUTINE
!+       A - ARRAY
!+
!+ CLASS \ PARAMETER
!+ ------+----------
!+ FVA   \ A     - NONZERO ENTRIES OF THE COEFFICIENT MATRIX M, STORED
!+       \           BY ROWS.
!+       \           SIZE = NUMBER OF NONZERO ENTRIES IN_SS1 M.
!+ FVA   \ B     - RIGHT-HAND SIDE B.
!+       \           SIZE = N.
!+ NVA   \ C     - ORDERING OF THE COLUMNS OF M.
!+       \           SIZE = N.
!+ FVRA  \ D     - RECIPROCALS OF THE DIAGONAL ENTRIES OF THE MATRIX D.
!+       \           SIZE = N.
!+ NR    \ FLAG  - ERROR FLAG;  VALUES AND THEIR MEANINGS ARE --
!+       \            0     NO ERRORS DETECTED
!+       \            N+K   NULL ROW IN_SS1 A  --  ROW = K
!+       \           2N+K   DUPLICATE ENTRY IN_SS1 A  --  ROW = K
!+       \           3N+K   INSUFFICIENT STORAGE FOR JL  --  ROW = K
!+       \           4N+1   INSUFFICIENT STORAGE FOR L
!+       \           5N+K   NULL PIVOT  --  ROW = K
!+       \           6N+K   INSUFFICIENT STORAGE FOR JU  --  ROW = K
!+       \           7N+1   INSUFFICIENT STORAGE FOR U
!+       \           8N+K   ZERO PIVOT  --  ROW = K
!+ NVA   \ IA    - POINTERS TO DELIMIT THE ROWS IN_SS1 A.
!+       \           SIZE = N+1.
!+ NVA   \ IC    - INVERSE OF THE ORDERING OF THE COLUMNS OF M;  I.E.,
!+       \           IC(C(I) = I  FOR I=1,...N.
!+       \           SIZE = N.
!+ NVRA  \ IL    - POINTERS TO DELIMIT THE ROWS IN_SS1 L.
!+       \           SIZE = N+1.
!+ NVRA  \ IU    - POINTERS TO DELIMIT THE ROWS IN_SS1 U.
!+       \           SIZE = N+1.
!+ NVA   \ JA    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A.
!+       \           SIZE = SIZE OF A.
!+ NVRA  \ JL    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF L.
!+       \           SIZE = JLMAX.
!+ NV    \ JLMAX - DECLARED DIMENSION OF JL;  JLMAX MUST BE LARGER THAN
!+       \           THE NUMBER OF NONZERO ENTRIES IN_SS1 THE STRICT LOWER
!+       \           TRIANGLE OF M PLUS FILLIN  (IL(N+1)-1 AFTER NSF).
!+ NVRA  \ JU    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF U.
!+       \           SIZE = JUMAX.
!+ NV    \ JUMAX - DECLARED DIMENSION OF JU;  JUMAX MUST BE LARGER THAN
!+       \           THE NUMBER OF NONZERO ENTRIES IN_SS1 THE STRICT UPPER
!+       \           TRIANGLE OF M PLUS FILLIN  (IU(N+1)-1 AFTER NSF).
!+ FVRA  \ L     - NONZERO ENTRIES IN_SS1 THE STRICT LOWER TRIANGULAR PORTION
!+       \           OF THE MATRIX L, STORED BY ROWS.
!+       \           SIZE = LMAX
!+ NV    \ LMAX  - DECLARED DIMENSION OF L;  LMAX MUST BE LARGER THAN
!+       \           THE NUMBER OF NONZERO ENTRIES IN_SS1 THE STRICT LOWER
!+       \           TRIANGLE OF M PLUS FILLIN  (IL(N+1)-1 AFTER NSF).
!+ NV    \ N     - NUMBER OF VARIABLES/EQUATIONS.
!+ NVA   \ R     - ORDERING OF THE ROWS OF M.
!+       \           SIZE = N.
!+ FVRA  \ U     - NONZERO ENTRIES IN_SS1 THE STRICT UPPER TRIANGULAR PORTION
!+       \           OF THE MATRIX U, STORED BY ROWS.
!+       \           SIZE = UMAX.
!+ NV    \ UMAX  - DECLARED DIMENSION OF U;  UMAX MUST BE LARGER THAN
!+       \           THE NUMBER OF NONZERO ENTRIES IN_SS1 THE STRICT UPPER
!+       \           TRIANGLE OF M PLUS FILLIN  (IU(N+1)-1 AFTER NSF).
!+ FRA   \ Z     - SOLUTION X.
!+       \           SIZE = N.
!
!history  E. RAZAFINDRAKOTO (LNH)
!+        20/11/06
!+        V5P7
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (LNHE)
!+        23/08/2012
!+        V6P2
!+   Size of ISP_SS1 doubled IN_SS1 non symmetric cases.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        18/04/2014
!+        V7P0
!+   Checking that 2*NSP is less than HUGE(1). Meshes with about
!+   2 millions of points will trigger memory allocations of numbers
!+   greater than the largest I4 integer.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DA             |-->| MATRIX DIAGONAL COEFFICIENTS
!| XA             |-->| OFF-DIAGONAL TERM OF MATRIX
!| GLOSEG         |-->| GLOBAL NUMBER OF SEGMENTS OF THE MATRIX
!| INFOGR         |-->| IF, YES INFORMATIONS ON LISTING
!| MAXSEG         |-->| MAXIMUM NUMBER OF SEGMENTS
!| NPOIN          |-->| NUMBER OF UNKNOWN
!| NSEGB          |-->| NUMBER OF SEGMENTS
!| RHS            |-->| SECOND MEMBER OF LINEAR EQUATION
!| TYPEXT         |---| = 'S' : SYMETRIC MATRIX
!| XINC           |<--| SOLUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_sd_solve_1 => sd_solve_1
      USE declarations_telemac, ONLY : in_ss1,ip_ss1,isegip_ss1,iw1_ss1,
     &                                 indtri_ss1,inx_ss1,ac_ss1,
     &                                 actri_ss1,isp_ss1,rsp_ss1,
     &                                 indtri_ss1,size_in,size_ip,
     &                                 size_isegip,size_iw1,size_indtri,
     &                                 size_inx,size_ac,size_actri,
     &                                 size_isp,size_rsp,size_ipx,
     &                                 ipx_ss1
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NPOIN,NSEGB,MAXSEG
      INTEGER, INTENT(IN)             :: GLOSEG(maxseg,2)
      LOGICAL, INTENT(IN)             :: INFOGR
      DOUBLE PRECISION, INTENT(IN)    :: XA(*),RHS(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: XINC(npoin),DA(npoin)
      CHARACTER(LEN=1), INTENT(IN)    :: TYPEXT
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER MEMFACTOR,IERR,NPBLK,NSEGBLK,NSP,ESP,INDIC,FLAG,I,IERRK
!
!
!
!-----------------------------------------------------------------------
!
!     CORRECTS DIAGONALS (TIDAL FLATS WITH MASKING)
!
      DO i=1,npoin
        IF(abs(da(i)).LT.1.d-15) da(i)=1.d0
      ENDDO
!
!-----------------------------------------------------------------------
!
      IF(infogr) THEN
        WRITE(lu,*) '                      DIRECT SYSTEM SOLVER'
      ENDIF
!
      npblk=npoin
      nsegblk=nsegb
!
!     1. MEMFACTOR: MEMORY FACTOR FOR SIZE ISP_SS1 AND RSP_SS1 IN_SS1 ODRV AND SDRV
!     =======================================================================
!
      IF(typext.EQ.'S') THEN
        memfactor = 5
      ELSE
        memfactor = 15
      ENDIF
      nsp=memfactor*(npblk+4*nsegblk)
      esp=memfactor*(npblk+4*nsegblk)
!
!     2. ALLOCATES ARRAYS (OR REALLOCATES IF TOO SMALL)
!     =======================================================================
!
      IF(size_in.EQ.0) THEN
        ALLOCATE(in_ss1(npblk+1))
        size_in=    npblk+1
      ELSEIF(       npblk+1.GT.size_in) THEN
        DEALLOCATE(in_ss1)
        ALLOCATE(in_ss1(npblk+1))
        size_in=    npblk+1
      ENDIF
!
      IF(size_ip.EQ.0) THEN
        ALLOCATE(ip_ss1(nsegblk*2))
        size_ip=    nsegblk*2
      ELSEIF(       nsegblk*2.GT.size_ip) THEN
        DEALLOCATE(ip_ss1)
        ALLOCATE(ip_ss1(nsegblk*2))
        size_ip=    nsegblk*2
      ENDIF
!
      IF(size_isegip.EQ.0) THEN
        ALLOCATE(isegip_ss1(nsegblk*2+1))
        size_isegip=    nsegblk*2+1
      ELSEIF(           nsegblk*2+1.GT.size_isegip) THEN
        DEALLOCATE(isegip_ss1)
        ALLOCATE(isegip_ss1(nsegblk*2+1))
        size_isegip=    nsegblk*2+1
      ENDIF
!
      IF(size_iw1.EQ.0) THEN
        ALLOCATE(iw1_ss1(npblk))
        size_iw1=    npblk
      ELSEIF(        npblk.GT.size_iw1) THEN
        DEALLOCATE(iw1_ss1)
        ALLOCATE(iw1_ss1(npblk))
        size_iw1=    npblk
      ENDIF
!
      IF(size_indtri.EQ.0) THEN
        ALLOCATE(indtri_ss1(npblk))
        size_indtri=    npblk
      ELSEIF(           npblk.GT.size_indtri) THEN
        DEALLOCATE(indtri_ss1)
        ALLOCATE(indtri_ss1(npblk))
        size_indtri=    npblk
      ENDIF
!
      IF(size_inx.EQ.0) THEN
        ALLOCATE(inx_ss1(npblk+1))
        size_inx=    npblk+1
      ELSEIF(        npblk+1.GT.size_inx) THEN
        DEALLOCATE(inx_ss1)
        ALLOCATE(inx_ss1(npblk+1))
        size_inx=    npblk+1
      ENDIF
!
      IF(size_ipx.EQ.0) THEN
        ALLOCATE(ipx_ss1(nsegblk*2+npblk+1))
        size_ipx=    nsegblk*2+npblk+1
      ELSEIF(        nsegblk*2+npblk+1.GT.size_ipx) THEN
        DEALLOCATE(ipx_ss1)
        ALLOCATE(ipx_ss1(nsegblk*2+npblk+1))
        size_ipx=    nsegblk*2+npblk+1
      ENDIF
!
      IF(size_ac.EQ.0) THEN
        ALLOCATE(ac_ss1(nsegblk*2+npblk+1))
        size_ac=    nsegblk*2+npblk+1
      ELSEIF(       nsegblk*2+npblk+1.GT.size_ac) THEN
        DEALLOCATE(ac_ss1)
        ALLOCATE(ac_ss1(nsegblk*2+npblk+1))
        size_ac=    nsegblk*2+npblk+1
      ENDIF
!
      IF(size_actri.EQ.0) THEN
        ALLOCATE(actri_ss1(npblk))
        size_actri=    npblk
      ELSEIF(          npblk.GT.size_actri) THEN
        DEALLOCATE(actri_ss1)
        ALLOCATE(actri_ss1(npblk))
        size_actri=    npblk
      ENDIF
!
      IF(typext.EQ.'S') THEN
!
        IF(size_isp.EQ.0) THEN
          ALLOCATE(isp_ss1(nsp))
        ELSEIF(        nsp.GT.size_isp) THEN
          DEALLOCATE(isp_ss1)
          ALLOCATE(isp_ss1(nsp))
        ENDIF
        size_isp=    nsp
!
      ELSE
!
        IF(2*nsp.LT.nsp) THEN
          WRITE(lu,*) 'SIZE OF LARGEST INTEGER  ',huge(1)
          WRITE(lu,*) 'TRESPASSED BY 2*NSP, NSP=',nsp
          CALL plante(1)
          stop
        ENDIF
!
        IF(size_isp.EQ.0) THEN
          ALLOCATE(isp_ss1(2*nsp))
        ELSEIF(      2*nsp.GT.size_isp) THEN
          DEALLOCATE(isp_ss1)
          ALLOCATE(isp_ss1(2*nsp))
        ENDIF
        size_isp=    2*nsp
!
      ENDIF
!
      IF(size_rsp.EQ.0) THEN
        ALLOCATE(rsp_ss1(esp))
        size_rsp=    esp
      ELSEIF(        esp.GT.size_rsp) THEN
        DEALLOCATE(rsp_ss1)
        ALLOCATE(rsp_ss1(esp))
        size_rsp=    esp
      ENDIF
!
!     3. BUILDS NONSYMMETRICAL COMPACT STORAGE (IN_SS1,IP_SS1)
!     WITHOUT THE DIAGONAL AND (INX_SS1,IPX_SS1) WITH THE DIAGONAL
!     =======================================================================
!
      CALL sd_strssd(npblk,nsegblk,gloseg(1,1),gloseg(1,2),
     &               in_ss1,ip_ss1,isegip_ss1,iw1_ss1)
!
      IF(typext.EQ.'S') THEN
!       XA IS THE OFF-DIAGONAL TERMS AND MAY COME DIRECTLY FROM TELEMAC
!       HENCE THE LOWER TRIANGULAR PART MAY NOT BE BUILT, WE GIVE TWICE XA
!       INSTEAD OF XA,XA(NSEGBLK+1)
        CALL sd_fabcad(npblk,nsegblk,in_ss1,ip_ss1,isegip_ss1,
     &                 indtri_ss1,iw1_ss1,inx_ss1,ipx_ss1,actri_ss1,xa,
     &                 xa,da,ac_ss1)
!                             ISTRI                !!
      ELSE
        CALL sd_fabcad(npblk,nsegblk,in_ss1,ip_ss1,isegip_ss1,
     &                 indtri_ss1,iw1_ss1,inx_ss1,ipx_ss1,actri_ss1,xa,
     &                 xa(nsegblk+1),da,ac_ss1)
      ENDIF
!
!     4. MINIMUM DEGREE PERMUTATION (YSMP PACKAGE)
!     =======================================================================
!
      indic=1
!
      CALL sd_odrv(npblk,inx_ss1,ipx_ss1,ac_ss1,in_ss1,iw1_ss1,nsp,
     &             isp_ss1,indic,flag)
!                                   PERM,INVP
!
      IF(flag.NE.0) THEN
        WRITE(lu,*) 'INCREASE THE MEMORY FACTOR (MEMFACTOR)',
     &              ' IN_SS1 THE ROUTINE SD_SOLVE_1'
        CALL plante(1)
        stop
      ENDIF
!
!---> SECOND MEMBER OF THE SYSTEM
!
!     5. LDLT DECOMPOSITION AND RESOLUTION (YSMP PACKAGE)
!
      IF(typext.EQ.'S') THEN
!                          PERM,INVP
        CALL sd_sdrv(npblk,in_ss1,iw1_ss1,inx_ss1,ipx_ss1,ac_ss1,rhs,
     &               xinc,nsp,isp_ss1,rsp_ss1,esp,indic,flag)
      ELSE
        CALL sd_cdrv(npblk,in_ss1,in_ss1,iw1_ss1,inx_ss1,ipx_ss1,
     &               ac_ss1,rhs,xinc,
     &               nsp,isp_ss1,rsp_ss1,esp,indic,flag)
      ENDIF
!
      IF(typext.EQ.'S') THEN
      IF(flag.NE.0) THEN
        ierr=flag-8*npblk
        IF(ierr.GT.0) THEN
          WRITE(lu,*) 'MATRIX WITH ZERO PIVOT AT ROW'
          WRITE(lu,*) ierr
        ELSE
          WRITE(lu,*) 'ADD 1 TO THE MEMORY FACTOR (MEMFACTOR)',
     &                ' IN_SS1 SUBROUTINE SD_SOLVE_1'
        ENDIF
        CALL plante(1)
        stop
      ENDIF
      ELSE
!
!---> COMMENTS THE ERROR: FLAG_SD_NDRV:
!
! FLAG  - ERROR FLAG;  VALUES AND THEIR MEANINGS ARE --
!             0     NO ERRORS DETECTED
!             N+K   NULL ROW IN_SS1 A  --  ROW = K
!            2N+K   DUPLICATE ENTRY IN_SS1 A  --  ROW = K
!            3N+K   INSUFFICIENT STORAGE IN_SS1 NSF  --  ROW = K
!            4N+1   INSUFFICIENT STORAGE IN_SS1 NNF
!            5N+K   NULL PIVOT  --  ROW = K
!            6N+K   INSUFFICIENT STORAGE IN_SS1 NSF  --  ROW = K
!            7N+1   INSUFFICIENT STORAGE IN_SS1 NNF
!            8N+K   ZERO PIVOT  --  ROW = K
!           10N+1   INSUFFICIENT STORAGE IN_SS1 NDRV
!           11N+1   ILLEGAL PATH SPECIFICATION (INDIC)
      IF(flag.NE.0) THEN
        ierr=int(flag/npblk)
        IF(ierr.EQ.3.OR.ierr.EQ.5.OR.ierr.EQ.8) THEN
          ierrk=flag-ierr*npblk
          WRITE(lu,*) 'MATRIX WITH ZERO PIVOT AT ROW'
          WRITE(lu,*) ierrk
        ELSE
          WRITE(lu,*) 'INCREASE THE MEMORY FACTOR (MEMFACTOR)',
     &                ' IN_SS1 SUBROUTINE SD_SOLVE_1'
        ENDIF
        CALL plante(1)
        stop
      ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END