* SUBROUTINE MXDPGB                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* OBTAINED BY THE SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
*         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
*         EQUATIONS.
*  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
*         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXDPGB(N,A,X,JOB)
      INTEGER          JOB,N
      DOUBLE PRECISION A(*),X(*)
      INTEGER          I,II,IJ,J
      IF (JOB.GE.0) THEN
*
*     PHASE 1 : X:=L**(-1)*X
*
          IJ = 0
          DO 20 I = 1,N
              DO 10 J = 1,I - 1
                  IJ = IJ + 1
                  X(I) = X(I) - A(IJ)*X(J)
   10         CONTINUE
              IJ = IJ + 1
   20     CONTINUE
      END IF
      IF (JOB.EQ.0) THEN
*
*     PHASE 2 : X:=D**(-1)*X
*
          II = 0
          DO 30 I = 1,N
              II = II + I
              X(I) = X(I)/A(II)
   30     CONTINUE
      END IF
      IF (JOB.LE.0) THEN
*
*     PHASE 3 : X:=TRANS(L)**(-1)*X
*
          II = N* (N-1)/2
          DO 50 I = N - 1,1,-1
              IJ = II
              DO 40 J = I + 1,N
                  IJ = IJ + J - 1
                  X(I) = X(I) - A(IJ)*X(J)
   40         CONTINUE
              II = II - I
   50     CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXDPGF                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* FACTORIZATION A+E=L*D*TRANS(L) OF A DENSE SYMMETRIC POSITIVE DEFINITE
* MATRIX A+E WHERE D AND E ARE DIAGONAL POSITIVE DEFINITE MATRICES AND
* L IS A LOWER TRIANGULAR MATRIX. IF A IS SUFFICIENTLY POSITIVE
* DEFINITE THEN E=0.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  ON INPUT A GIVEN DENSE SYMMETRIC (USUALLY POSITIVE
*         DEFINITE) MATRIX A STORED IN THE PACKED FORM. ON OUTPUT THE
*         COMPUTED FACTORIZATION A+E=L*D*TRANS(L).
*  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
*         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
*         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
*         INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
*         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
*         PROCESS.
*  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
*         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
*         FACTORIZATION PROCESS (IF INF>0).
*  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
*
* METHOD :
* P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
* LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
* PP. 311-350.
*
      SUBROUTINE MXDPGF(N,A,INF,ALF,TAU)
      DOUBLE PRECISION ALF,TAU
      INTEGER          INF,N
      DOUBLE PRECISION A(*)
      DOUBLE PRECISION BET,DEL,GAM,RHO,SIG,TOL
      INTEGER          I,IJ,IK,J,K,KJ,KK,L
      L = 0
      INF = 0
      TOL = ALF
*
*     ESTIMATION OF THE MATRIX NORM
*
      ALF = 0.0D0
      BET = 0.0D0
      GAM = 0.0D0
      TAU = 0.0D0
      KK = 0
      DO 20 K = 1,N
          KK = KK + K
          BET = MAX(BET,ABS(A(KK)))
          KJ = KK
          DO 10 J = K + 1,N
              KJ = KJ + J - 1
              GAM = MAX(GAM,ABS(A(KJ)))
   10     CONTINUE
   20 CONTINUE
      BET = MAX(TOL,BET,GAM/N)
*      DEL = TOL*BET
      DEL = TOL*MAX(BET,1.0D0)
      KK = 0
      DO 60 K = 1,N
          KK = KK + K
*
*     DETERMINATION OF A DIAGONAL CORRECTION
*
          SIG = A(KK)
          IF (ALF.GT.SIG) THEN
              ALF = SIG
              L = K
          END IF
          GAM = 0.0D0
          KJ = KK
          DO 30 J = K + 1,N
              KJ = KJ + J - 1
              GAM = MAX(GAM,ABS(A(KJ)))
   30     CONTINUE
          GAM = GAM*GAM
          RHO = MAX(ABS(SIG),GAM/BET,DEL)
          IF (TAU.LT.RHO-SIG) THEN
              TAU = RHO - SIG
              INF = -1
          END IF
*
*     GAUSSIAN ELIMINATION
*
          A(KK) = RHO
          KJ = KK
          DO 50 J = K + 1,N
              KJ = KJ + J - 1
              GAM = A(KJ)
              A(KJ) = GAM/RHO
              IK = KK
              IJ = KJ
              DO 40 I = K + 1,J
                  IK = IK + I - 1
                  IJ = IJ + 1
                  A(IJ) = A(IJ) - A(IK)*GAM
   40         CONTINUE
   50     CONTINUE
   60 CONTINUE
      IF (L.GT.0 .AND. ABS(ALF).GT.DEL) INF = L
      RETURN
      END
* SUBROUTINE MXSCMD               ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A DENSE RECTANGULAR MATRIX A BY A VECTOR X AND
* ADDITIOON OF THE SCALED VECTOR ALF*Y.
*
* PARAMETERS :
*  II  N  NUMBER OF ROWS OF THE MATRIX A.
*  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
*  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
*  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
*  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
*  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
*  RI  X(NA)  INPUT VECTOR.
*  RI  ALF  SCALING FACTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
*
* SUBPROGRAMS USED :
*  S   MXVSCL  SCALING OF A VECTOR.
*
      SUBROUTINE MXSCMD(N,NA,A,IA,JA,X,ALF,Y,Z)
      INTEGER N,NA,IA(*),JA(*)
      DOUBLE PRECISION A(*),X(*),ALF,Y(*),Z(*)
      INTEGER I,J,K,L,JP
      CALL MXVSCL(N,ALF,Y,Z)
      DO 2 I=1,NA
      K=IA(I)
      L=IA(I+1)-K
      DO 1 J=1,L
      JP=JA(K)
      IF (JP.GT.0) Z(JP)=Z(JP)+A(K)*X(I)
      K=K+1
    1 CONTINUE
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXSCMM               ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A DENSE RECTANGULAR MATRIX A BY A VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF ROWS OF THE MATRIX A.
*  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
*  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
*  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
*  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
*  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
*  RI  X(NA)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR EQUAL TO A*X.
*
* SUBPROGRAMS USED :
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE MXSCMM(N,NA,A,IA,JA,X,Y)
      INTEGER N,NA,IA(*),JA(*)
      DOUBLE PRECISION A(*),X(*),Y(*)
      INTEGER I,J,K,L,JP
      CALL MXVSET(N,0.0D 0,Y)
      DO 2 I=1,NA
      K=IA(I)
      L=IA(I+1)-K
      DO 1 J=1,L
      JP=JA(K)
      IF (JP.GT.0) Y(JP)=Y(JP)+A(K)*X(I)
      K=K+1
    1 CONTINUE
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXSPCA                 ALL SYSTEMS                93/12/01
C PORTABILITY : ALL SYSTEMS
C 93/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* REWRITE SYMMETRIC MATRIX INTO THE PERMUTED FACTORIZED COMPACT SCHEME.
* MOIDIFIED VERSION FOR THE USE WITH MXSPCJ.
*
* PARAMETERS:
*  II  N  SIZE OF THE SYSTEM SOLVED.
*  II  NB NUMBER OF NONZEROS IN THE UPPER TRIANGLE OF THE ORIGINAL
*         MATRIX.
*  II  ML SIZE OF THE COMPACT FACTOR.
*  RU  A(MMAX) NUMERICAL VALUES OF THE SPARSE HESSIAN APPROXIMATION
*              STORED AT THE POSITIONS 1, ...,NB.
*  IU  JA(MMAX) INDICES OF THE NONZERO ELEMENTS OF THE HESSIAN MATRIX IN
*             THE PACKED ROW FORM AT THE FIRST NB POSITIONS.
*             NEW POSITIONS
*             IN THE PERMUTED FACTOR STORED IN A(NB+1), ..., A(2*NB),
*             INDICES OF COMPACT SCHEME IN A(2*NB+1), ..., A(2*NB+ML).
*  II  PSL(N+1)  POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR
*             FACTOR OF THE HESSIAN APPROXIMATION.
*  RI  T  CORRECTION FACTOR THAT IS ADDED TO THE DIAGONAL.
*
*
      SUBROUTINE MXSPCA(N,NB,ML,A,IA,JA,T)
      INTEGER N,NB,ML,IA(*),JA(*)
      DOUBLE PRECISION A(*),T
      INTEGER I,J
      DO 100 I=1,N
        J=ABS(JA(IA(I)+NB+ML))
        A(NB+J)=A(NB+J)+T
 100  CONTINUE
      RETURN
      END
