C     ALGORITHM 408 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 04,
C     P. 265.
      FUNCTION IND(M,I,NM)
C     ********************
C UNPACKS I*TH COLUMN INDEX.ARRAY M CONTAINS TWO 4-DIGIT
C INDICES PER WORD, LOWER INDEX IN UPPER FOUR DIGITS.
C
      DIMENSION M(NM)
C J*TH WORD OF M CONTAINS I*TH INDEX.
      J      = (I+1)/2
C L IS 0 IF I EVEN, 1 IF I ODD.
      L      = I-(I/2)*2
C KT CONTAINS UPPER 4 DIGITS OF M(J).
      KT     = M(J)/10000
      IF (L) 1,1,2
    1 IND    = M(J)-KT*10000
      RETURN
    2 IND    = KT
      RETURN
      END
      SUBROUTINE IPK(K,M,I,NM)
C     ************************
C PACKS K (I*TH COLUMN INDEX) IN ARRAY M, WHICH WILL
C CONTAIN TWO 4-DIGIT INDICES PER WORD, LOWER INDEX
C UPPER 4 DIGITS.
C
      DIMENSION M(NM)
      J      = (I+1)/2
      L      = I-(I/2)*2
      IF (L) 1,1,2
    1 M(J)   = M(J)+K
      RETURN
    2 M(J)   = M(J)+K*10000
      RETURN
      END
C **********************************************************
      SUBROUTINE ACSPMX(A,M,AN,MN,R,IR,IT,NA,NM)
C     ******************************************
C
C ADDS R TIMES COL IR TO COL IT OF MATRIX STORED IN A,
C PLACING RESULT IN AN.  M,MN CONTAIN CONTROL DATA AND COL
C INDICES FOR A,AN.
C NA IS DIMENSION OF A,AN. NM IS DIMENSION OF M,MN.
      REAL A,AN,R,AR
      INTEGER M,MN,IR,IT,NA,NM,I,NRA,NCA,L,JF,NIR,NIRA,J2,K,
     * K1,IFL,LP,J
      DIMENSION A(NA),M(NM),AN(NA),MN(NM)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CHECK THAT PARAMETERS WITHIN RANGE.
      IF (R.EQ.0.0) WRITE (LP,14)
      IF (IR.LE.0.OR.IT.LE.0) GO TO 15
C CLEAR MN.
      DO 1 I = 1,NM
    1  MN(I) = 0
C CHECK THAT AN DOES NOT OVER-WRITE A.
      IF (M(1).EQ.0) GO TO 17
C UNPACK AND TRANSFER CONTROL DATA. NRA,NCA ARE NUMBERS OF
C ROWS,COLS IN A.
      NRA    = IND(M,1,NM)
      NCA    = IND(M,2,NM)
      MN(1)  = M(1)
C CHECK PARAMETERS WITHIN RANGE.
      IF (IR.GT.NCA.OR.IT.GT.NCA) GO TO 15
C L COUNTS ELEMENTS OF AN.
      L      = 1
      JF     = 4+NRA
C J COUNTS ELEMENTS TO END OF ROW (I-1).
      J      = 0
C I COUNTS ROWS OF A.
      DO 13 I  = 1,NRA
C NIR IS NUMBER IN NEW ROW.
       NIR     = 0
C NIRA IS NUMBER IN ROW I OF A.
       NIRA    = IND(M,4+I,NM)
       IF (NIRA.EQ.0) GO TO 12
C J2 COUNTS ELEMENTS TO END OF CURRENT ROW.
       J2      = J+NIRA
C J1 COUNTS ELEMENTS UP TO FIRST ONE IN CURRENT ROW.
       J1      = J+1
C PICK OUT ELEMENT IN COLUMN IR.
       AR      = 0.
       DO 2 K  = J1,J2
        K1     = IND(M,K+JF,NM)
        IF (K1.NE.IR) GO TO 2
        AR     = A(K)
    2  CONTINUE
C PICK OUT AND ALTER IF NECESSARY ELEMENT IN COL IT.
C TRANSFER REST OF ROW TO AN,MN. IFL SET TO 1 WHEN ELEMENT
C IN COL IT FOUND OR CREATED.
       IFL     = 0
       K       = J1-1
    3  K      = K+1
       K1     = IND(M,K+JF,NM)
       A1     = A(K)
       IF (K1.GE.IT) GO TO 7
C CHECK IF ARRAYS FULL.
    4  IF (L.LE.NA.AND.JF+L.LE.2*NM) GO TO 6
       WRITE (LP,5)
    5  FORMAT(21H IN ACSPMX ARRAY FULL)
       CALL EXIT
C COLUMN IT NOT YET REACHED.
    6  AN(L)  = A1
       CALL IPK(K1,MN,JF+L,NM)
       L      = L+1
       NIR    = NIR+1
       GO TO 10
    7  IF (K1.GT.IT) GO TO 8
C K1 EQUALS IT, I.E. THERE IS A NON-ZERO ELEMENT IN COL. IT
C OF ROW I.
       IFL    = 1
       A1     = AR*R+A(K)
       IF (A1.NE.0.0) GO TO 4
       GO TO 10
    8  IF (IFL.NE.0) GO TO 9
C K1 GREATER THAN IT AND ELEMENT IN COL IT HAS NOT YET
C BEEN FOUND, THUS COL IT HAS A ZERO ELEMENT.
       IFL    = 1
       K      = K-1
       A1     = AR*R
C A NEW ELEMENT IN COL IT IS CREATED IF A1 NOT ZERO.
       K1     = IT
       IF (A1.NE.0.0) GO TO 4
       GO TO 10
C K GREATER THAN IT AND ELEMENT IN COL IT ALREADY FOUND OR
C CREATED, JUST TRANSFER TO NEW ARRAY.
    9  A1     = A(K)
       GO TO 4
   10  IF (K.LT.J2) GO TO 3
       IF (IFL.NE.0.OR.AR.EQ.0.0) GO TO 12
       IF (L.LE.NA.AND.JF+L.LE.2*NM) GO TO 11
       WRITE (LP,5)
       CALL EXIT
   11  AN(L) = AR*R
       CALL IPK(IT,MN,JF+L,NM)
       L     = L+1
       NIR   = NIR+1
C END OF ROW.
   12  CALL IPK(NIR,MN,4+I,NM)
       J     = J2
   13 CONTINUE
C END OF LAST ROW.
      MN(2)  = L-1
      RETURN
