      SUBROUTINE SLVSPI(A,B,X,MDIAG,NEQ,KEY,IERR)
C     CALL SLVSPI (A,MAVAIL,B,X,MDIAG,NEQ,KEY,IERR)
C     USAGE : IN-CORE SYMMETRIC POSITIVE DEFINITE SYSTEM.
C     =====
C     PURPOSE
C     -------
C     TO SOLVE A LINEAR SYSTEM OF EQUATIONS [A] [X] = [B] , WHERE A
C     IS TOTAL STIFFNESS MATRIX
C     INPUT ARGUMENTS
C     ---------------
C     A       = VECTOR OF LENGTH MAVAIL CONTAINING TOTAL STIFFNESS
C               MATRIX.(AS OUTPUTTED BY ADDASI)
C     MAVAIL  = AVAILABLE MEMORY ; MUST BE AS OUTPUTTED BY FSTSPI.
C     B       = LOAD VECTOR OF LENGTH NEQ (SEE NOTE 1).
C     MDIAG   = VECTOR OF LENGTH NEQ WHICH CONTAINS THE POSITION
C               OF THE DIAGONAL ELEMENT IN THE TOTAL STIFFNESS
C               MATRIX.(OUTPUTTED BY FSTSPI)
C     NEQ     = NUMBER OF EQUATIONS.
C     KEY     = OPERATION CONTROL KEY :
C               KEY = 1  PERFORM LDLT FACTORISATION OF A AND SOLVE
C                        AX = B.
C               KEY = 2  PERFORM LDLT FACTORISATION OF A ONLY.
C               KEY = 3  REDUCTION OF LOAD VECTOR B AND BACK
C                        SUBSTITUTION TO OBTAIN SOLUTION X (SEE NOTE 2).
C     OUTPUT ARGUMENTS
C     ----------------
C     X       = SOLUTION VECTOR OF LENGTH NEQ (SEE NOTE 1).
C     IERR    = ERROR FLAG (SEE NOTE 5).
C     NOTES
C     -----
C     1.  THE VECTORS B AND X CAN OCCUPY THE SAME STORAGE SPACE I.E. ON
C         OUTPUT B WILL CONTAIN THE SOLUTION.
C     2.  IF KEY = 3 , SLVSPI MUST HAVE PREVIOUSLY BEEN CALLED WITH
C         KEY = 1 OR 2. SLVSPI CAN BE CALLED REPEATEDLY WITH KEY = 3,
C         PROVIDING A AND MDIAG ARE NOT CHANGED.
C     3.  SLVSPI USES THE LIBRARY-SUPPLIED FUNCTION VCSLPR(X,Y,N) WHICH
C         CALCULATES THE INNER (DOT,SCALAR) PRODUCT OF THE VECTORS X,Y
C         OF LENGTH N.
C     4.  THE SUBROUTINE SLVSPI CAN BE USED INDEPENDENTLY OF THE PACKAGE
C         PRESENTED HERE. THE SOLE REQUIREMENT OF SLVSPI IS THE VECTORS
C         A,B AND MDIAG WRITTEN IN THE PRESCRIBED FORMAT.
C     5.  IERR = 0  EXECUTION SUCCESSFULL
C              = N  A ZERO PIVOT WAS ENCOUNTERED ON ROW N DURING LDLT
C                   FACTORISATION.FACTORISATION ABORTED.
C              = -1 A NEGATIVE PIVOT WAS ENCOUNTERED FACTORISATION
C                   CONTINUED.
      DIMENSION A(*),B(NEQ),X(NEQ),MDIAG(NEQ)
      IERR = 0
      IF (KEY.EQ.3)GO TO 60
C     LOOP ON ROWS OF MATRIX.
      MDI = 0
      DO 50  I = 1,NEQ
         MDIM  = MDI
         MDI   = MDIAG(I)
         LOCI0 = MDI - I
         LOWI  = MDIM + 1 - LOCI0
         IM    = I - 1
         IFJ   = LOWI
         ILJ   = IM
         IF (IFJ.GT.ILJ) GO TO 20
         JJ    = IFJ - 1
         MDJ   = 0
         IF (IFJ.GT.1) MDJ = MDIAG(JJ)
C        LOOP ON COLUMNS OF MATRIX.
         DO 10  J = IFJ,ILJ
            MDJM  =  MDJ
            MDJ   =  MDIAG(J)
            LOCJ0 = MDJ - J
            LOWJ  = MDJM + 1 - LOCJ0