* SUBROUTINE MXSPCB                ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* STORED IN THE COMPACT SCHEME. THIS FACTORIZATION CAN BE OBTAINED
* USING THE SUBROUTINE MXSSFN.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(MMAX)  FACTORS L,D OF THE FACTORIZATION A+E=L*D*TRANS(L)
*                STORED USING THE COMPACT SCHEME OF STORING.
*  II  PSL(N+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX
*  II  SL(MMAX)  ARRAY OF COLUMN INDICES OF THE FACTORS L AND D
*         STORED USING THE COMPACT SCHEME.
*  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
*         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
*         EQUATIONS.
*  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
*         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXSPCB(N,A,PSL,SL,X,JOB)
      INTEGER N
      DOUBLE PRECISION  A(*),X(*)
      INTEGER PSL(*),SL(*),JOB
      INTEGER I,J,IS
C
C     FIRST PHASE
C
      IF (JOB.GE.0) THEN
                    DO 300 I=1,N
                    IS=SL(I)+N+1
                    DO 200 J=PSL(I)+I,PSL(I+1)+I-1
                    X(SL(IS))=X(SL(IS)) - A(J)*X(I)
                    IS=IS+1
200                 CONTINUE
300                 CONTINUE
      END IF
C
C     SECOND PHASE
C
      IF (JOB.EQ.0) THEN
                    DO 400 I=1,N
                    X(I) = X(I)/A(PSL(I)+I-1)
400   CONTINUE
      END IF
C
C     THIRD PHASE
C
      IF (JOB.LE.0) THEN
                    DO 600 I=N,1,-1
                    IS=SL(I)+N+1
                    DO 500 J=PSL(I)+I,PSL(I+1)+I-1
                    X(I)=X(I)-A(J)*X(SL(IS))
                    IS=IS+1
500                 CONTINUE
600                 CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXSPCC                ALL SYSTEMS               92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* SPARSE MATRIX REORDER, SYMBOLIC FACTORIZATION, DATA STRUCTURES
* TRANSFORMATION - INITIATION OF THE DIRECT SPARSE SOLVER.
* MODIFIED VERSION WITH CHANGED DATA STRUCTURES.
*
*  PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NJA  NUMBER OF NONZERO ELEMENTS OF THE MATRIX.
*  IO  ML  SIZE OF THE COMPACT STRUCTURE OF THE TRIANGULAR FACTOR
*         OF THE HESSIAN APPROXIMATION.
*  II  MMAX  SIZE OF THE ARRAYS JA,A.
*  RO  A(MMAX)   NUMERICAL VALUES OF THE SPARSE HESSIAN APPROXIMATION
*         STORED AT THE POSITIONS 1, ...,NJA. LOWER TRIANGULAR
*         FACTOR OF THE HESSIAN APPROXIMATION STORED AT THE
*         POSITIONS NJA+1, ..., NJA+MB.
*  II  IA(N) POINTERS OF THE DIAGONAL ELEMENTS OF THE HESSIAN MATRIX.
*  II  JA(MMAX)  INDICES OF THE NONZERO ELEMENTS OF THE HESSIAN MATRIX IN
*         THE PACKED ROW FORM AT THE FIRST NJA POSITIONS. COMPACT
*         STRUCTURE OF INDICES OF ITS TRIANGULAR FACTOR IS ROWWISE
*         STORED.
*  II  PSL(N+1)  POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR
*         FACTOR OF THE HESSIAN APPROXIMATION.
*  IO  PERM(N) PERMUTATION VECTOR.
*  IO  INVP(N) INVERSE PERMUTATION VECTOR.
*  IA  WN11(N) AUXILIARY VECTOR.
*  IA  WN12(N) AUXILIARY VECTOR.
*  IA  WN13(N) AUXILIARY VECTOR.
*  IA  WN14(N) AUXILIARY VECTOR.
*  IO  ITERM  TERMINATION INDICATOR. TERMINATION IF ITERM .NE. 0.
*
* COMMON DATA :
*
* SUBPROGRAMS USED :
*  S   MXSTG1  WIDTHENING OF THE PACKED FORM OF THE SPARSE MATRIX.
*  S   MXSSMN  SPARSE MATRIX REORDERING.
*  S   MXVSIP  INVERSE PERMUTATION COMPUTING.
*  S   MXSPCI  SYMBOLIC FACTORIZATION.
*  S   MXSTL1  PACKING OF THE WIDTHENED FORM OF THE SPARSE MATRIX.
*
      SUBROUTINE MXSPCC(N,NJA,ML,MMAX,A,IA,JA,PSL,PERM,INVP,WN11,WN12,
     *                  WN13,WN14,ITERM)
      INTEGER N,NJA,MMAX,ML,ITERM
      INTEGER PERM(*),INVP(*),WN11(*),WN12(*),WN13(*),WN14(*)
      INTEGER PSL(*),IA(*),JA(*)
      INTEGER JSTRT,JSTOP,I,J,K,L,NJASAVE
      INTEGER LL,LL1,NJABIG,KSTRT
      DOUBLE PRECISION A(*)
      IF (ML.GT.0) RETURN
      IF(2*NJA.GE.MMAX) THEN
        ITERM=-41
        GO TO 1000
      END IF
C
C WIDTHENING OF THE PACKED FORM
C
      NJASAVE=NJA
      CALL MXSTG1(N,NJA,IA,JA,WN12,WN11)
      NJABIG=NJA
C
C REORDERING OF THE SPARSE MATRIX
C
      CALL MXSSMN(N,IA,JA,PERM,WN11,WN12,WN13)
C
C FIND THE INVERSE PERMUTATION VECTOR INVP
C
      CALL MXVSIP(N,PERM,INVP)
C
C SHRINK THE STRUCTURE
C
      CALL MXSTL1(N,NJA,IA,JA,WN11)
C
      DO 40 I=1,N
        WN11(I)=0
        WN12(I)=0
 40   CONTINUE
C
C WN11 CONTAINS BEGINNINGS OF THE FACTOR ROWS
C
      DO 50 I=1,N
        K=PERM(I)
        JSTRT=IA(K)
        JSTOP=IA(K+1)-1
        DO 60 J=JSTRT,JSTOP
          L=JA(J)
          L=INVP(L)
          IF(L.GE.I) THEN
            WN12(I)=WN12(I)+1
          ELSE
            WN12(L)=WN12(L)+1
          END IF
 60     CONTINUE
 50   CONTINUE
      WN11(1)=1
      DO 69 I=1,N-1
       WN11(I+1)=WN11(I)+WN12(I)
 69   CONTINUE
C
C CREATE UPPER TRIANGULAR STRUCTURE NECESSARY FOR THE TRANSFER
C
      DO 300 I=1,N
        K=PERM(I)
        JSTRT=IA(K)
        JSTOP=IA(K+1)-1
        DO 200 J=JSTRT,JSTOP
          L=JA(J)
          L=INVP(L)
          IF(L.GE.I) THEN
            LL1=WN11(I)
            WN11(I)=LL1+1
            JA(NJABIG+LL1)=L
            A(J)=LL1
            A(NJA+LL1)=J
          ELSE
            LL1=WN11(L)
            WN11(L)=LL1+1
            JA(NJABIG+LL1)=I
            A(J)=LL1
            A(NJA+LL1)=J
          END IF
 200    CONTINUE
 300  CONTINUE
C
C SORT INDICES IN THE PERMUTED UPPER TRIANGLE
C
      DO 599 I=1,N
        WN11(I)=0
 599  CONTINUE
      WN11(1)=1
      WN14(1)=1
      DO 67 I=2,N+1
        WN11(I)=WN11(I-1)+WN12(I-1)
        WN14(I)=WN11(I)
   67 CONTINUE
 80   CONTINUE
      DO 602 I=1,N
        WN12(I)=0
 602  CONTINUE
       JSTOP=WN11(N+1)
       DO 499 I=N,1,-1
         JSTRT=WN11(I)
         CALL MXVSR5(JSTOP-JSTRT,JSTRT-1,JA(NJABIG+JSTRT),
     *               A,A(NJASAVE+JSTRT))
         JSTOP=JSTRT
 499   CONTINUE
C
C WIDTHENING OF THE PERMUTED PACKED FORM.
C
      NJASAVE=NJA
C
      CALL MXSTG1(N,NJA,IA,JA,WN12,WN11)
C
      NJABIG=NJA
C
C
C SYMBOLIC FACTORIZATION.
C
      CALL MXSPCI(N,ML,MMAX-2*NJA,IA,JA,PSL,A(2*NJASAVE+1),PERM,
     *            INVP,WN11,WN12,WN13,ITERM)
      IF(ITERM.NE.0) THEN
        ITERM=-42
        GO TO 1000
      END IF
C
C
C RETRIEVE PARAMETERS
C
      CALL MXSTL1(N,NJA,IA,JA,WN11)
C
C SHIFT PERMUTED UPPER TRIANGLE.
C
      DO 73 I=1,NJASAVE
        JA(NJA+I)=JA(NJABIG+I)
  73  CONTINUE
C
C SHIFT STRUCTURE SL.
C
      IF(2*NJASAVE+ML.GE.MMAX) THEN
        ITERM=-41
        GO TO 1000
      END IF
      DO 70 I=1,ML
        JA(2*NJASAVE+I)=A(2*NJASAVE+I)
  70  CONTINUE
C
C SET POINTERS
C
      DO 600 I=1,N
        WN12(I)=0
 600  CONTINUE
       LL1=PSL(N)+N-1
       JSTOP=WN14(N+1)
       DO 500 I=N,1,-1
         JSTRT=WN14(I)
         DO 700 J=JSTRT,JSTOP-1
           K=JA(NJASAVE+J)
           WN12(K)=J
           LL=A(NJASAVE+J)
           WN13(K)=LL
 700     CONTINUE
         JSTOP=JSTRT
         KSTRT=JA(2*NJASAVE+I)+N+1+2*NJASAVE
         DO 800 J=KSTRT+PSL(I+1)-PSL(I)-1,KSTRT,-1
           L=JA(J)
           IF(WN12(L).NE.0) THEN
             LL=WN13(L)
             A(LL)=LL1
             WN12(L)=0
           END IF
           LL1=LL1-1
 800     CONTINUE
         K=WN12(I)
         WN12(I)=0
         LL=WN13(I)
         A(LL)=LL1
         LL1=LL1-1
 500   CONTINUE
C
      DO 76 I=1,ML
        JA(NJASAVE+I)=JA(2*NJASAVE+I)
  76  CONTINUE
C
      DO 72 I=1,NJASAVE
        JA(ML+NJASAVE+I)=A(I)
  72  CONTINUE
C
 1000 CONTINUE
      RETURN
      END
* SUBROUTINE MXSPCF                ALL SYSTEMS                90/12/01
C PORTABILITY : ALL SYSTEMS
C 90/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* NUMERICAL  FACTORIZATION A+E=L*D*TRANS(L) OF A SPARSE
* SYMMETRIC POSITIVE DEFINITE MATRIX A+E WHERE D AND E ARE DIAGONAL
* POSITIVE DEFINITE MATRICES AND L IS A LOWER TRIANGULAR MATRIX. IF
* A IS SUFFICIENTLY POSITIVE DEFINITE THEN E=0. THE STRUCTURE ON
* INPUT WAS OBTAINED BY THE SYMBOLIC FACTORIZATION AND IT MAKES
* USE OF THE COMPACT SCHEME OF STORING THE SPARSE MATRIX IN THE
* POINTER ARRAY PSL ,INDEX ARRAY SL. NUMERICAL VALUES OF THE FACTOR
* CAN BE FOUND IN THE ARRAY A.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(MMAX) ON INPUT NUMERICAL VALUES OF THE LOWER HALF OF THE
*      MATRIX THAT IS BEEING FACTORIZED(INCLUDING THE DIAGONAL
*      ELEMENTS. ON OUTPUT IT CONTAINS FACTORS L AND D AS IF THEY
*      FORM THE LOWER HALF OF THE MATRIX.STRUCTURE INFORMATION
*      IS SAVED IN THE COMPACT SCHEME IN THE PAIR OF VECTORS PSL
*      AND SL.
*  II  PSL(NF+1) POINTER VECTOR OF THE FACTOR
*  II  SL(MMAX) STRUCTURE OF THE FACTOR IN THE COMPACT FORM
*  IA  WN11(NF+1) AUXILIARY VECTOR.
*  IA  WN12(NF+1) AUXILIARY VECTOR.
*  RA  RN01(NF+1) AUXILIARY VECTOR.
*  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
*         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
*         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
*         INF>0 THEN THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
*         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
*         PROCESS.
*  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
*         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
*         FACTORIZATION PROCESS (IF INF>0).
*  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
*
* METHOD :
* S.C.EISENSTAT,M.C.GURSKY,M.H.SCHULTZ,A.H.SHERMAN:YALE SPARSE MATRIX
* PACKAGE I. THE SYMMETRIC CODES,YALE UNIV. RES. REPT.
* NO.112,1977.
*
      SUBROUTINE MXSPCF(N,A,PSL,SL,WN11,WN12,RN01,INF,ALF,TAU)
      INTEGER N,PSL(*),SL(*),WN11(*),WN12(*),INF
      DOUBLE PRECISION A(*),RN01(*),ALF
      DOUBLE PRECISION BET,GAM,DEL,RHO,SIG,TOL,TADD,TBDD,TAU
      INTEGER  I, J, K, L, II
      INTEGER ISTRT,ISTOP,NEWK,KPB,ISUB
      L = 0
      INF = 0
      TOL = ALF
      ALF = 0.0D 0
      BET = 0.0D 0
      GAM = 0.0D 0
      TAU = 0.0D 0
      DO 60 I=1,N
      BET=MAX(BET,ABS(A(PSL(I)+I-1)))
      DO 50 J=PSL(I)+I,PSL(I+1)+I-1
      GAM = MAX( GAM, ABS(A(J)) )
50    CONTINUE
60    CONTINUE
      BET = MAX(TOL,BET,GAM/N)
      DEL = TOL*BET
      DO 100 I=1,N
      WN11(I)=0
      RN01(I)=0.0D 0
100   CONTINUE
      DO 600 J=1,N
C
C     DETERMINATION OF A DIAGONAL CORRECTION
C
      SIG=A(PSL(J)+J-1)
      RHO=0.0D 0
      NEWK=WN11(J)
200   K=NEWK
      IF(K.EQ.0) GO TO 400
      NEWK=WN11(K)
      KPB=WN12(K)
      TADD=A(KPB+K)
      TBDD=TADD*A(PSL(K)+K-1)
      RHO=RHO+TADD*TBDD
      ISTRT=KPB+1
      ISTOP=PSL(K+1)-1
      IF(ISTOP.LT.ISTRT) GO TO 200
C
      WN12(K)=ISTRT
      I=SL(K)+(KPB-PSL(K))+1
      ISUB=SL(N+1+I)
      WN11(K)=WN11(ISUB)
      WN11(ISUB)=K
      DO 300 II=ISTRT,ISTOP
      ISUB=SL(N+1+I)
      RN01(ISUB)=RN01(ISUB)+A(II+K)*TBDD
      I=I+1
300   CONTINUE
      GO TO 200
400   SIG=A(PSL(J)+J-1)-RHO
      IF(ALF.GT.SIG) THEN
        ALF=SIG
        L=J
      END IF
      GAM=0.0D 0
      ISTRT=PSL(J)
      ISTOP=PSL(J+1)-1
      IF(ISTOP.LT.ISTRT) GO TO 370
      WN12(J)=ISTRT
      I=SL(J)
      ISUB=SL(N+1+I)
      WN11(J)=WN11(ISUB)
      WN11(ISUB)=J
      DO 500 II=ISTRT,ISTOP
      ISUB=SL(N+1+I)
      A(II+J)=(A(II+J)-RN01(ISUB))
      RN01(ISUB)=0.0D 0
      I=I+1
500   CONTINUE
      DO 350 K=PSL(J)+J,PSL(J+1)+J-1
      GAM=MAX(GAM,ABS(A(K)))
350   CONTINUE
      GAM=GAM*GAM
370   RHO=MAX(ABS(SIG),GAM/BET,DEL)
      IF(TAU.LT.RHO-SIG) THEN
                         TAU=RHO-SIG
                         INF=-1
      END IF
C
C     GAUSSIAN ELIMINATION
C
      A(PSL(J)+J-1)=RHO
      DO 550 II=ISTRT,ISTOP
      A(II+J)=A(II+J)/RHO
550   CONTINUE
600   CONTINUE
      IF (L.NE.0.AND.ABS(ALF).GT.DEL) INF=L
      RETURN
      END
* SUBROUTINE MXSPCI                ALL SYSTEMS                89/12/01
C PORTABILITY : ALL SYSTEMS
C 89/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* SYMBOLIC FACTORIZATION OF A SPARSE SYMMETRIC MATRIX GIVEN IN THE
* NORMAL SCHEME PA,SA. ON OUTPUT WE HAVE POINTER VECTOR OF THE FACTOR
* PSL AND VECTOR OF COLUMN INDICES SL. ML IS THE NUMBER OF THE INDICES
* USED FOR THE VECTOR SL, WHERE WE HAVE FACTOR IN THE COMPACT FORM.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  IO  ML NUMBER OF THE NONZERO ELEMENTS IN THE FACTOR'S COMPACT SCHEME
*  II  MMAX  LENGTH OF THE ARRAY SL. IN THE CASE OF THE
*            INSUFFICIENT SPACE IT IS TO BE INCREASED.
*  II  PA(N+1) POINTER VECTOR OF THE INPUT MATRIX
*  II  SA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX
*  IO  PSL(N+1) POINTER VECTOR OF THE FACTOR
*  RO  SL(MMAX) COMPACT SCHEME OF THE INDICES OF THE FACTOR
*  II  PERM(N) PERMUTATION VECTOR
*  II  INVP(N) INVERSE PERMUTATION VECTOR
*  IA  WN11(N+1) WORK VECTOR OF THE LENGTH N+1
*  IA  WN12(N+1) WORK VECTOR OF THE LENGTH N+1
*  IA  WN13(N+1) WORK VECTOR OF THE LENGTH N+1
*  IO  ISPACE AN INFORMATION ON SPACE OBTAINED DURING THE PROCESS
*       OF THE FACTORIZATION.
*      ISPACE=0    THE FACTORIZATION HAS TERMINATED NORMALLY
*      ISPACE=1    INSUFFICIENT SPACE AVAILABLE
*
* METHOD :
* S.C.EISENSTAT,M.C.GURSKY,M.H.SCHULTZ,A.H.SHERMAN:YALE SPARSE MATRIX
* PACKAGE I. THE SYMMETRIC CODES,ACM TRANS. ON MATH. SOFTWARE.
*
* NOTE: TYPE OF SL CHANGED FOR THE UFO APPLICATION.
*
      SUBROUTINE MXSPCI(N,ML,MMAX,PA,SA,PSL,SL,PERM,INVP,
     *            WN11,WN12,WN13,ISPACE)
      INTEGER N,MMAX,PA(*),SA(*),PSL(*)
      INTEGER PERM(*),INVP(*),WN11(*),WN12(*),WN13(*)
      INTEGER ISPACE,I,INZ,J,JSTOP,JSTRT,K,KNZ,KXSUB,MRGK,LMAX,ML
      INTEGER NABOR,NODE,NP1,NZBEG,NZEND,RCHM,MRKFLG,M
      DOUBLE PRECISION SL(*)
      NZBEG=1
      NZEND=0
      PSL(1)=1
      DO 100 K=1,N
      WN11(K)=0
      WN13(K)=0
100   CONTINUE
      NP1=N+1
      DO 1500 K=1,N
      KNZ=0
      MRGK=WN11(K)
      MRKFLG=0
      WN13(K)=K
      IF(MRGK.NE.0) WN13(K)=WN13(MRGK)
      SL(K)=NZEND
      NODE=PERM(K)
      JSTRT=PA(NODE)
      JSTOP=PA(NODE+1)-1
      IF(JSTRT.GT.JSTOP) GO TO 1500
      WN12(K)=NP1
      DO 300 J=JSTRT,JSTOP
      NABOR=SA(J)
      IF(NABOR.EQ.NODE) GO TO 300
      NABOR=INVP(NABOR)
      IF(NABOR.LE.K) GO TO 300
      RCHM=K
200   M=RCHM
      RCHM=WN12(M)
      IF(RCHM.LE.NABOR) GO TO 200
      KNZ=KNZ+1
      WN12(M)=NABOR
      WN12(NABOR)=RCHM
      IF(WN13(NABOR).NE.WN13(K)) MRKFLG=1
300   CONTINUE
      LMAX=0
      IF(MRKFLG.NE.0.OR.MRGK.EQ.0) GO TO 350
      IF(WN11(MRGK).NE.0) GO TO 350
      SL(K)=SL(MRGK)+1
      KNZ=PSL(MRGK+1)-(PSL(MRGK)+1)
      GO TO 1400
350   I=K
400   I=WN11(I)
      IF(I.EQ.0) GO TO 800
      INZ=PSL(I+1)-(PSL(I)+1)
      JSTRT=SL(I)+1
      JSTOP=SL(I)+INZ
      IF(INZ.LE.LMAX) GO TO 500
      LMAX=INZ
      SL(K)=JSTRT
500   RCHM=K
      DO 700 J=JSTRT,JSTOP
      NABOR=SL(N+1+J)
600   M=RCHM
      RCHM=WN12(M)
      IF(RCHM.LT.NABOR) GO TO 600
      IF(RCHM.EQ.NABOR) GO TO 700
      KNZ=KNZ+1
      WN12(M)=NABOR
      WN12(NABOR)=RCHM
      RCHM=NABOR
700   CONTINUE
      GO TO 400
800   IF(KNZ.EQ.LMAX) GO TO 1400
      IF(NZBEG.GT.NZEND) GO TO 1200
      I=WN12(K)
      DO 900 JSTRT=NZBEG,NZEND
      IF(SL(N+1+JSTRT)-I .GE.0) THEN
        IF(SL(N+1+JSTRT).EQ.I) THEN
          GO TO 1000
        ELSE
          GO TO 1200
        END IF
      END IF
900   CONTINUE
      GO TO 1200
1000  SL(K)=JSTRT
      DO 1100 J=JSTRT,NZEND
      IF (SL(N+1+J).NE.I) GO TO 1200
      I=WN12(I)
      IF(I.GT.N) GO TO 1400
1100  CONTINUE
      NZEND=JSTRT-1
1200  NZBEG=NZEND+1
      NZEND=NZEND+KNZ
C
C A VARIANT IS USED WHEN CALLED SO THAT SL(X)=A(NB+X)
C
      IF(NZEND.GE.MMAX-N-1) GO TO 1600
      I=K
      DO 1300 J=NZBEG,NZEND
      I=WN12(I)
      SL(N+1+J)=I
      WN13(I)=K
1300  CONTINUE
      SL(K)=NZBEG
      WN13(K)=K
1400     IF(KNZ.GT.1) THEN
            KXSUB=SL(K)
            I=SL(N+1+KXSUB)
            WN11(K)=WN11(I)
            WN11(I)=K
      END IF
      PSL(K+1)=PSL(K)+KNZ
1500  CONTINUE
      SL(N)=SL(N)+1
      SL(N+1)=SL(N)
      ML=N+SL(N+1)
      ISPACE=0
      RETURN
1600  ISPACE=1
      RETURN
      END
* SUBROUTINE MXSPCM                ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A GIVEN VECTOR X BY A SPARSE SYMMETRIC POSITIVE
* DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED
* BY THE SUBROUTINE MXSPGN. FACTORS ARE STORED IN THE COMPACT FORM.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(MMAX)      FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXSPGN.IT CONTAINS THE NUMERICAL VALUES OF THE
*         FACTORS STORED IN THE COMPACT FORM ACCORDING TO THE
*         INFORMATION IN THE VECTORS PSL,SL.
*  II  PSL(N+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX
*  II  SL(MMAX) INDICES OF THE COMPACT SCHEME OF THE FACTORS.
*  RU  X(N)  ON INPUT THE GIVEN VECTOR, ON OUTPUT THE RESULT
*         OF MULTIPLICATION.
*  RA  RN01(N) AUXILIARY VECTOR.
*  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)*X. IF JOB>0 THEN
*         X:=TRANS(L)*X. IF JOB<0 THEN X:=L*X.
*
      SUBROUTINE MXSPCM(N,A,PSL,SL,X,RN01,JOB)
      INTEGER N
      INTEGER PSL(*),SL(*),JOB
      DOUBLE PRECISION  A(*),X(*),RN01(*),ZERO
      PARAMETER(ZERO=0.0D0)
      INTEGER I,J,IS
      DO 50 I=1,N
      RN01(I)=ZERO
 50   CONTINUE
C
C     FIRST PHASE:X=TRANS(L)*X
C
      IF (JOB.GE.0) THEN
                    DO 300 I=1,N
                    IS=SL(I)+N+1
                    DO 200 J=PSL(I)+I,PSL(I+1)+I-1
                    X(I)=X(I)+A(J)*X(SL(IS))
                    IS=IS+1
200                 CONTINUE
300                 CONTINUE
      END IF
C
C     SECOND PHASE:X=D*X
C
      IF (JOB.EQ.0) THEN
                    DO 400 I=1,N
                    X(I) = X(I)*A(PSL(I)+I-1)
400   CONTINUE
      END IF
C
C     THIRD PHASE:X=L*X
C
      IF (JOB.LE.0) THEN
                    DO 600 I=N,1,-1
                    IS=SL(I)+N+1
                    DO 500 J=PSL(I)+I,PSL(I+1)+I-1
                    RN01(SL(IS))=RN01(SL(IS))+A(J)*X(I)
                    IS=IS+1
500                 CONTINUE
600                 CONTINUE
                    DO 700 I=1,N
                    X(I)=RN01(I)+X(I)
700                 CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXSPCT                 ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* REWRITE SYMMETRIC MATRIX INTO THE PERMUTED FACTORIZED COMPACT SCHEME.
* MOIDIFIED VERSION FOR THE USE WITH MXSPCJ.
*
* PARAMETERS:
*  II  N  SIZE OF THE SYSTEM SOLVED.
*  II  NB NUMBER OF NONZEROS IN THE UPPER TRIANGLE OF THE ORIGINAL
*         MATRIX.
*  II  ML SIZE OF THE COMPACT FACTOR.
*  II  MMAX DECLARED LENGTH OF THE ARRAYS JA,A.
*  RU  A(MMAX) NUMERICAL VALUES OF THE SPARSE HESSIAN APPROXIMATION
*              STORED AT THE POSITIONS 1, ...,NB.
*  IU  JA(MMAX) INDICES OF THE NONZERO ELEMENTS OF THE HESSIAN MATRIX IN
*             THE PACKED ROW FORM AT THE FIRST NB POSITIONS.
*             NEW POSITIONS
*             IN THE PERMUTED FACTOR STORED IN A(NB+1), ..., A(2*NB),
*             INDICES OF COMPACT SCHEME IN A(2*NB+1), ..., A(2*NB+ML).
*  II  PSL(N+1)  POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR
*             FACTOR OF THE HESSIAN APPROXIMATION.
*  IO  ITERM  ERROR FLAG. IF ITERM < 0 - LACK OF SPACE IN JA.
*
*
      SUBROUTINE MXSPCT(N,NB,ML,MMAX,A,JA,PSL,ITERM)
      INTEGER N,NB,ML,MMAX,JA(*)
      INTEGER PSL(*),ITERM
      DOUBLE PRECISION A(*)
      INTEGER I,J
C
C WN11 CONTAINS BEGINNINGS OF THE FACTOR ROWS
C
      ITERM=0
C
C LACK OF SPACE
C
      IF(MMAX.LE.NB+PSL(N+1)+N-1) THEN
        ITERM=-43
        RETURN
      END IF
      IF(MMAX.LE.2*NB+ML) THEN
        ITERM=-44
        RETURN
      END IF
C
      DO 50 I=NB+1,NB+PSL(N+1)+N-1
        A(I)=0.0D 0
 50   CONTINUE
      DO 100 I=NB+ML+1,2*NB+ML
        J=JA(I)
        A(NB+J)=A(I-NB-ML)
 100  CONTINUE
      RETURN
      END
* SUBROUTINE MXSPCX                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* STORED IN THE COMPACT SCHEME. THIS FACTORIZATION CAN BE OBTAINED
* USING THE SUBROUTINE MXSSFN.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  II  ND  NUMBER OF DENSE ROWS.
*  RI  A(MMAX)  FACTORS L,D OF THE FACTORIZATION A+E=L*D*TRANS(L)
*                STORED USING THE COMPACT SCHEME OF STORING.
*  II  PSL(N+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX
*  II  SL(MMAX)  ARRAY OF COLUMN INDICES OF THE FACTORS L AND D
*         STORED USING THE COMPACT SCHEME.
*  II  PERM(N)  PERMUTATION VECTOR
*  RI  U(N*ND)  DENSE ROWS.
*  RI  D(ND)  SCALING MATRIX.
*  RO  B(ND*(ND+1)/2))  DECOMPOSITION OF THE UPDATED MATRIX.
*  RA  X(N)  AUXILIARY ARRAY.
*  RA  Y(N)  AUXILIARY ARRAY.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXSPCX(N,ND,A,PSL,SL,PERM,U,D,B,X,Y)
      INTEGER N,ND
      DOUBLE PRECISION A(*),U(*),D(*),B(*),X(*),Y(*)
      INTEGER PSL(*),SL(*),PERM(*)
      INTEGER I,J,INF,INDII,INDIJ
      DOUBLE PRECISION ALF,BET,MXVDOT
      DO 2 I=1,ND
      INDII=I*(I+1)/2
      CALL MXVCOP(N,U((I-1)*N+1),Y)
      CALL MXVSFP(N,PERM,Y,X)
      CALL MXSPCB(N,A,PSL,SL,Y,0)
      CALL MXVSBP(N,PERM,Y,X)
      DO 1 J=1,I
      INDIJ=I*(I-1)/2+J
      B(INDIJ)=MXVDOT(N,U((J-1)*N+1),Y)
    1 CONTINUE
      B(INDII)=B(INDII)+D(I)
    2 CONTINUE
      ALF=1.0D-16
      CALL MXDPGF(ND,B,INF,ALF,BET)
      RETURN
      END
* SUBROUTINE MXSPCZ                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* STORED IN THE COMPACT SCHEME. THIS FACTORIZATION CAN BE OBTAINED
* USING THE SUBROUTINE MXSSFN.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  II  ND  NUMBER OF DENSE ROWS.
*  RI  A(MMAX)  FACTORS L,D OF THE FACTORIZATION A+E=L*D*TRANS(L)
*                STORED USING THE COMPACT SCHEME OF STORING.
*  II  PSL(N+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX
*  II  SL(MMAX)  ARRAY OF COLUMN INDICES OF THE FACTORS L AND D
*         STORED USING THE COMPACT SCHEME.
*  II  PERM(N)  PERMUTATION VECTOR
*  RI  U(N*ND)  DENSE ROWS.
*  RA  D(ND)  AUXILIARY ARRAY.
*  RO  B(ND*(ND+1)/2))  DECOMPOSITION OF THE UPDATED MATRIX.
*  RA  X(N)  AUXILIARY ARRAY.
*  RA  Y(N)  VECTOR MULTIPLIED BY THE INVERSE OF THE UPDATED MATRIX.
*  RA  Z(N)  AUXILIARY ARRAY.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXSPCZ(N,ND,A,PSL,SL,PERM,U,D,B,X,Y,Z)
      INTEGER N,ND
      DOUBLE PRECISION A(*),U(*),D(*),B(*),X(*),Y(*),Z(*)
      INTEGER PSL(*),SL(*),PERM(*)
      INTEGER I
      DOUBLE PRECISION MXVDOT
      CALL MXVCOP(N,Y,Z)
      CALL MXVSFP(N,PERM,Y,X)
      CALL MXSPCB(N,A,PSL,SL,Y,0)
      CALL MXVSBP(N,PERM,Y,X)
      IF (ND.LE.0) RETURN
      DO 1 I=1,ND
      D(I)=MXVDOT(N,U((I-1)*N+1),Y)
    1 CONTINUE
      CALL MXDPGB(ND,B,D,0)
      CALL MXVCOP(N,Z,Y)
      DO 2 I=1,ND
      CALL MXVDIR(N,-D(I),U((I-1)*N+1),Y,Y)
    2 CONTINUE
      CALL MXVSFP(N,PERM,Y,X)
      CALL MXSPCB(N,A,PSL,SL,Y,0)
      CALL MXVSBP(N,PERM,Y,X)
      RETURN
      END
* SUBROUTINE MXSPIB                ALL SYSTEMS                94/12/01
C PORTABILITY : ALL SYSTEMS
C 94/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING INCOMPLETE FACTORIZATION
* A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXSPIF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(M) INCOMPLETE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXSPGF.
*  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
*  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
*  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
*         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
*         EQUATIONS.
*  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
*         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXSPIB(N,A,IA,JA,X,JOB)
      INTEGER N,IA(*),JA(*),JOB
      DOUBLE PRECISION A(*),X(*)
      INTEGER I,K,IK,KK,LL
      DOUBLE PRECISION CON
      PARAMETER (CON=1.0D 60)
      IF (JOB .GE. 0) THEN
C
C     PHASE 1 : X:=L**(-1)*X
C
      KK=1
      DO 3 K=1,N
      LL=IA(K+1)
      DO 2 IK=KK+1,LL-1
      I=JA(IK)
      X(I)=X(I)-A(IK)*X(K)
    2 CONTINUE
      KK=LL
    3 CONTINUE
      END IF
      IF (JOB .EQ. 0) THEN
C
C     PHASE 2 : X:=D**(-1)*X
C
      DO 4 K=1,N
      X(K)=X(K)/A(IA(K))
    4 CONTINUE
      END IF
      IF (JOB .LE. 0) THEN
C
C     PHASE 3 : X:=TRANS(L)**(-1)*X
C
      LL=IA(N+1)
      DO 6 K=N,1,-1
      KK=IA(K)
      DO 5 IK=KK+1,LL-1
      I=JA(IK)
      X(K)=X(K)-A(IK)*X(I)
    5 CONTINUE
      LL=KK
    6 CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXSPIF                ALL SYSTEMS                94/12/01
C PORTABILITY : ALL SYSTEMS
C 94/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* INCOMPLETE FACTORIZATION A+E=L*D*TRANS(L) OF A SPARSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E WHERE D AND E ARE DIAGONAL POSITIVE
* DEFINITE MATRICES AND L IS A LOWER TRIANGULAR MATRIX.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(M)  ON INPUT A GIVEN SPARSE SYMMETRIC (USUALLY POSITIVE
*         DEFINITE) MATRIX. ON OUTPUT THE COMPUTED INCOMPLETE
*         FACTORIZATION A+E=L*D*TRANS(L).
*  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
*  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
*  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
*         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
*         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
*         INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
*         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
*         PROCESS.
*  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
*         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
*         FACTORIZATION PROCESS (IF INF>0).
*  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
*
* METHOD :
* P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
* LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
* PP. 311-350.
*
      SUBROUTINE MXSPIF(N,A,IA,JA,INF,ALF,TAU,X)
      INTEGER N,IA(*),JA(*),INF
      DOUBLE PRECISION A(*),ALF,TAU,X(*)
      DOUBLE PRECISION BET,GAM,DEL,RHO,SIG,TOL,POM
      INTEGER I,J,K,L,IJ,KJ,KK,KN
      L=0
      INF=0
      TOL=ALF
C
C     ESTIMATION OF THE MATRIX NORM
C
      ALF=0.0D 0
      BET=0.0D 0
      GAM=0.0D 0
      TAU=0.0D 0
      CALL MXVSET(N,0.0D 0,X)
      DO 1 K=1,N
      BET=MAX(BET,ABS(A(IA(K))))
    1 CONTINUE
      DO 2 K=1,IA(N+1)-1
      GAM=MAX(GAM,ABS(A(K)))
    2 CONTINUE
      BET=MAX(TOL,BET,GAM/DBLE(N))
      DEL=TOL*BET
      DO 7 K=1,N
      KK=IA(K)
      KN=IA(K+1)-1
C
C     DETERMINATION OF A DIAGONAL CORRECTION
C
      SIG=A(KK)
      IF (ALF.GT.SIG) THEN
      ALF=SIG
      L=K
      END IF
      GAM=0.0D 0
      DO 3 KJ=KK+1,KN
      J=JA(KJ)
      X(J)=A(KJ)
      GAM=MAX(GAM,ABS(A(KJ)))
    3 CONTINUE
      GAM=GAM*GAM
      RHO=MAX(ABS(SIG),GAM/BET,DEL)
      IF (TAU.LT.RHO-SIG) THEN
      TAU=RHO-SIG
      INF=-1
      END IF
C
C     GAUSSIAN ELIMINATION
C
      A(KK)=RHO
      DO 5 I=K+1,N
      IF (X(I).EQ.0.0D 0) GO TO 5
      POM=X(I)/RHO
      DO 4 IJ=IA(I),IA(I+1)-1
      J=JA(IJ)
      IF (X(J).EQ.0.0D 0) GO TO 4
      A(IJ)=A(IJ)-POM*X(J)
    4 CONTINUE
      X(I)=POM
    5 CONTINUE
      DO 6 KJ=KK+1,KN
      J=JA(KJ)
      A(KJ)=X(J)
      X(J)=0.0D 0
    6 CONTINUE
    7 CONTINUE
      IF (L.GT.0.AND.ABS(ALF).GT.DEL) INF = L
      RETURN
      END
* SUBROUTINE MXSPIX                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* STORED IN THE COMPACT SCHEME. THIS FACTORIZATION CAN BE OBTAINED
* USING THE SUBROUTINE MXSSFN.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  II  ND  NUMBER OF DENSE ROWS.
*  RI  A(M) INCOMPLETE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXSPIF.
*  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
*  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
*  RI  U(N*ND)  DENSE ROWS.
*  RI  D(ND)  SCALING MATRIX.
*  RO  B(ND*(ND+1)/2))  DECOMPOSITION OF THE UPDATED MATRIX.
*  RA  Y(N)  AUXILIARY ARRAY.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXSPIX(N,ND,A,IA,JA,U,D,B,Y)
      INTEGER N,ND
      DOUBLE PRECISION A(*),U(*),D(*),B(*),Y(*)
      INTEGER IA(*),JA(*)
      INTEGER I,J,INF,INDII,INDIJ
      DOUBLE PRECISION ALF,BET,MXVDOT
      INDII=I*(I+1)/2
      DO 2 I=1,ND
      CALL MXVCOP(N,U((I-1)*N+1),Y)
      CALL MXSPIB(N,A,IA,JA,Y,0)
      DO 1 J=1,I
      INDIJ=I*(I-1)/2+J
      B(INDIJ)=MXVDOT(N,U((J-1)*N+1),Y)
    1 CONTINUE
      B(INDII)=B(INDII)+D(I)
    2 CONTINUE
      ALF=1.0D-16
      CALL MXDPGF(ND,B,INF,ALF,BET)
      RETURN
      END
* SUBROUTINE MXSPIZ                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* STORED IN THE COMPACT SCHEME. THIS FACTORIZATION CAN BE OBTAINED
* USING THE SUBROUTINE MXSSFN.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  II  ND  NUMBER OF DENSE ROWS.
*  RI  A(M) INCOMPLETE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXSPIF.
*  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
*  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
*  RI  U(N*ND)  DENSE ROWS.
*  RA  D(ND)  AUXILIARY ARRAY.
*  RO  B(ND*(ND+1)/2))  DECOMPOSITION OF THE UPDATED MATRIX.
*  RA  Y(N)  VECTOR MULTIPLIED BY THE INVERSE OF THE UPDATED MATRIX.
*  RA  Z(N)  AUXILIARY ARRAY.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXSPIZ(N,ND,A,IA,JA,U,D,B,Y,Z)
      INTEGER N,ND
      DOUBLE PRECISION A(*),U(*),D(*),B(*),Y(*),Z(*)
      INTEGER IA(*),JA(*)
      INTEGER I
      DOUBLE PRECISION MXVDOT
      CALL MXVCOP(N,Y,Z)
      CALL MXSPIB(N,A,IA,JA,Y,0)
      IF (ND.LE.0) RETURN
      DO 1 I=1,ND
      D(I)=MXVDOT(N,U((I-1)*N+1),Y)
    1 CONTINUE
      CALL MXDPGB(ND,B,D,0)
      CALL MXVCOP(N,Z,Y)
      DO 2 I=1,ND
      CALL MXVDIR(N,-D(I),U((I-1)*N+1),Y,Y)
    2 CONTINUE
      CALL MXSPIB(N,A,IA,JA,Y,0)
      RETURN
      END