C ERROR MESSAGES.
   14 FORMAT(20H IN ACSPMX R IS ZERO)
   15 WRITE (LP,16)
   16 FORMAT(32H IN ACSPMX IR OR IT OUT OF RANGE)
      CALL EXIT
   17 WRITE (LP,18)
   18 FORMAT(35H IN ACSPMX OUTPUT OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE ADSPMX(A,MA,B,MB,C,MC,NA,NM)
C     ***************************************
C
C ADD TWO SPARSE MATRICES.
C
C A,B,C CONTAIN ELEMENTS OF FIRST,SECOND AND SUM MATRICES.
C MA,MB,MC, CONTAIN CONTROL DATA AND COL. INDICES FOR A,B,C.
C NA IS DIMENSION OF A,B,C. NM IS DIMENSION OF MA,MB,MC.
C
      REAL A,B,C
      INTEGER MA,MB,MC,LP,NRA,NCA,NRB,NCB,JC,KA,KB,JB,KF,
     * I,KA1,JA,J1,J2,NOLD,NA,NM
      DIMENSION A(NA),MA(NM),B(NA),MB(NM),C(NA),MC(NM)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CLEAR MC.
      DO 1 I = 1,NM
    1  MC(I) = 0
C CHECK THAT C DOES NOT OVER-WRITE A OR B.
      IF (MA(1).EQ.0.OR.MB(1).EQ.0) GO TO 18
C UNPACK CONTROL DATA.  NRA,NCA ARE NUMBER OF ROWS,COLUMNS
C IN A.  NRB, NCB ARE ROWS, COLUMNS IN B.
      NRA    = IND(MA,1,NM)
      NCA    = IND(MA,2,NM)
      NRB    = IND(MB,1,NM)
      NCB    = IND(MB,2,NM)
C TEST FOR COMPATIBILITY.
      IF (NRA.EQ.NRB) GO TO 3
      WRITE (LP,2) NRA,NRB
    2 FORMAT(32H IN ADSPMX NUMBER OF ROWS IN A (,I4,
     * 37H) DOES NOT EQUAL NUMBER OF ROWS IN B(,I4,2H).)
      CALL EXIT
    3 IF (NCA.EQ.NCB) GO TO 5
      WRITE (LP,4) NCA,NCB
    4 FORMAT(31H IN ADSPMX NUMBER OF COLS IN A(,I4,
     * 36H) DOES NOT EQUAL NUM. OF COLS. IN B(,I4,2H).)
      CALL EXIT
C JC COUNTS ELEMENTS OF C.
    5 JC     = 1
C KA,KB ARE NUMBERS IN FIRST I ROWS OF A,B.
      KA     = 0
      KB     = 0
C KF IS NUMBER OF CONTROL DATA IN A,B OR C.
      KF     = 4+NRA
C JB COUNTS ELEMENTS OF B.
      JB     = 1
C I COUNTS ROWS OF A,B,C.
      DO 15 I=1,NRA
       KB      = KB+IND(MB,4+I,NM)
C NIRA IS NUMBER IN ROW I OF A.
       NIRA    = IND(MA,4+I,NM)
       IF (NIRA.EQ.0) GO TO 12
       KA1     = KA+1
       KA      = KA+NIRA
C JA COUNTS ELEMENTS OF A.
       DO 11 JA= KA1,KA
    6   J1     = IND(MA,JA+KF,NM)
C AT END OF B-ROW TRANSFER REST OF A-ROW.
        IF (JB.GT.KB) GO TO 7
        J2     = IND(MB,JB+KF,NM)
        IF (J1-J2) 7,9,10
C IF A-INDEX LESS THAN B-INDEX TRANSFER A-ELEMENT TO C.
    7   IF (JC.GT.NA) GO TO 16
        C(JC)  = A(JA)
    8   IF (JC+KF.GT.2*NM) GO TO 16
        CALL IPK(J1,MC,JC+KF,NM)
        JC     = JC+1
        GO TO 11
C IF A-INDEX EQUALS B-INDEX ADD ELEMENTS ,PLACE SUM IN C.
    9   T      = A(JA)+B(JB)
C IGNORE SUM ELEMENT IF ZERO.
        IF (T.EQ.0.0) GO TO 11
        IF (JC.GT.NA) GO TO 16
        C(JC)  = T
        JB     = JB+1
        GO TO 8
C IF A-INDEX GREATER THAN B-INDEX TRANSFER B-ELEMENT TO C.
   10   IF (JC.GT.NA) GO TO 16
        C(JC)  = B(JB)
        IF (JC+KF.GT.2*NM) GO TO 16
        CALL IPK(J2,MC,JC+KF,NM)
        JB     = JB+1
        JC     = JC+1
        GO TO 6
   11  CONTINUE
C END OF ROW OF A.  TRANSFER REST OF ROW OF B.
   12  IF (JB.GT.KB) GO TO 13
       IF (JC.GT.NA) GO TO 16
       C(JC) = B(JB)
       J2    = IND(MB,JB+KF,NM)
       IF (JC+KF.GT.2*NM) GO TO 16
       CALL IPK(J2,MC,JC+KF,NM)
       JC    = JC+1
       JB      = JB+1
       GO TO 12
   13  IF (I.GT.1) GO TO 14
       NOLD  = JC-1
C NIRC IS NUMBER IN ROW I OF C.
       NIRC  = JC-1
       GO TO 15
   14  NIRC   = JC-1-NOLD
       NOLD  = JC-1
   15  CALL IPK(NIRC,MC,4+I,NM)
C
C LAST ROW.   STORE CONTROL DATA IN MC.
      CALL IPK(NRA,MC,1,NM)
      CALL IPK(NCA,MC,2,NM)
      MC(2)  = JC-1
      RETURN
C ERROR MESSAGES.
   16 WRITE (LP,17)
   17 FORMAT(41H IN ADSPMX SPACE FOR SUM MATRIX EXCEEDED.)
      CALL EXIT
   18 WRITE (LP,19)
   19 FORMAT(33H IN ADSPMX SUM OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE ANTIP(N,P,AP)
C     ************************
C
C INVERT PERMUTATION IN P, PLACING RESULT IN AP.
C N IS NUMBER OF ELEMENTS IN P,AP.
C
      INTEGER P(N),AP(N)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CHECK THAT OUTPUT DOES NOT OVER-WRITE INPUT.
      AP(1)  = 0
      IF (P(1).EQ.0) GO TO 4
      DO 3 I = 1,N
       J     = P(I)
       IF (J) 1,1,3
    1  WRITE (LP,2) I
    2  FORMAT(37H IN ANTIP PERMUTATION CONTAINS A NON-,
     * 28HPOSITIVE NUMBER IN POSITION ,I5)
       CALL EXIT
    3  AP(J) = I
      RETURN
C ERROR MESSAGES.
    4 WRITE (LP,5)
    5 FORMAT(42H IN ANTIP OUTPUT OVER-WRITES INPUT OR P(1),
     * 8H IS ZERO)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE ARSPMX(A,M,AN,MN,R,IR,IT,NA,NM)
C     ******************************************
C
C ADD R TIMES ROW IR OF SPARSE MATRIX TO ROW IT.
C
C A,M CONTAIN ELEMENTS, COLUMN INDICES OF INPUT MATRIX.
C AN,MN CONTAIN ELEMENTS, COLUMN INDICES OF NEW MATRIX.
C NA IS DIMENSION OF A,AN.  NM IS DIMENSION OF M,MN.
C
      REAL A,AN,R
      INTEGER M,MN,IR,IT,I,NRA,NCA,NEA,NIRA,JR,I1,JT,J,JF,
     * J1,JT2,JN,J2,KR,KT,JT1,IT1,K,NA,NM,JT0,LP
      DIMENSION A(NA),M(NM),AN(NA),MN(NM)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CHECK PARAMETERS WITHIN RANGE.
      IF (IR.LE.0.OR.IT.LE.0) GO TO 23
C CLEAR MN.
      DO 1 I = 1,NM
    1  MN(I) = 0
C CHECK THAT AN DOES NOT OVER-WRITE A.
      IF (M(1).EQ.0) GO TO 25
C UNPACK CONTROL DATA, STORE IN MN.
      NRA    = IND(M,1,NM)
      NEA    = M(2)
      MN(1)  = M(1)
      DO 2 I = 1,NRA
       IF (I.EQ.IT) GO TO 2
       K     = IND(M,4+I,NM)
       CALL IPK(K,MN,4+I,NM)
    2  CONTINUE
C CHECK PARAMETERS WITHIN RANGE.
      IF (IR.GT.NRA.OR.IT.GT.NRA) GO TO 23
C JR,JT ARE NUMBERS OF ELEMENTS BELOW ROWS IR,IT.
      JR     = 0
      IF (IR.LE.1) GO TO 4
      IR1    = IR-1
      DO 3 I = 1,IR1
    3  JR    = JR+IND(M,4+I,NM)
    4 JT     = 0
      JF     = 4+NRA
      IF (IT.LE.1) GO TO 7
      IT1    = IT-1
      DO 5 I = 1,IT1
    5  JT    = JT+IND(M,4+I,NM)
C TRANSFER ELEMENTS BELOW ROW IT.
      DO 6 I = 1,JT
       AN(I) = A(I)
       J     = IND(M,JF+I,NM)
    6  CALL IPK(J,MN,JF+I,NM)
    7 JT0    = JT
C ADD R TIMES ROW IR TO ROW IT.
      JT2    = JT+IND(M,4+IT,NM)
      JT     = JT+1
C JN COUNTS ELEMENTS OF NEW MATRIX.
      JN     = JT
C NIRR IS NUMBER OF ELEMENTS IN ROW IR OF A.
      NIRR   = IND(M,4+IR,NM)
      IF (NIRR.EQ.0) GO TO 14
      J1     = JR+1
      J2     = JR+NIRR
      DO 13 I= J1,J2
C CHECK ARRAY LIMIT.
       IF (JN.LE.NA.AND.(JN+JF).LE.2*NM) GO TO 8
       WRITE (LP,21)
       CALL EXIT
    8  KR    = IND(M,JF+I,NM)
    9  IF (JT.GT.JT2) GO TO 12
       KT    = IND(M,JF+JT,NM)
       IF (KT.GE.KR) GO TO 10
       AN(JN)=A(JT)
       CALL IPK(KT,MN,JN+JF,NM)
       JN    = JN+1
       JT    = JT+1
       GO TO 9
   10  IF (KT.GT.KR) GO TO 12
       S     = A(JT)+R*A(I)
       IF (S.EQ.0.0) GO TO 11
       AN(JN)=S
       CALL IPK(KT,MN,JN+JF,NM)
       JN    = JN+1
   11  JT    = JT+1
       GO TO 13
   12  S     = R*A(I)
       IF (S.EQ.0.0) GO TO 13
       AN(JN)=S
       CALL IPK(KR,MN,JN+JF,NM)
       JN    = JN+1
   13 CONTINUE
C TRANSFER REST OF ROW IT.
   14 IF (JT.GT.JT2) GO TO 16
      IF (JN.LE.NA.AND.(JN+JF).LE.2*NM) GO TO 15
      WRITE (LP,21)
      CALL EXIT
   15 AN(JN) = A(JT)
      KT     = IND(M,JF+JT,NM)
      CALL IPK(KT,MN,JF+JN,NM)
      JN     = JN+1
      JT     = JT+1
      GO TO 14
C JT1 IS NUMBER IN ROW IT OF AN.
   16 JT1    = JN-1-JT0
      CALL IPK(JT1,MN,4+IT,NM)
C STORE ROWS ABOVE IT.
      IF (IT.EQ.NRA) GO TO 20
      IT1    = IT+1
      DO 19 I=IT1,NRA
C K IS NUMBER IN ROW I.
       K     = IND(M,4+I,NM)
       IF (K.EQ.0) GO TO 19
       JT    = JT2
       JT2   = JT+K
       JT1   = JT+1
       DO 18 J = JT1,JT2
        IF (JN.LE.NA.AND.(JN+JF).LE.2*NM) GO TO 17
        WRITE (LP,21)
        CALL EXIT
   17   AN(JN) = A(J)
        K      = IND(M,J+JF,NM)
        CALL IPK(K,MN,JN+JF,NM)
   18   JN     = JN+1
   19 CONTINUE
   20 MN(2)  = JN-1
      RETURN
C ERROR MESSAGES.
   21 FORMAT(21H0ARRAY FULL IN ARSPMX)
   22 FORMAT(20H IN ARSPMX R IS ZERO)
   23 WRITE (LP,24)
   24 FORMAT(33H IN ARSPMX IR OR IT OUT OF RANGE.)
      CALL EXIT
   25 WRITE (LP,26)
   26 FORMAT(35H IN ARSPMX OUTPUT OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE CVSPMX(A,M,IC,V,AT,MT,N,NA,NM,IP)
C     ********************************************
C EXTRACTS COL IC OF SPARSE MATRIX IN A,STORING RESULT IN V
C IN EXTENDED FORM, I.E. ALL ELEMENTS INCLUDING ZEROS ARE
C REPRESENTED.
C M CONTAINS COLUMN INDICES OF A.
C AT,MT ARE USED INTERNALLY.
C N IS DIMENSION OF V. NA IS DIMENSION OF A,AT.
C NM IS DIMENSION OF M,MT.
C IP US USED INTERNALLY BY TRSPMX.
C
      REAL A,AT,V
      INTEGER M,IC,MT,N,NA,NM,IP
      DIMENSION A(NA),M(NM),AT(NA),MT(NM),V(N),IP(N)
      CALL TRSPMX(A,M,AT,MT,NA,NM,IP,N)
      CALL RVSPMX(AT,MT,IC,V,N,NA,NM)
      RETURN
      END
C **********************************************************
      SUBROUTINE ECSPMX(A,M,AN,MN,IR,JR,J,AT,MT,NA,NM,NP)
C     ***************************************************
C
C RESULT IN AN.
C EXCHANGE COLUMNS IR,JR OF SPARSE MATRIX IN A, STORING
C M,MN CONTAIN COLUMN INDICES OF A,AN. J IS USED INTERNALLY.
C AT,MT ARE USED BY PERCOL.
C NA IS DIMENSION OF A,AN. NM IS DIMENSION OF M,MN.
C NP IS DIMENSION OF J,AT,MT.
C
      REAL A,AN,AT
      INTEGER M,MN,IR,JR,J,NCA,NA,NM,NP,MT
      DIMENSION A(NA),M(NM),AN(NA),
     * MN(NM),J(NP),AT(NP),MT(NP)
C SET UP PERMUTATION ARRAY J WITH IR,JR INTERCHANGED.
C NCA IS NUMBER OF COLUMNS IN A.
      NCA    = IND(M,2,NM)
      DO 1 I = 1,NCA
    1  J(I)  = I
      J(IR)  = JR
      J(JR)  = IR
C PERMUTE COLS OF A.
      CALL PERCOL(A,M,AN,MN,J,AT,MT,NA,NM,NP)
      RETURN
      END
C **********************************************************
      SUBROUTINE ERSPMX(A,M,AN,MN,IR,JR,J,NA,NM,NP)
C     *********************************************
C
C EXCHANGE ROWS IR,JR OF SPARSE MATRIX IN A,STORING RESULT
C IN AN. M, MN CONTAIN COLUMN INDICES OF A, AN.
C NA IS DIMENSION OF A,AN.  NM IS DIMENSION OF M,MN.
C NP IS DIMENSION OF J,WHICH IS USED INTERNALLY.
C
      REAL A,AN
      INTEGER M,MN,IR,JR,NRA,NA,NM,J,NP
      DIMENSION A(NA),M(NM),AN(NA),MN(NM),J(NP)
C SET UP PERMUTATION ARRAY WITH IR,JR INTERCHANGED AND ALL
C  OTHER INTEGERS IN NATURAL ORDER.
C NRA IS NUMBER OF ROWS IN A.
      NRA    = IND(M,1,NM)
      DO 1 I = 1,NRA
    1  J(I)  = I
      J(IR)  = JR
      J(JR)  = IR
C PERMUTE ROWS OF A.
      CALL PERROW(A,M,AN,MN,J,NA,NM,NP)
      RETURN
C
      END
C **********************************************************
      SUBROUTINE OTSPMX(A,M,N,NA,NM)
C     ******************************
C WRITE SPARSE MATRIX IN A,M ON FORTRAN UNIT N (MASS STORAGE
C DEVICE).
C NA,NM ARE DIMENSIONS OF A,M.
C
      REAL A
      INTEGER M,N,NEA,NRA,NEM,NA,NM
      DIMENSION A(NA),M(NM)
C NEA IS NUMBER OF ELEMENTS IN A.
      NEA    = M(2)
C NRA IS NUMBER OF ROWS IN A.
      NRA    = IND(M,1,NM)
C NEM IS NUMBER OF WORDS IN M.
      NEM    = (5+NRA+NEA)/2
      WRITE (N) NEM
      WRITE (N) (M(I),I=1,NEM)
      WRITE (N) (A(I),I=1,NEA)
      REWIND N
      RETURN
      END
C **********************************************************
      SUBROUTINE INSPMX(A,M,N,NA,NM)
C     *******************************
C READ SPARSE MATRIX FROM FORTRAN UNIT N (MASS STORAGE
C DEVICE).  STORE IN A, WITH COLUMN INDEX ARRAY IN M.
C NA,NM ARE DIMENSIONS OF A,M.
C
      REAL A
      INTEGER M,N,NEM,NEA
      DIMENSION A(NA),M(NM)
C
C NEM IS NUMBER OF WORDS IN M.
      READ (N) NEM
      READ (N) (M(I),I=1,NEM)
C NEA IS NUMBER OF ELEMENTS IN A.
      NEA    = M(2)
      READ (N) (A(I),I=1,NEA)
      REWIND N
      RETURN
      END
C     ******************************************************
      SUBROUTINE MCSPMX(A,M,AN,MN,R,IC,NA,NM)
C     ***************************************
C
C MULTIPLIES COL IC OF SPARSE MATRIX IN A BY R, STORING
C RESULT IN AN.  M,MN CONTAIN COLUMN INDICES OF A,AN.
C NA IS DIMENSION OF A,AN. NM IS DIMENSION OF M,MN.
C
      INTEGER M,MN,IC,NRA,NCA,NEA,I1,I,J,JF,NIRA,J2,J1,K,IT,
     * NA,NM,L,NC
      DIMENSION A(NA),M(NM),AN(NA),MN(NM)
      REAL A,AN,R
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CLEAR MN.
      DO 1 I = 1,NM
    1  MN(I) = 0
C CHECK THAT OUTPUT DOES NOT OVER-WRITE INPUT.
      IF (M(1).EQ.0) GO TO 9
C UNPACK AND TRANSFER CONTROL DATA. NRA,NCA,NEA ARE NUMBERS
C  OF ROWS, COLS, ELEMENTS IN A.
      NRA    = IND(M,1,NM)
      NCA    = IND(M,2,NM)
      NEA    = M(2)
      DO 2 I = 1,2
    2  MN(I) = M(I)
C CHECK IC WITHIN RANGE.
      IF (IC.GT.NCA.OR.IC.LE.0) GO TO 7
C J COUNTS ELEMENTS TO END OF ROW (I-1) OF A.
      J      = 0
      JF     = 4+NRA
C L COUNTS ELEMENTS OF NEW MATRIX.
      L      = 1
C I COUNTS ROW OF A.
      DO 5 I= 1,NRA
C NIRA IS NUMBER OF ELEMENTS IN ROW I OF A.
       NIRA  = IND(M,4+I,NM)
C NIRAN IS NUMBER OF ELEMENTS IN ROW OF AN.
      NIRAN = NIRA
       IF (NIRA.EQ.0) GO TO 5
C J2 COUNTS ELEMENTS TO END OF ROW I OF A.
       J2    = J+NIRA
C J1 COUNTS ELEMENTS UP TO FIRST ONE IN ROW I OF A.
       J1    = J+1
C PROCESS ROW I OF A.
       DO 4 K  = J1,J2
        IT     = IND(M,JF+K,NM)
        IF (IT.EQ.IC) GO TO 3
C TRANSFER COLUMNS OTHER THAN IC.
        AN(L)  = A(K)
        I1     = IND(M,JF+K,NM)
        CALL IPK(I1,MN,JF+L,NM)
        L      = L+1
        GO TO 4
C MULTIPLY COL IC BY R.
    3   IF (R.EQ.0.0) NIRAN= NIRA-1
        IF (R.EQ.0.0) GO TO 4
        AN(L)  = R*A(K)
        I1     = IND(M,JF+K,NM)
        CALL IPK(I1,MN,JF+L,NM)
        L      = L+1
    4  CONTINUE
C END OF ROW I.
       J     = J2
       CALL IPK(NIRAN,MN,4+I,NM)
    5 CONTINUE
C END OF LAST ROW.
      IF (R.NE.0.0) GO TO 6
      M(2)   = L-1
    6 RETURN
C ERROR MESSAGES.
    7 WRITE (LP,8)
    8 FORMAT(26H IN MCSPMX IC OUT OF RANGE)
      CALL EXIT
    9 WRITE (LP,10)
   10 FORMAT(35H IN MCSPMX OUTPUT OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C     ******************************************************
      SUBROUTINE MRSPMX(A,M,AN,MN,R,IR,NA,NM)
C     ***************************************
C
C MULTIPLIES ROW IR OF SPARSE MATRIX A BY R, STORING
C RESULT IN AN.  M,MN CONTAIN COLUMN INDICES OF A,AN.
C NA IS DIMENSION OF A,AN.  NM IS DIMENSION OF M,MN.
      REAL A,AN,R
      INTEGER M,MN,IR,NM,NRA,NEA,I1,I,
     * J,NIRA,J2,J1,K,NA,L,JF
      DIMENSION A(NA),M(NM),AN(NA),MN(NM)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CLEAR MN.
      DO 1 I = 1,NM
    1  MN(I) =0
C CHECK THAT OUTPUT DOES NOT OVER-WRITE INPUT.
      IF (M(1).EQ.0) GO TO 9
C UNPACK AND TRANSFER CONTROL DATA.
      NRA    = IND(M,1,NM)
      NEA    = M(2)
      DO 2 I = 1,2
    2  MN(I) = M(I)
C CHECK THAT IR IS WITHIN RANGE.
      IF (IR.GT.NRA.OR.IR.LT.1) GO TO 11
C J COUNTS ELEMENTS TO END OF ROW (I-1) OF A.
      J      = 0
C L COUNTS ELEMENTS OF NEW MATRIX.
      L      = 1
      JF     = 4+NRA
C I COUNTS ROWS OF A.
      DO 7 I = 1,NRA
C NIRA IS NUMBER IN ROW I OF A.
       NIRA  = IND(M,4+I,NM)
       CALL IPK(NIRA,MN,4+I,NM)
       IF (NIRA.EQ.0) GO TO 7
C J2 COUNTS ELEMENTS TO END OF ROW I OF A.
       J2    = J+NIRA
C J1 COUNTS ELEMENTS TO FIRST ONE IN ROW I OF A.
       J1    = J+1
       IF (I.EQ.IR) GO TO 4
C TRANSFER ROWS OTHER THAN I.
       DO 3 K  = J1,J2
        I1     = IND(M,JF+K,NM)
        CALL IPK(I1,MN,JF+L,NM)
        AN(L)  = A(K)
    3   L      = L+1
       GO TO 6
C MULTIPLY ROW IR BY R.
    4  DO 5 K  = J1,J2
        IF (R.EQ.0.0) GO TO 5
        I1     = IND(M,JF+K,NM)
        CALL IPK(I1,MN,JF+L,NM)
        AN(L)  = R*A(K)
        L      = L+1
    5   CONTINUE
    6  J       = J2
C END OF ROW I.
    7 CONTINUE
C END OF LAST ROW.
      IF (R.NE.0.0) GO TO 8
      M(2)   = L-1
      CALL IPK(0,MN,4+IR,NM)
    8 RETURN
C ERROR MESSAGES.
    9 WRITE (LP,10)
   10 FORMAT(36H IN MRSPMX OUTPUT OVER-WRITES INPUT ,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
   11 WRITE (LP,12)
   12 FORMAT(26H IN MRSPMX IR OUT OF RANGE)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE MUSPMX(A,MA,B,MB,C,MC,NA,NM)
C     ***************************************
C
C MULTIPLY TWO SPARSE MATRICES.
C
C B MUST BE STORED BY COLUMNS, I.E. WE FORM C =
C A*(B TRANSPOSED).
C A,B,C CONTAIN FIRST,SECOND AND PRODUCT MATRICES RESPECT-
C IVELY.
C NA,MB,MC CONTAIN COLUMN INDICES OF FIRST SECOND AND
C PRODUCT MATRICES RESPECTIVELY.
C NA IS DIMENSION OF A,B,C.  NM IS DIMENSION OF MA,MB,MC.
C
      REAL A,B,C,S
      INTEGER MA,MB,MC,NRA,NCA,NRB,NCB,LC,KA,KAF,KBF,KB,KA1,
     * KB1,JB,JAM,JBM,J1,J2,LCM,NEC,K,IT,I,J,NA,NM,LP
      DIMENSION A(NA),B(NA),C(NA),MA(NM),MB(NM),MC(NM)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CLEAR MC.
      DO 1 I = 1,NM
    1  MC(I) = 0
C CHECK THAT C DOES NOT OVER-WRITE A OR B.
      IF (MA(1).EQ.0.OR.MB(1).EQ.0) GO TO 16
C UNPACK CONTROL INFORMATION.  NRA IS NUMBER OF ROWS IN A.
      NRA     = IND(MA,1,NM)
C NCA IS NUMBER OF COLS IN A.
      NCA     = IND(MA,2,NM)
C NRB,NCB ARE NUMBER OF ROWS AND COLUMNS IN B.
      NRB     = IND(MB,1,NM)
      NCB     = IND(MB,2,NM)
C TEST FOR COMPATIBILITY.
      IF (NCA.EQ.NRB) GO TO 3
      WRITE (LP,2)
    2 FORMAT(31H A AND B INCOMPATIBLE IN MUSPMX)
      CALL EXIT
C LC IS NUMBER OF ELEMENTS IN C.
    3 LC      = 1
C KAF,KBF ARE NUMBERS OF CONTROL DATA IN MA,MB.
      KAF     = 4+NRA
      KBF     = 4+NRB
C KA,KB ARE NUMBERS OF ELEMENTS IN  FIRST I ROWS OF A,B.
      KA      = 0
C NEC IS NUMBER OF ELEMENTS IN C.
      NEC    = 0
      DO 15 I = 1,NRA
       KB     = 0
C NIRA IS NUMBER IN ROW I OF A.
       NIRA  = IND(MA,4+I,NM)
C NIRC IS NUMBER IN ROW I OF C.
       NIRC  = 0
       IF (NIRA.EQ.0) GO TO 15
       KA1    = KA+1
       KA     = KA+NIRA
       DO 14 J = 1,NCB
C NIRB IS NUMBER IN ROW I OF B.
        NIRB   = IND(MB,4+J,NM)
        IF (NIRB.EQ.0) GO TO 14
        KB1    = KB+1
        KB     = KB+NIRB
C S WILL CONTAIN I,J ELEMENT OF C.
        S      = 0.
C JB COUNTS ELEMENTS IN B.
        JB     = KB1
        DO 8 JA=KA1,KA
         JAM    = JA+KAF
    5    JBM    = JB+KBF
         J1     = IND(MA,JAM,NM)
         J2     = IND(MB,JBM,NM)
         IF (J1-J2) 8,6,7
    6    S      = S+A(JA)*B(JB)
         IF (JB.EQ.KB) GO TO 9
         JB     = JB+1
         GO TO 8
    7    IF (JB.EQ.KB) GO TO 9
         JB     = JB+1
         GO TO 5
    8   CONTINUE
C IF ELEMENT ZERO DO NOT STORE.
    9   IF (S.EQ.0.0) GO TO 14
        IF (LC.LE.NA) GO TO 11
        WRITE  (LP,10)
   10   FORMAT(17H0SIZE OF PRODUCT ,
     *   25HMATRIX EXCEEDED IN MUSPMX)
        CALL EXIT
C STORE ELEMENT AND INDEX IN C,MC.
   11   C(LC)  = S
        LCM    = LC+KAF
        IF (LCM.LE.2*NM) GO TO 13
        WRITE (LP,12)
   12   FORMAT(40H SIZE OF INDEX MATRIX EXCEEDED IN MUSPMX)
        CALL EXIT
   13   CALL IPK(J,MC,LCM,NM)
        LC     = LC+1
        NIRC   = NIRC+1
   14  CONTINUE
       NEC   = NEC+NIRC
       CALL IPK(NIRC,MC,4+I,NM)
   15 CONTINUE
C STORE CONTROL DATA IN MC.
      CALL IPK(NRA,MC,1,NM)
      CALL IPK(NCB,MC,2,NM)
      MC(2)  = NEC
      RETURN
C ERROR MESSAGE.
   16 WRITE (LP,17)
   17 FORMAT(36H IN MUSPMX PRODUCT OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE MVSPMX(A,M,AN,MN,NA,NM)
C     **********************************
C
C MOVE SPARSE MATRIX IN A TO AN.
C M,MN CONTAIN COLUMN INDICES FOR A,AN.
C NA IS DIMENSION OF A,AN.  NM IS DIMENSION OF M,MN.
C
      REAL A,AN
      INTEGER M,MN,NEA,I,NRA,N,NA,NM
      DIMENSION A(NA),M(NM),AN(NA),MN(NM)
C NEA IS NUMBER OF ELEMENTS IN A.
      NEA    = M(2)
C MOVE A.
      DO 1 I = 1,NEA
    1  AN(I) = A(I)
C NRA IS NUMBER OF ROWS IN A.
      NRA    = IND(M,1,NM)
C N IS NUMBER OF WORDS IN M.
      N      = (5+NRA+NEA)/2
C MOVE M.
      DO 2 I = 1,N
    2  MN(I) = M(I)
      RETURN
      END
C     **********************************
      SUBROUTINE PERCOL(A,M,AP,MP,IP,AT,MT,NA,NM,NP)
C     **********************************************
C PERMUTE COLUMNS OF A SPARSE MATRIX STORED IN A.
C M CONTAINS COLUMN INDICES OF A.
C AP,MP WILL CONTAIN ELEMENTS AND COLUMN INDICES OF RESULT.
C IP CONTAINS PERMUTATION--THAT IS OLD COLUMN I BECOMES NEW
C COLUMN IP(I).
C AT,MT WILL CONTAIN ROW OF A,M.
C NA IS DIMENSION OF A,AP. NM IS DIMENSION OF M,MP.
C NP IS DIMENSION OF AT,MT,IP.
C
      REAL A,AP,AT,A1
      INTEGER M,MP,IP,NR,NC,I,I1,NIR,K,
     *L,N,J,J1,LJ,N1,IFL,M1,NA,NM,NP,LP
      DIMENSION A(NA),M(NM),AP(NA),
     * MP(NM),IP(NP),AT(NP),MT(NP)
C
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CLEAR MP.
      DO 1 I = 1,NM
    1  MP(I) = 0
C CHECK THAT OUTPUT DOES NOT OVER-WRITE INPUT.
      IF (M(1).EQ.0) GO TO 15
C UNPACK CONTROL INFORMATION.
      NR      = IND(M,1,NM)
      NC      = IND(M,2,NM)
C CHECK THAT IP(K) IS WITHIN RANGE.
      DO 2 K = 1,NC
       IF (IP(K).LE.0.OR.IP(K).GT.NC) GO TO 13
    2  CONTINUE
C TRANSFER CONTROL DATA TO MP.
      I2      = 4+NR
      DO 3 I  = 1,I2
       K      = IND(M,I,NM)
    3  CALL IPK(K,MP,I,NM)
C L COUNTS ELEMENTS ALREADY PERMUTED.
      L       = 0
C I IS ROW COUNTER.
      DO 12 I = 1,NR
       N      = IND(M,4+I,NM)
       IF (N.EQ.0) GO TO 12
C STORE ROW I IN AT WITH COLUMN INDICES IN MT.
       DO 4 J  = 1,N
        J1     = 4+NR+L+J
        K      = IND(M,J1,NM)
        MT(J)  = IP(K)
        LJ     = L+J
    4   AT(J)  = A(LJ)
       IF (N.EQ.1) GO TO 10
       N1      = N-1
C IFL WILL REMAIN 0 WHEN SORTING OF ROW I COMPLETE.
    5  IFL     = 0
C SORT ELEMENTS OF ROW I IN ORDER OF INCREASING COLUMN INDEX
       DO 9 J  = 1,N1
        IF (MT(J)-MT(J+1)) 9,6,8
C ERROR MESSAGE.
    6   WRITE (LP,7)
    7   FORMAT(26H IN PERCOL 2 INDICES EQUAL)
        CALL EXIT
    8   M1     = MT(J)
        MT(J)  = MT(J+1)
        MT(J+1)= M1
        A1     = AT(J)
        AT(J)  = AT(J+1)
        AT(J+1)= A1
        IFL    = 1
    9  CONTINUE
       IF (IFL.EQ.1) GO TO 5
C TRANSFER ROW I.
   10  DO 11 J = 1,N
        LJ     = L+J
        AP(LJ) = AT(J)
        J1     = LJ+4+NR
        K      = MT(J)
   11   CALL IPK(K,MP,J1,NM)
   12  L      = L+N
      RETURN
C ERROR MESSAGES.
   13 WRITE (LP,14)
   14 FORMAT(43H IN PERCOL PERM CONTAINS INDEX OUT OF RANGE)
      CALL EXIT
   15 WRITE (LP,16)
   16 FORMAT(35H IN PERCOL OUTPUT OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE PERROW(A,M,AP,MP,IP,NA,NM,NP)
C     ****************************************
C PERMUTE ROWS OF A SPARSE MATRIX STORED IN A.
C M CONTAINS COLUMN INDICES OF INPUT MATRIX A.
C AP CONTAINS ELEMENTS OF OUTPUT MATRIX .
C MP CONTAINS COLUMN INDICES OF OUTPUT MATRIX.
C IP CONTAINS PERMUTATION--I.E. OLD ROW IP(I) BECOMES NEW
C ROW I.
C NA IS DIMENSION OF A,AP. NM IS DIMENSION OF M,MP.
C NP IS DIMENSION OF IP.
C
      REAL A,AP
      INTEGER M,MP,IP,NR,NC,I,I1,NIR,I2,
     *K,LA,LM,N1,J,J1,I3,N2,NM,NA,NP,LP
      DIMENSION A(NA),M(NM),AP(NA),MP(NM),IP(NP)
C
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CLEAR MP.
      DO 1 I = 1,NM
    1  MP(I) = 0
C CHECK THAT OUTPUT DOES NOT OVER-WRITE INPUT.
      IF (M(1).EQ.0) GO TO 10
C UNPACK CONTROL INFORMATION.
      NR     =  IND(M,1,NM)
      NC     =  IND(M,2,NM)
C RECORD NUMBERS OF ROWS ,COLUMNS AND ELEMENTS IN MP.
      DO 2 I = 1,2
    2  MP(I) = M(I)
C LA,LM COUNT ELEMENTS IN AP,MP.
      LA     = 1
      LM     = 5+NR
C PERMUTE ROWS.
      DO  9 I= 1,NR
       N1    = 0
C J IS OLD NUMBER OF NEW ROW I.
       J     = IP(I)
       K     = IND(M,4+J,NM)
       CALL IPK(K,MP,4+I,NM)
       IF (J.GT.NR.OR.J.LE.0) GO TO 3
       GO TO 5
C J OUT OF RANGE--GIVE ERROR MESSAGE.
    3  WRITE (LP,4) I
    4  FORMAT(38H IN PERROW PERM CONTAINS INDEX OUT OF ,
     *  17HRANGE IN POSITION,I3)
       CALL EXIT
C PICK OUT START AND END OF ROW J.
    5  IF (J.EQ.1) GO TO 7
       J1     = J-1
       DO 6 I3= 1,J1
    6   N1    =N1+IND(M,4+I3,NM)
C NIRJ IS NUMBER IN ROW J OF A.
    7  NIRJ   = IND(M,4+J,NM)
       IF (NIRJ.EQ.0) GO TO 9
       N2     = N1+NIRJ
       N1     = N1+1
C TRANSFER ROW J OF A,M TO ROW I OF AP,MP.
       DO 8 I3= N1,N2
        AP(LA)=  A(I3)
        K     = IND(M,4+NR+I3,NM)
        CALL IPK(K,MP,LM,NM)
        LA    = LA+1
    8   LM    = LM+1
C END OF LOOP ON I(ROW NUMBER).
    9 CONTINUE
      RETURN
C ERROR MESSAGE.
   10 WRITE (LP,11)
   11 FORMAT(35H IN PERROW OUTPUT OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE RDSPMX(A,M,NA,NM)
C     ****************************
C
C READS A SPARSE MATRIX FROM CARDS INTO ARRAY A,STORING
C COLUMN INDICES AND CONTROL DATA IN M.
C NA IS DIMENSION OF A, NM IS DIMENSION OF M.
C
      REAL A
      INTEGER NR,NC,NE,JE,JR,JF,K,MIN,M,
     * I,NIR,J,J1,J2,L1,L2,NA,NM,LIMF,IN,LP
      DIMENSION A(NA),M(NM),AIN(4),MIN(4)
C IN IS UNIT NUMBER OF CARD READER.
      IN     = 5
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C NR,NC,NE ARE NUMBERS OF ROWS,COLS,AND ELEMENTS IN A.
      READ (IN,1) NR,NC,NE
    1 FORMAT(3I5)
C JE,JR COUNT NUMBER OF ELEMENTS,ROWS.
      JE     = 1
      JR     = 1
      DO 2 I = 1,NM
    2  M(I)  = 0
      IERR   = 0
      LIMF   = 0
C JF IS NUMBER OF CONTROL DATA.
      JF     = 4+NR
C K COUNTS ELEMENTS WITHIN ROW.
      K      = 0
C AIN,MIN ARE ELEMENTS AND INDICES AS READ FROM CARD.
    3 READ (IN,4) (AIN(I),MIN(I),I=1,4)
    4 FORMAT(4(E15.8,I5))
      DO 10 I = 1,4
C CHECK FOR ROW-SENTINEL.
       IF (MIN(I).GE.90000) GO TO 9
C CHECK VALIDITY OF COLUMN-INDEX.
       IF (MIN(I).LE.NC) GO TO 6
       WRITE (LP,5) MIN(I),JE
    5  FORMAT(15H COLUMN INDEX (,I5,20H)GREATER THAN NUMBER,
     *  24H OF COL IN ELEMENT NUM. ,I5)
       IERR  = 1
C STORE ELEMENT.
    6  IF (JE.LE.NA.AND.(JE+JF).LE.2*NM) GO TO 8
       IF (LIMF.EQ.1) GO TO 10
       LIMF  = 1
       WRITE (LP,7)
    7  FORMAT(21H0ARRAY FULL IN RDSPMX)
       GO TO 10
    8  A(JE) = AIN(I)
       CALL IPK(MIN(I),M,JE+JF,NM)
       JE    = JE+1
       K     = K+1
       GO TO 10
C CHECK FOR END-OF-MATRIX SENTINEL.
    9  IF (MIN(I).EQ.99999) GO TO 11
C RECORD NUMBER OF ELEMENTS IN ROW JR OF A.
       CALL IPK(K,M,4+JR,NM)
       K       = 0
       JR      = JR+1
   10 CONTINUE
C READ NEW CARD.
      GO TO 3
C AT END OF MATRIX CHECK NUMBER OF ROWS IS AS STATED.
   11 JR     = JR-1
      IF (JR.EQ.NR) GO TO 13
      WRITE (LP,12) JR,NR
   12 FORMAT(17H NUMBER OF ROWS (,I5,17H) DOES NOT EQUAL ,
     * 15HSTATED NUMBER (,I5,1H))
      IERR   = 1
C CHECK NUMBER OF ELEMENTS.
   13 JE     = JE-1
      IF (JE.EQ.NE) GO TO 15
      WRITE (LP,14) JE,NE
   14 FORMAT(21H NUMBER OF ELEMENTS (,I5,11H) DOES NOT ,
     * 21HEQUAL STATED NUMBER (,I5,1H))
      IERR   = 1
C CHECK ASCENDING ORDER OF INDICES.
   15 J      = JF
      DO 19 I= 1,JR
C NIR IS NUMBER IN ROW I OF A.
       NIR   = IND(M,4+I,NM)
       J2    = NIR+J
       J1    = J+2
       IF (NIR.LE.1) GO TO 18
       DO 17 K = J1,J2
        L1     = IND(M,K-1,NM)
        L2     = IND(M,K,NM)
        IF (L1.LT.L2) GO TO 17
        K1     = K-J1+2
        WRITE (LP,16) K1,I
   16   FORMAT( 9H ELEMENT ,I5,8H IN ROW ,I5,11H HAS WRONG ,
     *   12HCOLUMN INDEX)
        IERR   = 1
   17  CONTINUE
   18  J     = J2
   19 CONTINUE
C STORE CONTROL DATA.
      CALL IPK(NR,M,1,NM)
      CALL IPK(NC,M,2,NM)
      M(2)   = NE
      IF (IERR.GE.1) CALL EXIT
C
      RETURN
      END
C     ********************************
      SUBROUTINE RVSPMX(A,M,IR,V,N,NA,NM)
C     ***********************************
C EXTRACTS ROW IR OF SPARSE MATRIX IN A,STORING RESULT IN
C VECTOR V IN EXTENDED FORM, I.E. INCLUDING ZERO ELEMENTS.
C M CONTAINS COLUMN INDICES OF A.
C N,NA,NM ARE DIMENSIONS OF V,A,M.
C
      REAL A,V
      INTEGER M,IR,N,NRA,I,NIRS,IR1,JM,K,J,NIRA,NA,NM,LP
      DIMENSION A(NA),M(NM),V(N)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C NRA IS NUMBER OF ROWS IN A.
      NRA    = IND(M,1,NM)
      IF (IR.GT.NRA) GO TO 2
C N IS NUMBER OF COLS IN A (EQUALS NUMBER OF ELEMENTS IN V).
      N      = IND(M,2,NM)
C CLEAR V.
      DO 1 I = 1,N
    1  V(I)  = 0.0
C NIRS WILL BE NUMBER OF ELEMENTS IN ROWS PRIOR TO IR.
      NIRS   = 0
      IF (IR-1) 2,6,4
    2 WRITE (LP,3)
    3 FORMAT(34H IN RVSPMX ROW NUMBER OUT OF RANGE)
      CALL EXIT
    4 IR1    = IR-1
      DO 5 I = 1,IR1
    5  NIRS  = NIRS+IND(M,4+I,NM)
    6 JM     = 4+NRA+NIRS
C NIRA IS NUMBER IN ROW IR.
      NIRA   = IND(M,4+IR,NM)
      IF (NIRA.LE.0) GO TO 8
C TRANSFER ELEMENTS OF ROW IR.
      DO 7 I = 1,NIRA
       K     = IND(M,JM+I,NM)
       J     = NIRS+I
    7  V(K)  = A(J)
    8 RETURN
      END
C     **********************************
      SUBROUTINE SMSPMX(A,M,AN,MN,S,NA,NM)
C     ************************************
C MULTIPLY SPARSE MATRIX IN A BY SCALAR S, STORING RESULT IN
C AN.
C M,MN CONTAIN COLUMN INDICES FOR A,AN.
C NA IS DIMENSION OF A,AN. NM IS DIMENSION OF M,MN.
C
      REAL A,AN,S
      INTEGER M,MN,NEA,I,NRA,N,NA,NM,LP
      DIMENSION A(NA),M(NM),AN(NA),MN(NM)
C
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CHECK THAT OUTPUT DOES NOT OVER-WRITE INPUT.
      MN(1)  = 0
      IF (M(1).EQ.0) GO TO 5
      IF (S.EQ.0.0) GO TO 3
C NEA IS NUMBER OF ELEMENTS IN A.
      NEA    = M(2)
C MULTIPLY A BY S.
      DO 1 I = 1,NEA
    1  AN(I) = A(I)*S
C NRA IS NUMBER OF ROWS IN A.
      NRA    = IND(M,1,NM)
C N IS NUMBER OF WORDS IN M.
      N      = (5+NRA+NEA)/2
C MOVE M.
      DO 2 I = 1,N
    2  MN(I) = M(I)
      RETURN
    3 MN(1)  = M(1)
      NRA    = IND(M,1,NM)
      K      =(5+NRA)/2
      DO 4 I = 2,K
    4  MN(I) = 0
      RETURN
C ERROR MESSAGE.
    5 WRITE (LP,6)
    6 FORMAT(35H IN SMSPMX OUTPUT OVER-WRITES INPUT,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE TRSPMX(A,M,AT,MT,NA,NM,IP,NP)
C     ****************************************
C
C TRANSPOSE A SPARSE MATRIX IN A, STORING THE RESULT IN AT.
C M,MT CONTAIN COLUMN INDICES OF A,AT.
C NA IS DIMENSION OF A,AT. NM IS DIMENSION OF M,MT.
C IP(I) WILL BE NUMBER OF ELEMENTS IN COLUMN  I OF A, ALSO
C IP(I) WILL BE POINTER TO FIRST ELEMENT IN ROW I OF AT.
C NP IS DIMENSION OF IP.
C
      INTEGER I,IC,IFC,IFR,IT,J1,IP
      REAL A,AT
      DIMENSION A(NA),M(NM),AT(NA),MT(NM),IP(NP)
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C CLEAR MT.
      DO 1 I = 1,NM
    1  MT(I) = 0
C CHECK THAT AT DOES NOT OVER-WRITE A.
      IF (M(1).EQ.0) GO TO 8
C UNPACK CONTROL INFORMATION. NR,NC,NE ARE NUMBERS OF ROWS,
C  COLUMNS AND ELEMENTS IN A.
      NR     = IND(M,1,NM)
      NC     = IND(M,2,NM)
      NE     = M(2)
C CHECK FOR POSSIBLE OVERFLOW OF MT.
      L      = 4+NC+NE
      IF (L.GT.2*NM) GO TO 10
C PACK NUMBER OF ROWS(NC), COLUMNS(NR) AND ELEMENTS OF AT.
      CALL IPK(NC,MT,1,NM)
      CALL IPK(NR,MT,2,NM)
      MT(2)  = M(2)
C IFR,IFC ARE NUMBER OF CONTROL DATA IN A,AT.
      IFR    = 4+NR
      IFC    = 4+NC
C CLEAR IP.
      DO 2 I = 1,NC
    2  IP(I) = 0
C COUNT NUMBER OF ELEMENTS IN EACH COLUMN OF A(ROW OF AT).
      DO 3 I = 1,NE
       K     = IND(M,IFR+I,NM)
    3  IP(K) = IP(K)+1
C PACK NUMBERS OF ELEMENTS IN ROWS OF AT.
      DO 4 J = 1,NC
       K     = IP(J)
    4  CALL IPK(K,MT,4+J,NM)
C SET UP POINTER TO FIRST ELEMENT IN ROW I OF AT.
       NICA1 = IP(1)
      IP(1)  = 1
      DO 5 I = 2,NC
       NICA  = IP(I)
       IP(I) = IP(I-1)+NICA1
    5  NICA1 = NICA
C PROCESS ROWS OF A. J1 IS POSITION OF FIRST ELEMENT OF
C CURRENT ROW OF A.
      J1     = 1
C I COUNTS ROWS OF A.
      DO 7 I = 1,NR
C NIRA IS NUMBER OF ELEMENTS IN CURRENT ROW OF A.
       NIRA  = IND(M,4+I,NM)
       IF (NIRA.EQ.0) GO TO 7
       J2    = J1+NIRA-1
       DO 6 J = J1,J2
C K IS COLUMN NUMBER OF J*TH ELEMENT IN A, I.E.ROW NUMBER
C  IN AT.
        K      = IND(M,IFR+J,NM)
C IT IS POSITION OF CURRENT ELEMENT IN AT.
        IT     = IP(K)
        AT(IT) = A(J)
        CALL IPK(I,MT,IT+IFC,NM)
    6   IP(K)  = IP(K)+1
    7  J1    = J2+1
      RETURN
C ERROR MESSAGES.
    8 WRITE (LP,9)
    9 FORMAT(27H IN TRSPMX AT OVER-WRITES A,
     * 34H OR INPUT HAS NO ROWS AND COLUMNS.)
      CALL EXIT
   10 WRITE (LP,11)
   11 FORMAT(27H IN TRSPMX MT WILL OVERFLOW)
      CALL EXIT
      END
C **********************************************************
      SUBROUTINE WRSPMX(A,M,TIT,NA,NM)
C     ********************************
C WRITE SPARSE MATRIX A.
C M CONTAINS COLUMN INDICES OF A.
C TIT CONTAINS DESCRIPTION OF A.
C NA,NM ARE DIMENSIONS OF A,M.
C
      REAL A,AOUT,TIT
      INTEGER M,P,I,L,NRA,JF,NIRA,J,J2,K2,K,KJ,MOUT,NA,NM,LP
      DIMENSION A(NA),M(NM),TIT(10),AOUT(5),MOUT(5)
C
C LP IS UNIT NUMBER OF LINE PRINTER.
      LP     = 6
C P IS PAGE COUNTER.
      P      = 1
C HEADING AND DESCRIPTION.
      WRITE (LP,1) (TIT(I),I=1,10),P
    1 FORMAT(23H1PRINTOUT OF SPARSE MAT,
     * 4HRIX.,13X,10A4,10X,5HPAGE ,I5//)
C L IS LINE COUNTER.
      L      = 2
C NRA IS NUMBER OF ROWS IN A.
      NRA    = IND(M,1,NM)
      JF     = 4+NRA
C NIRS IS NUMBER IN ROWS ALREADY WRITTEN.
      NIRS   = 0
C I IS ROW COUNTER.
      DO 8 I = 1,NRA
C NIRA IS NUMBER IN ROW OF A.
       NIRA  = IND(M,4+I,NM)
       IF (NIRA.EQ.0) GO TO 8
C J COUNTS ELEMENTS WRITTEN.
       J     = 0
       IF (L.LT.51) GO TO 3
C AT END OF PAGE WRITE NEW HEADING ON NEXT PAGE,UPDATING
C PAGE NUMBER.
    2  P     = P+1
       WRITE (LP,1) (TIT(K),K=1,10),P
       L     = 2
    3  WRITE (LP,4) I
    4  FORMAT(12H0ROW NUMBER ,I5//
     * 1X,5(4X,7HELEMENT,5X,3HCOL,5X))
       L     = L+4
C EXTRACT NEXT LINE OF OUTPUT.
    5  J2    = MIN0(NIRA,J+5)
       K2    = J2-J
       DO 6 K  = 1,K2
        KJ     = K+J+NIRS
        MOUT(K)= IND(M,KJ+JF,NM)
    6   AOUT(K)= A(KJ)
       WRITE (LP,7) (AOUT(K),MOUT(K),K=1,K2)
    7  FORMAT(1X,5(E15.7,I5,4X))
       L     = L+1
       J     = J+5
       IF (J.GE.NIRA) GO TO 8
       IF (L-55) 5,2,2
    8 NIRS   = NIRS+NIRA
C LAST ROW WRITTEN.
      RETURN
      END