C           COMPUTE L(I,J)*D(J,J) IN A(LOCIJ)
            K     = MAX0(LOWI,LOWJ)
            LOCIJ = LOCI0 + J
            LOCIK = LOCI0 + K
            LOCJK = LOCJ0 + K
            A(LOCIJ)= A(LOCIJ)-VCSLPR(A(LOCIK),A(LOCJK),J-K)
   10    CONTINUE
C           COMPUTE D(I,I) AND STORE ITS RECIPORICAL IN A(MDI)
C           COMPUTE L(I,J)
   20    AA    = A(MDI)
         IF (LOWI.GT.IM)GO TO 40
         DO 30 K = LOWI,IM
            MDK  = MDIAG(K)
            LOCIK= LOCI0 + K
            BB   = A(LOCIK)*A(MDK)
            AA   = AA - BB*A(LOCIK)
         A(LOCIK)= BB
   30    CONTINUE
   40    IF (AA.EQ.0.) GO TO 120
         IF (AA.LT.0.) IERR = -1
         A(MDI)  = 1./AA
         IF (KEY.EQ.2) GO TO 50
C        REDUCE B IF KEY = 1
         LOCIK = LOCI0 + LOWI
         X(I)  = B(I) - VCSLPR(A(LOCIK),X(LOWI),I-LOWI)
   50 CONTINUE
      IF (KEY.EQ.2) RETURN
      GO TO 80
C     REDUCE B IF KEY = 3
   60 MDI   = 0
      DO 70 I = 1,NEQ
         MDIM = MDI
         MDI  = MDIAG(I)
         LOCI0= MDI - I
         LOWI = MDIM + 1 - LOCI0
         LOCIK= MDIM + 1
         X(I) = B(I) - VCSLPR(A(LOCIK),X(LOWI),I-LOWI)
   70 CONTINUE
   80 DO 90 I = 1,NEQ
         MDI  = MDIAG(I)
         X(I) = X(I)*A(MDI)
   90 CONTINUE
C     BACK SUBSTITUTION
      MDIM = MDIAG(NEQ)
      IM   = NEQ
      DO 110 LL = 1,NEQ
         I      = IM
         IM     = IM - 1
         MDI    = MDIM
         MDIM   = 0
         IF (I.NE.1) MDIM = MDIAG(IM)
         LOCI0  = MDI - I
         LOWI   = MDIM + 1 - LOCI0
         IF (LOWI.GT.IM) GO TO 110
         LOCIJ  = MDIM + 1
         DO 100 J = LOWI,IM
            X(J)  = X(J) - A(LOCIJ)*X(I)
            LOCIJ = LOCIJ + 1
  100    CONTINUE
  110 CONTINUE
CCC      IF (IERR.EQ.-1) CALL ERFAC(4HSLVS,IERR)
      RETURN
C     ZERO PIVOT
  120 IERR = I
C      MERF=100*IERR
C      CALL ERFAC(4HSLVS,MERF)
C      RETURN
      END
C===================================================================
      real function VCSLPR (X,Y,N)
      DIMENSION X(*),Y(*)
      vcslpr = 0.
      DO 100 I=1,N 
        VCSLPR=VCSLPR + X(I)*Y(I)
  100 CONTINUE
      END
C===================================================================
      SUBROUTINE ERFAC(ERSUB,MERR)
C     CALL ERFAC(ERSUB,MERR)
C     PURPOSE
C     -------
C     TO OUTPUT A PRINTED MESSAGE FOR ANY ERROR ENCOUNTED.
C     INPUT ARGUMENTS
C     ---------------
C     ERSUB  = NAME OF THE ROUTINE IN WHICH ERROR WAS FOUND.
C     MERR   = ERROR FLAG (SEE NOTE 1)
C     OUTPUT ARGUMENTS
C     ----------------
C     PRINTS ERRORS AND SUBROUTINE NAME.
C     NOTES
C     -----
C     1. MERR .LT. 100 NON-FATAL ERROR.
C             .GE. 100 FATAL ERROR.
      IF(MERR.GE.100.OR.MERR.LE.-100) GO TO 100
      PRINT 1000,ERSUB,MERR
      RETURN
 100  MERR=MERR/100
      PRINT 1010, ERSUB,MERR
      STOP
C     FORMATS
 1000 FORMAT(/
     &' *** WARNING *** NON-FATAL ERROR IN SUBROUTINE ',A6/
     &' ERROR FLAG ',I5/
     &' PROGRAM CONTINUING')
 1010 FORMAT(/
     &' *** FATAL ERROR *** IN SUBROUTINE ',A6/
     &' ERROR FLAG ',I5)
      END