* SUBROUTINE MXSRMD               ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF TRANSPOSE OF A DENSE RECTANGULAR MATRIX A BY
* A VECTOR X AND ADDITION OF A SCALED VECTOR ALF*Y.
*
* PARAMETERS :
*  II  N  NUMBER OF ROWS OF THE MATRIX A.
*  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
*  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
*  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
*  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
*  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
*  RI  X(N)  INPUT VECTOR.
*  RI  ALF  SCALING FACTOR.
*  RI  Y(NA)  INPUT VECTOR.
*  RO  Z(NA)  OUTPUT VECTOR EQUAL TO TRANS(A)*X+ALF*Y.
*
      SUBROUTINE MXSRMD(NA,A,IA,JA,X,ALF,Y,Z)
      INTEGER NA,IA(*),JA(*)
      DOUBLE PRECISION A(*),X(*),ALF,Y(*),Z(*)
      DOUBLE PRECISION TEMP
      INTEGER I,J,K,L,JP
      DO 2 I=1,NA
      K=IA(I)
      L=IA(I+1)-K
      TEMP=ALF*Y(I)
      DO 1 J=1,L
      JP=JA(K)
      IF (JP.GT.0) TEMP=TEMP+A(K)*X(JP)
      K=K+1
    1 CONTINUE
      Z(I)=TEMP
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXSRMM               ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF TRANSPOSE OF A DENSE RECTANGULAR MATRIX A BY
* A VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF ROWS OF THE MATRIX A.
*  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
*  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
*  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
*  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
*  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(M)  OUTPUT VECTOR EQUAL TO TRANS(A)*X.
*
      SUBROUTINE MXSRMM(NA,A,IA,JA,X,Y)
      INTEGER NA,IA(*),JA(*)
      DOUBLE PRECISION A(*),X(*),Y(*)
      DOUBLE PRECISION TEMP
      INTEGER I,J,K,L,JP
      DO 2 I=1,NA
      K=IA(I)
      L=IA(I+1)-K
      TEMP=0.0D 0
      DO 1 J=1,L
      JP=JA(K)
      IF (JP.GT.0) TEMP=TEMP+A(K)*X(JP)
      K=K+1
    1 CONTINUE
      Y(I)=TEMP
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXSSMD                ALL SYSTEMS                93/12/01
C PORTABILITY : ALL SYSTEMS
C 93/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A SPARSE SYMMETRIC MATRIX A BY A VECTOR X
* AND ADDITION OF A SCALED VECTOR ALF*Y.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
*      PACKED FORM.
*  II  IA(N)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
*  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
*  RI  X(N)  INPUT VECTOR.
*  RI  ALF  SCALING FACTOR.
*  RI  Y(NA)  INPUT VECTOR.
*  RO  Z(NA)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
*
      SUBROUTINE MXSSMD(N,A,IA,JA,X,ALF,Y,Z)
      INTEGER N,IA(*),JA(*)
      DOUBLE PRECISION  A(*),X(*),ALF,Y(*),Z(*)
      INTEGER I,J,K,JSTRT,JSTOP
      DO 100 I=1,N
      Z(I)=ALF*Y(I)
100   CONTINUE
      JSTOP=0
      DO 300 I=1,N
        JSTRT=JSTOP+1
        JSTOP=IA(I+1)-1
        IF (JA(JSTRT).GT.0) THEN
          DO 200 J=JSTRT,JSTOP
            K=JA(J)
            IF(J.EQ.JSTRT) THEN
              Z(I)=Z(I)+A(J)*X(I)
            ELSE IF (K.GT.0) THEN
              Z(K)=Z(K)+A(J)*X(I)
              Z(I)=Z(I)+A(J)*X(K)
            END IF
200       CONTINUE
        END IF
300   CONTINUE
      RETURN
      END
* SUBROUTINE MXSSMM                ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A SPARSE SYMMETRIC MATRIX BY A VECTOR X.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  II  M  NUMBER OF NONZERO ELEMENTS OF THE MATRIX A.
*  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
*      PACKED FORM.
*  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
*  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y := A * X.
*
      SUBROUTINE MXSSMM(N,A,IA,JA,X,Y)
      INTEGER N,IA(*),JA(*)
      DOUBLE PRECISION  A(*),X(*),Y(*)
      INTEGER I,J,K,JSTRT,JSTOP
      DO 100 I=1,N
      Y(I)=0.0D 0
100   CONTINUE
      JSTOP=0
      DO 300 I=1,N
        JSTRT=JSTOP+1
        JSTOP=IA(I+1)-1
        IF (JA(JSTRT).GT.0) THEN
          DO 200 J=JSTRT,JSTOP
            K=JA(J)
            IF(J.EQ.JSTRT) THEN
              Y(I)=Y(I)+A(J)*X(I)
            ELSE IF (K.GT.0) THEN
              Y(K)=Y(K)+A(J)*X(I)
              Y(I)=Y(I)+A(J)*X(K)
            END IF
200       CONTINUE
        END IF
300   CONTINUE
      RETURN
      END
* SUBROUTINE MXSSMN                ALL SYSTEMS                89/12/01
C PORTABILITY : ALL SYSTEMS
C 89/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* THIS SUBROUTINE FINDS THE PERMUTATION VECTOR PERM FOR THE
* SPARSE SYMMETRIC MATRIX GIVEN IN THE VECTOR PAIR PA,SA.IT USES
* THE SO-CALLED NESTED DISSECTION METHOD.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  II  MMAX  LENGTH OF THE PRINCIPAL MATRIX VECTORS.
*  II  PA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
*  II  SA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
*  IO  PERM(N) PERMUTATION VECTOR.
*  IA  WN11(N+1) AUXILIARY VECTOR.
*  IA  WN12(N+1) AUXILIARY VECTOR.
*  IA  WN13(N+1) AUXILIARY VECTOR.
*
* METHOD :
* NESTED DISSECTION METHOD
*
      SUBROUTINE MXSSMN(N,PA,SA,PERM,WN11,WN12,WN13)
      INTEGER N
      INTEGER PA(*),SA(*),PERM(*)
      INTEGER WN11(*),WN12(*),WN13(*)
      INTEGER I,J,K,NUM,ROOT,NLVL,LVLEND,LBEGIN,ICS
      INTEGER NN,N1,MINDEG,N2,LVSIZE,NDEG,NUNLVL,MIDLVL
      INTEGER TEMP,NPUL,NSEP,I1,I2,I3,I4,J1,J2
      DO 100 I = 1, N
      WN11(I) = 1
  100 CONTINUE
      NUM= 0
      DO 2000 I = 1, N
  200 IF ( WN11(I) .EQ. 0 ) GO TO 2000
      ROOT = I
      WN11(ROOT) = 0
      WN13(1) = ROOT
      NLVL = 0
      LVLEND = 0
      ICS = 1
  300 LBEGIN = LVLEND + 1
      LVLEND = ICS
      NLVL = NLVL + 1
      WN12(NLVL) = LBEGIN
      DO 500 K = LBEGIN, LVLEND
      NN = WN13(K)
      DO 400 J=PA(NN),PA(NN+1)-1
      N2 = SA(J)
      IF(N2.EQ.NN) GO TO 400
      IF  (WN11(N2).EQ.0) GO TO 400
      ICS = ICS + 1
      WN13(ICS) = N2
      WN11(N2) = 0
  400 CONTINUE
  500 CONTINUE
      LVSIZE=ICS-LVLEND
      IF(LVSIZE.GT.0) GO TO 300
      WN12(NLVL+1) = LVLEND + 1
      DO 600 K = 1, ICS
      NN = WN13(K)
      WN11(NN) = 1
  600 CONTINUE
      ICS = WN12(NLVL+1) - 1
      IF  ( NLVL .EQ. 1 .OR. NLVL .EQ. ICS )GO TO 1470
  700 J1 = WN12(NLVL)
      MINDEG = ICS
      ROOT = WN13(J1)
      IF  ( ICS .EQ. J1 ) GO TO 1000
      DO 900 J = J1, ICS
      NN = WN13(J)
      NDEG = 0
      DO 800 K=PA(NN),PA(NN+1)-1
      N1 = SA(K)
      IF(N1.EQ.NN) GO TO 800
      IF  ( WN11(N1)  .GT. 0 ) NDEG = NDEG + 1
  800 CONTINUE
      IF  ( NDEG .GE. MINDEG  ) GO TO 900
      ROOT = NN
      MINDEG = NDEG
  900 CONTINUE
 1000 CONTINUE
      WN11(ROOT) = 0
      WN13(1) = ROOT
      NUNLVL = 0
      LVLEND = 0
      ICS = 1
 1100 LBEGIN = LVLEND + 1
      LVLEND = ICS
      NUNLVL = NUNLVL + 1
      WN12(NUNLVL) = LBEGIN
      DO 1300 K = LBEGIN, LVLEND
      NN = WN13(K)
      DO 1200 J=PA(NN),PA(NN+1)-1
      N2 = SA(J)
      IF(N2.EQ.NN) GO TO 1200
      IF  (WN11(N2).EQ.0) GO TO 1200
      ICS = ICS + 1
      WN13(ICS) = N2
      WN11(N2) = 0
 1200 CONTINUE
 1300 CONTINUE
      LVSIZE=ICS-LVLEND
      IF(LVSIZE.GT.0) GO TO 1100
      WN12(NUNLVL+1) = LVLEND + 1
      DO 1400 K = 1, ICS
      NN = WN13(K)
      WN11(NN) = 1
 1400 CONTINUE
 1450 CONTINUE
      IF(NUNLVL .LE.NLVL) GO TO 1470
      NLVL=NUNLVL
      IF(NLVL.LT.ICS) GO TO 700
1470  CONTINUE
      IF ( NLVL .GE. 3 ) GO TO 1600
      NSEP = WN12(NLVL+1) - 1
      DO 1500 K = 1, NSEP
      NN = WN13(K)
      PERM(NUM+K) = NN
      WN11(NN) = 0
 1500 CONTINUE
      GO TO 1950
 1600 MIDLVL = (NLVL+2)/2
      I3 = WN12(MIDLVL)
      I1 = WN12(MIDLVL + 1)
      I4 = I1 - 1
      I2 = WN12(MIDLVL+2) - 1
      DO 1700 K = I1, I2
      NN = WN13(K)
      PA(NN) = - PA(NN)
 1700 CONTINUE
      NSEP = 0
      DO 1800 K = I3, I4
      NN = WN13(K)
      J1 = PA(NN)
      J2 = IABS(PA(NN+1)) - 1
      DO 1750 J = J1, J2
      N2 = SA(J)
      IF(N2.EQ.NN) GO TO 1750
      IF ( PA(N2) .GT. 0 ) GO TO 1750
      NSEP = NSEP + 1
      PERM(NSEP+NUM) = NN
      WN11(NN) = 0
      GO TO 1800
 1750 CONTINUE
 1800 CONTINUE
      DO 1900 K = I1, I2
      NN = WN13(K)
      PA(NN) = - PA(NN)
 1900 CONTINUE
 1950 CONTINUE
      NUM = NUM + NSEP
      IF ( NUM .GE. N ) GO TO 2100
      GO TO 200
 2000 CONTINUE
 2100 CONTINUE
      IF (N.LT.2) GO TO 2300
      NPUL = N/2
      DO 2200 I=1,NPUL
      TEMP=PERM(I)
      PERM(I)=PERM(N-I+1)
      PERM(N-I+1)=TEMP
2200  CONTINUE
2300  CONTINUE
      RETURN
      END
* SUBROUTINE MXSTG1                ALL SYSTEMS                89/12/01
C PORTABILITY : ALL SYSTEMS
C 89/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* WIDTHENING THE PACKED FORM OF THE VECTORS IA, JA OF THE SPARSE MATRIX
*
* PARAMETERS :
*  II  N ORDER OF THE SPARSE MATRIX.
*  IU  M NUMBER OF NONZERO ELEMENTS IN THE MATRIX.
*  II  MMAX LENGTH OF THE ARRAY JA.
*  II  IA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
*  II  JA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
*  IA  PD(N+1) AUXILIARY VECTOR.
*  IA  WN11(N+1) AUXILIARY VECTOR.
*
      SUBROUTINE MXSTG1(N,M,IA,JA,PD,WN11)
      INTEGER N,M
      INTEGER IA(*),PD(*),JA(*),WN11(*)
      INTEGER I,J,L1,L,K
C
C UPPER TRIANGULAR INFORMATION TO THE AUXILIARY ARRAY
C
      L1=IA(1)
      DO 100 I=1,N
      L=L1
      L1=IA(I+1)
      WN11(I)=L1-L
100   CONTINUE
C
C LOWER TRIANGULAR INFORMATION TO THE AUXILIARY ARRAY
C
      DO 300  I=1,N
        DO 200  J=IA(I)+1,IA(I+1)-1
        K=ABS(JA(J))
        WN11(K)=WN11(K)+1
200     CONTINUE
300   CONTINUE
C
C BY PARTIAL SUMMING WE GET POINTERS OF THE WIDE STRUCTURE
C WN11(I) POINTS AT THE END OF THE ROW I
C
      L=0
      DO 400  I=2,N
      WN11(I)=WN11(I)+WN11(I-1)
400   CONTINUE
C
C DEFINE LENGTH OF THE WITHENED STRUCTURE
C
      M=WN11(N)
C
C SHIFT OF UPPER TRIANGULAR ROWS
C
      PD(1)=1
      DO 600  I=N,1,-1
      L=WN11(I)
      PD(I+1)=L+1
        DO 500  J=IA(I+1)-1,IA(I),-1
        JA(L)=JA(J)
        L=L-1
500     CONTINUE
600   CONTINUE
C
C FORMING OF THE LOWER TRIANGULAR PART
C
      DO 800  I=1,N
        DO 700  J=WN11(I)+IA(I)+2-IA(I+1),WN11(I)
        K=ABS(JA(J))
        JA(PD(K))=I
        PD(K)=PD(K)+1
700     CONTINUE
800   CONTINUE
      DO 900  I=1,N
      IA(I+1)=WN11(I)+1
900   CONTINUE
      RETURN
      END
* SUBROUTINE MXSTL1                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* PACKING OF THE WIDTHENED FORM OF THE VECTORS IA, JA OF THE SPARSE
* MATRIX
*
* PARAMETERS :
*  II  N ORDER OF THE SPARSE MATRIX.
*  IU  M NUMBER OF NONZERO ELEMENTS IN THE MATRIX.
*  II  MMAX LENGTH OF THE ARRAY JA.
*  IU  IA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
*  IU  JA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
*  IA  PD(N+1) AUXILIARY VECTOR.
*
      SUBROUTINE MXSTL1(N,M,IA,JA,PD)
      INTEGER N,M
      INTEGER IA(*),PD(*),JA(*)
      INTEGER I,J,L,JSTRT,JSTOP
      L=1
C
C PD DEFINITION
C
      JSTOP=0
      DO 60 I=1,N
        JSTRT=JSTOP+1
        JSTOP=IA(I+1)-1
        DO 50 J=JSTRT,JSTOP
          IF(ABS(JA(J)).EQ.I) THEN
            PD(I)=J
            GO TO 60
          END IF
   50    CONTINUE
   60 CONTINUE
C
C REWRITE THE STRUCTURE
C
      DO 200 I=1,N
        DO 100 J=PD(I),IA(I+1)-1
        JA(L)=JA(J)
        L=L+1
100     CONTINUE
      IA(I+1)=L
200   CONTINUE
      IA(1)=1
C
C SET THE LENGTH OF THE PACKED STRUCTURE
C
      M=L-1
      RETURN
      END
* SUBROUTINE MXSTL2                ALL SYSTEMS                90/12/01
C PORTABILITY : ALL SYSTEMS
C 90/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* PACKING OF THE WIDTHENED FORM OF THE VECTORS A,IA,JA OF THE SPARSE
* MATRIX
*
* PARAMETERS :
*  II  N ORDER OF THE SPARSE MATRIX.
*  IU  M NUMBER OF NONZERO ELEMENTS IN THE MATRIX.
*  II  MMAX LENGTH OF THE ARRAY JA.
*  RU  A(MMAX) VECTOR OF NUMERICAL VALUES OF THE MATRIX BEING SHRINKED.
*  IU  IA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
*  IU  JA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
*  IA  PD(N+1) AUXILIARY VECTOR.
*
      SUBROUTINE MXSTL2(N,M,A,IA,JA,PD)
      INTEGER N,M
      INTEGER IA(*),PD(*),JA(*)
      DOUBLE PRECISION A(*)
      INTEGER I,J,L,JSTRT,JSTOP
      L=1
C
C PD DEFINITION
C
      JSTOP=0
      DO 60 I=1,N
        JSTRT=JSTOP+1
        JSTOP=IA(I+1)-1
        DO 50 J=JSTRT,JSTOP
          IF(ABS(JA(J)).EQ.I) THEN
            PD(I)=J
            GO TO 60
          END IF
   50    CONTINUE
   60 CONTINUE
C
C REWRITE THE STRUCTURE
C
      DO 200 I=1,N
        DO 100 J=PD(I),IA(I+1)-1
        JA(L)=JA(J)
        A(L)=A(J)
        L=L+1
100     CONTINUE
      IA(I+1)=L
200   CONTINUE
      IA(1)=1
C
C SET THE LENGTH OF THE PACKED STRUCTURE
C
      M=L-1
      RETURN
      END
* SUBROUTINE MXVCOP                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COPYING OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
*
      SUBROUTINE MXVCOP(N,X,Y)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*)
      INTEGER          I
      DO 10 I = 1,N
          Y(I) = X(I)
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVDIF                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR DIFFERENCE.
*
* PARAMETERS :
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X - Y.
*
      SUBROUTINE MXVDIF(N,X,Y,Z)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*),Z(*)
      INTEGER          I
      DO 10 I = 1,N
          Z(I) = X(I) - Y(I)
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVDIR                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR AUGMENTED BY THE SCALED VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  A  SCALING FACTOR.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= Y + A*X.
*
      SUBROUTINE MXVDIR(N,A,X,Y,Z)
      DOUBLE PRECISION A
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*),Z(*)
      INTEGER          I
      DO 10 I = 1,N
          Z(I) = Y(I) + A*X(I)
   10 CONTINUE
      RETURN
      END
* FUNCTION MXVDOT                  ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DOT PRODUCT OF TWO VECTORS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RR  MXVDOT  VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y.
*
      FUNCTION MXVDOT(N,X,Y)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*),MXVDOT
      DOUBLE PRECISION TEMP
      INTEGER          I
      TEMP = 0.0D0
      DO 10 I = 1,N
          TEMP = TEMP + X(I)*Y(I)
   10 CONTINUE
      MXVDOT = TEMP
      RETURN
      END
* SUBROUTINE MXVICP             ALL SYSTEMS                   93/12/01
* PORTABILITY : ALL SYSTEMS
* 93/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COPYING OF AN INTEGER VECTOR.
*
* PARAMETERS :
*  II  N DIMENSION OF THE INTEGER VECTOR.
*  II  IX(N)  INPUT INTEGER VECTOR.
*  IO  IY(N)  OUTPUT INTEGER VECTOR SUCH THAT IY(I):= IX(I) FOR ALL I.
*
      SUBROUTINE MXVICP(N,IX,IY)
      INTEGER N,IX(N),IY(N)
      INTEGER I
      DO 1 I=1,N
      IY(I)=IX(I)
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXVINS             ALL SYSTEMS                   90/12/01
* PORTABILITY : ALL SYSTEMS
* 90/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* INITIATION OF THE INTEGER VECTOR.
*
* PARAMETERS :
*  II  N DIMENSION OF THE INTEGER VECTOR.
*  II  IP  INTEGER PARAMETER.
*  IO  IX(N)  INTEGER VECTOR SUCH THAT IX(I)=IP FOR ALL I.
*
      SUBROUTINE MXVINS(N,IP,IX)
      INTEGER          IP,N
      INTEGER          IX(*)
      INTEGER          I
      DO 10 I = 1,N
          IX(I) = IP
   10 CONTINUE
      RETURN
      END
* FUNCTION MXVMAX             ALL SYSTEMS                   91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* L-INFINITY NORM OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RR  MXVMAX  L-INFINITY NORM OF THE VECTOR X.
*
      FUNCTION MXVMAX(N,X)
      INTEGER N
      DOUBLE PRECISION X(*),MXVMAX
      INTEGER I
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D 0)
      MXVMAX=ZERO
      DO 1 I=1,N
      MXVMAX=MAX(MXVMAX,ABS(X(I)))
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXVMUL             ALL SYSTEMS                   89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR IS PREMULTIPLIED BY THE POWER OF A DIAGONAL MATRIX.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  D(N)  DIAGONAL MATRIX STORED AS A VECTOR WITH N ELEMENTS.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:=(D**K)*X.
*  II  K  INTEGER EXPONENT.
*
      SUBROUTINE MXVMUL(N,D,X,Y,K)
      INTEGER          K,N
      DOUBLE PRECISION D(*),X(*),Y(*)
      INTEGER          I
      IF (K.EQ.0) THEN
          CALL MXVCOP(N,X,Y)
      ELSE IF (K.EQ.1) THEN
          DO 10 I = 1,N
              Y(I) = X(I)*D(I)
   10     CONTINUE
      ELSE IF (K.EQ.-1) THEN
          DO 20 I = 1,N
              Y(I) = X(I)/D(I)
   20     CONTINUE
      ELSE
          DO 30 I = 1,N
              Y(I) = X(I)*D(I)**K
   30     CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE MXVNEG                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CHANGE THE SIGNS OF VECTOR ELEMENTS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= - X.
*
      SUBROUTINE MXVNEG(N,X,Y)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*)
      INTEGER          I
      DO 10 I = 1,N
          Y(I) = -X(I)
   10 CONTINUE
      RETURN
      END
* FUNCTION  MXVNOR               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EUCLIDEAN NORM OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RR  MXVNOR  EUCLIDEAN NORM OF X.
*
      FUNCTION MXVNOR(N,X)
      INTEGER N
      DOUBLE PRECISION X(*),MXVNOR
      DOUBLE PRECISION DEN,POM
      INTEGER I
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
      DEN = ZERO
      DO 10 I = 1,N
          DEN = MAX(DEN,ABS(X(I)))
   10 CONTINUE
      POM = ZERO
      IF (DEN.GT.ZERO) THEN
          DO 20 I = 1,N
              POM = POM + (X(I)/DEN)**2
   20     CONTINUE
      END IF
      MXVNOR = DEN*SQRT(POM)
      RETURN
      END
* SUBROUTINE MXVSAV                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RU  Y(N)  UPDATE VECTOR WHERE Y:= X - Y.
*
      SUBROUTINE MXVSAV(N,X,Y)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*)
      DOUBLE PRECISION TEMP
      INTEGER          I
      DO 10 I = 1,N
          TEMP = Y(I)
          Y(I) = X(I) - Y(I)
          X(I) = TEMP
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVSBP                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR X(N) IS PERMUTED ACCORDING TO THE FORMULA
* X(PERM(I)):=X(I).
*
* PARAMETERS :
*  II  N  LENGTH OF VECTORS.
*  II  PERM(N)  INPUT PERMUTATION VECTOR.
*  RU  X(N)  VECTOR THAT IS TO BE PERMUTED.
*  RA  RN01(N)  AUXILIARY VECTOR.
*
      SUBROUTINE MXVSBP(N,PERM,X,RN01)
      INTEGER N,PERM(*),I
      DOUBLE PRECISION RN01(*),X(*)
C
      DO 100 I=1,N
      RN01(PERM(I))=X(I)
100   CONTINUE
      DO 200 I=1,N
      X(I)=RN01(I)
200   CONTINUE
      RETURN
      END
* SUBROUTINE MXVSCL                ALL SYSTEMS                88/12/01
C PORTABILITY : ALL SYSTEMS
C 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALING OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  A  SCALING FACTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= A*X.
*
      SUBROUTINE MXVSCL(N,A,X,Y)
      INTEGER N
      DOUBLE PRECISION  A, X(*), Y(*)
      INTEGER I
      DO 1 I=1,N
      Y(I) = A*X(I)
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXVSET                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  A  INITIAL VALUE.
*  RO  X(N)  OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
*
      SUBROUTINE MXVSET(N,A,X)
      DOUBLE PRECISION A
      INTEGER          N
      DOUBLE PRECISION X(*)
      INTEGER          I
      DO 10 I = 1,N
          X(I) = A
   10 CONTINUE
      RETURN
      END
* SUBROUTINE MXVSFP                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR X(N) IS PERMUTED ACCORDING TO THE FORMULA
* X(I)=X(PERM(I)).
*
* PARAMETERS :
*  II  N  LENGTH OF VECTORS.
*  II  PERM(N)  INPUT PERMUTATION VECTOR.
*  RU  X(N)  VECTOR THAT IS TO BE PERMUTED.
*  RA  RN01(N)  AUXILIARY VECTOR.
*
      SUBROUTINE MXVSFP(N,PERM,X,RN01)
      INTEGER N,PERM(*),I
      DOUBLE PRECISION RN01(*),X(*)
C
      DO 100 I=1,N
      RN01(I)=X(PERM(I))
100   CONTINUE
      DO 200 I=1,N
      X(I)=RN01(I)
200   CONTINUE
      RETURN
      END
* SUBROUTINE MXVSIP                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* THE VECTOR OF THE INVERSE PERMUTATION IS COMPUTED.
*
* PARAMETERS :
*  II  N  LENGTH OF VECTORS.
*  II  PERM(N)  INPUT PERMUTATION VECTOR.
*  IO  INVP(N)  INVERSE PERMUTATION VECTOR.
*
      SUBROUTINE MXVSIP(N,PERM,INVP)
      INTEGER N,PERM(*),INVP(*)
      INTEGER I,J
CC
      DO 100 I=1,N
      J=PERM(I)
      INVP(J)=I
100   CONTINUE
      RETURN
      END
* SUBROUTINE MXVSR2                ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* RADIXSORT.
*
* PARAMETERS :
*  II  MCOLS  NUMBER OF INTEGER VALUES OF THE SORTED ARRAY.
*  RI  DEG(MCOLS)  VALUES OF THE SORTED ARRAY.
*  RO  ORD(MCOLS)  SORTED OUTPUT.
*  RA  RADIX(MCOLS+1)  AUXILIARY ARRAY.
*  II  WN01(MCOLS)  INDICES OF THE SORTED ARRAY.
*  II  LENGTH   NUMBER OF SORTED PIECES.
*
      SUBROUTINE MXVSR2(MCOLS,DEG,ORD,RADIX,WN01,LENGTH)
      INTEGER MCOLS,WN01(*)
      DOUBLE PRECISION DEG(*),ORD(*),RADIX(*)
      INTEGER LENGTH,I,L,L1,L2
C
C RADIX IS SHIFTED : 0-(MCOLS-1) --- 1-MCOLS
C
      DO 50 I=1,MCOLS+1
        RADIX(I)=0
50    CONTINUE
      DO 100 I=1,LENGTH
        L2=WN01(I)
        L=DEG(L2)
        RADIX(L+1)=RADIX(L+1)+1
100   CONTINUE
C
C     RADIX COUNTS THE NUMBER OF VERTICES WITH DEG(I)>=L
C
      L=0
      DO 200 I=MCOLS,0,-1
        L=RADIX(I+1)+L
        RADIX(I+1)=L
200   CONTINUE
C
C     ARRAY ORD IS FILLED
C
      DO 300 I=1,LENGTH
        L2=WN01(I)
        L=DEG(L2)
        L1=RADIX(L+1)
        ORD(L1)=L2
        RADIX(L+1)=L1-1
300   CONTINUE
      RETURN
      END
* SUBROUTINE MXVSR5                ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 TU : ORIGINAL VERSION
*
* PURPOSE :
* SHELLSORT.
*
* PARAMETERS :
*  II  K  NUMBER OF INTEGER VALUES OF THE SORTED ARRAY.
*  II  L  CORRECTION FOR THE ABSOLUTE INDEX IN THE SORTED ARRAY
*  IU  ARRAY(K)  INTEGER SORTED ARRAY.
*  RO  ARRAYC(K)  REAL OUTPUT ARRAY.
*  RU  ARRAYD(K)  REAL ARRAY WHICH IS PERMUTED IN THE SAME WAY
*         AS THE INTEGER SORTED ARRAY.
*
      SUBROUTINE MXVSR5(K,L,ARRAY,ARRAYC,ARRAYD)
      INTEGER K,L
      INTEGER ARRAY(*)
      DOUBLE PRECISION ARRAYC(*),ARRAYD(*)
      INTEGER IS,LA,LT,LS,LLS,I,J,JS,KHALF
      DOUBLE PRECISION LD
C
C NOTHING TO BE SORTED
C
      IF(K.LE.1) GO TO 400
C
C SHELLSORT
C
C L - CORRECTION FOR THE ABSOLUTE INDEX IN THE SORTED ARRAY
C
      LS=131071
      KHALF=K/2
      DO 300 LT=1,17
        IF(LS.GT.KHALF) THEN
          LS=LS/2
          GO TO 300
        END IF
        LLS=K-LS
        DO 200 I=1,LLS
          IS=I+LS
          LA=ARRAY(IS)
          LD=ARRAYD(IS)
          J=I
          JS=IS
 100      IF(LA.GE.ARRAY(J)) THEN
            ARRAY(JS)=LA
            ARRAYD(JS)=LD
            ARRAYC(INT(LD))=JS+L
            GO TO 200
          ELSE
            ARRAY(JS)=ARRAY(J)
            ARRAYD(JS)=ARRAYD(J)
            ARRAYC(INT(ARRAYD(J)))=JS+L
            JS=J
            J=J-LS
          END IF
          IF(J.GE.1) GO TO 100
          ARRAY(JS)=LA
          ARRAYD(JS)=LD
          ARRAYC(INT(LD))=JS+L
 200    CONTINUE
        LS=LS/2
 300  CONTINUE
 400  CONTINUE
      RETURN
      END
* SUBROUTINE MXVSUM                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SUM OF TWO VECTORS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X + Y.
*
      SUBROUTINE MXVSUM(N,X,Y,Z)
      INTEGER          N
      DOUBLE PRECISION X(*),Y(*),Z(*)
      INTEGER          I
      DO 10 I = 1,N
          Z(I) = X(I) + Y(I)
   10 CONTINUE
      RETURN
      END
* FUNCTION MXWDOT                  ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DOT PRODUCT OF TWO VECTORS IN THE PACKED CASE.
*
* PARAMETERS :
*  II  L  PACKED OR UNPACKED VECTOR DIMENSION.
*  II  N  UNPACKED VECTOR DIMENSION.
*  II  JBL(L)  INDICES OF PACKED VECTOR.
*  RI  X(L)  UNPACKED OR PACKED INPUT VECTOR.
*  RI  Y(N)  UNPACKED INPUT VECTOR.
*  II  JOB  FORM OF THE VECTOR X. JOB=1-UNPACKED FORM. JOB=2-PACKED
*          FORM.
*  RR  MXWDOT  VALUE OF DOT PRODUCT MXWDOT=TRANS(X)*Y.
*
      FUNCTION MXWDOT(L,JBL,X,Y,JOB)
      INTEGER L,JBL(*),JOB
      DOUBLE PRECISION  X(*), Y(*), MXWDOT
      DOUBLE PRECISION TEMP
      INTEGER I,IP
      TEMP=0.0D0
      IF (JOB.EQ.1) THEN
      DO 1 I=1,L
      IP=JBL(I)
      IF (IP.GT.0) TEMP=TEMP+X(IP)*Y(IP)
    1 CONTINUE
      ELSE
      DO 2 I=1,L
      IP=JBL(I)
      IF (IP.GT.0) TEMP=TEMP+X(I)*Y(IP)
    2 CONTINUE
      END IF
      MXWDOT=TEMP
      RETURN
      END