From K193301%CZHRZU1A.BitNet@wigate.nic.wisc.edu Wed Nov 17 17:38:55 1993
Received: from cs.wisc.edu by tilsit.cs.wisc.edu; Wed, 17 Nov 93 17:38:53 -0600
Received: from vms3.macc.wisc.edu by cs.wisc.edu; Wed, 17 Nov 93 17:38:50 -0600
Received: from CEARN by WISCMACC.BitNet; Wed, 17 Nov 93 17:39 CDT
Received: from CZHRZU1A.BITNET (NJE origin K193301@CZHRZU1A)
          by CEARN.CERN.CH (LMail V1.1d/1.7f) with RFC822 id 9075;
          Thu, 18 Nov 1993 00:36:41 +0100
Message-Id: <23111717394158@WISCMACC.BitNet>
Date: Thu, 18 Nov 93 00:36:58 GMT
From: "P. Kall" <K193301%CZHRZU1A.BitNet@wigate.nic.wisc.edu>
Subject: sparse.for
To: Spyros Kontogiorgis <kontogio>
*Destination: CS.WISC.EDU
Status: R


C=====================================================================
C                           S P A R S E
C=====================================================================
C     This file contains procedures for manipulating sparse matrices
C=====================================================================
      SUBROUTINE GETCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK1,INDIC1,NONZR1,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Retrieves column with index NPOS1 from the data
C              structure into  WORK1,INDIC1,NONZR1.
C
C     Method : straightforward
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8   ABYCOL(MAXNZ),WORK1(MAXVC)
      INTEGER  ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC1(MAXVC)
      INTEGER  NN,NONZR1,NPOS1,ISTAT,I

      NONZR1 = 0
      IF (NPOS1 .LT. 1 .OR. NPOS1 .GT. NN) THEN
        ISTAT=-1
        RETURN
      END IF
      IF (COLPNT(NPOS1) .EQ. COLPNT(NPOS1+1)) THEN
        NONZR1 = 0
        ISTAT=1
        RETURN
      END IF
      DO 10 I=COLPNT(NPOS1)+1,COLPNT(NPOS1+1)
         WORK1(I-COLPNT(NPOS1)) =ABYCOL(I)
         INDIC1(I-COLPNT(NPOS1))=ROWNOS(I)
         NONZR1=NONZR1+1
   10 CONTINUE
      ISTAT = 0
      END
C=====================================================================
      SUBROUTINE GETROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK2,INDIC2,NONZR2,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Retrieves row with index NPOS2 from the data structure
C              into  WORK2,INDIC2,NONZR2.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false row index;
C                    1  -- row does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8   ABYCOL(MAXNZ),WORK2(MAXVC)
      INTEGER  ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC2(MAXVC)
      INTEGER  MM,NN,NONZR2,NPOS2,ISTAT,I,J

      ISTAT = 0
      NONZR2 = 0
      IF (NPOS2 .LT. 1 .OR. NPOS2 .GT. MM) THEN
        ISTAT=-1
        RETURN
      END IF
      DO 10 I=1,NN
        IF ((COLPNT(I+1)-COLPNT(I)) .NE. 0) THEN
          DO 15  J=COLPNT(I)+1,COLPNT(I+1)
            IF (ROWNOS(J) .EQ. NPOS2) THEN
              NONZR2 = NONZR2+1
              WORK2(NONZR2)  = ABYCOL(J)
              INDIC2(NONZR2) = I
            ENDIF
   15     CONTINUE
        ENDIF
   10 CONTINUE
      IF (NONZR2 .EQ. 0) ISTAT=1
      RETURN
      END
C=====================================================================
      SUBROUTINE DELCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Deletes a column with index NPOS1.
C              (the number of columns is not changed)
C
C     Method:  straigthforward
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8   ABYCOL(MAXNZ)
      INTEGER  ROWNOS(MAXNZ),COLPNT(MAXCOP)
      INTEGER  NN,NPOS1,ISTAT,I
      INTEGER  ISHIFT

      IF (NPOS1 .LT. 1 .OR. NPOS1 .GT. NN) THEN
        ISTAT=-1
        RETURN
      END IF
      IF (COLPNT(NPOS1) .EQ. COLPNT(NPOS1+1)) THEN
        ISTAT=1
        RETURN
      END IF
C shift remainder of the structure to the left; the shift distance
C is the column length.
      ISHIFT=COLPNT(NPOS1+1)-COLPNT(NPOS1)
      IF (COLPNT(NN+1) .GT. COLPNT(NPOS1+1)) THEN
        DO 10 I=COLPNT(NPOS1+1)+1,COLPNT(NN+1)
           ABYCOL(I-ISHIFT)=ABYCOL(I)
           ROWNOS(I-ISHIFT)=ROWNOS(I)
   10   CONTINUE
      END IF
      IF (NN .GT. NPOS1) THEN
        DO 20 I=NPOS1+2,NN+1
          IF (COLPNT(I) .GT. 0)
     *       COLPNT(I)=COLPNT(I)-ISHIFT
   20   CONTINUE
      END IF
      COLPNT(NPOS1+1)=COLPNT(NPOS1)
      ISTAT=0
      END
C=====================================================================
      SUBROUTINE DELROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Deletes a row with index NPOS2.
C              (The number of rows is not changed)
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK3(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC3(MAXVC)
      INTEGER NONZR3,NPOS2,MM,NN,ISTAT,I,J
      INTEGER IND

      IF (NPOS2 .LT. 1 .OR. NPOS2 .GT. MM) THEN
        ISTAT=-1
        RETURN
      END IF
      CALL GETROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      IF (ISTAT .EQ. 1) RETURN
      DO 10 I=1,NONZR3
         IND=0
         DO 20 J=COLPNT(INDIC3(I))+1,COLPNT(INDIC3(I)+1)
   20    IF (ROWNOS(J) .EQ. NPOS2) IND=J
         IF (IND .EQ. 0) THEN
           WRITE(*,'(A)')' Program error... (1)'
           RETURN
         ENDIF
         DO 30 J=IND+1,COLPNT(NN+1)
            ABYCOL(J-1)=ABYCOL(J)
            ROWNOS(J-1)=ROWNOS(J)
   30    CONTINUE
         DO 40 J=INDIC3(I),NN
   40    COLPNT(J+1)=COLPNT(J+1)-1
   10 CONTINUE
      ISTAT=0
      END
C=====================================================================
      SUBROUTINE INSCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK1,INDIC1,NONZR1,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Inserts for column index NPOS1 a column in the matrix,
C              the column to be inserted is given by WORK1,INDIC1,NONZR
C              If a column with index NPOS1 already exists in the
C              matrix, then it will be replaced by the new column.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8   ABYCOL(MAXNZ),WORK1(MAXVC)
      INTEGER  ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC1(MAXVC)
      INTEGER  NONZR1,NN,NPOS1,ISTAT
      INTEGER  ISHIFT,INX,I

      IF (NPOS1 .LT. 1 .OR. NPOS1 .GT. NN) THEN
        ISTAT=-1
        RETURN
      END IF
      ISHIFT=NONZR1
      IF (COLPNT(NPOS1) .NE. COLPNT(NPOS1+1))
     *    CALL DELCOL
     *         (ABYCOL,ROWNOS,COLPNT,NN,
     *          NPOS1,ISTAT,
     *          MAXNZ,MAXCOP,MAXVC)
      IF (COLPNT(NPOS1+1) .LT. COLPNT(NN+1)) THEN
         DO 10 I=COLPNT(NPOS1+1)+1,COLPNT(NN+1)
            INX=COLPNT(NN+1)-I+COLPNT(NPOS1+1)+1
            ABYCOL(INX+ISHIFT)=ABYCOL(INX)
            ROWNOS(INX+ISHIFT)=ROWNOS(INX)
   10    CONTINUE
      ENDIF
      DO 30 I=1,ISHIFT
         ABYCOL(COLPNT(NPOS1)+I)=WORK1(I)
         ROWNOS(COLPNT(NPOS1)+I)=INDIC1(I)
   30 CONTINUE
      DO 20 I=NPOS1+1,NN+1
   20 COLPNT(I)=COLPNT(I)+ISHIFT
      ISTAT=0
      END
C=====================================================================
      SUBROUTINE INSROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK2,INDIC2,NONZR2,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Inserts for row index NPOS2 a row in the matrix,
C              the row to be inserted is given by WORK2,INDIC2,NONZR2.
C              If a row with index NPOS2 already exists in the matrix,
C              then it will be replaced by the new row.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false row index;
C                    1  -- row does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK2(MAXVC),WORK3(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC2(MAXVC),INDIC3(MAXVC)
      INTEGER NONZR2,NONZR3,MM,NN,NPOS2,ISTAT
      INTEGER IND,I,J

      ISTAT = 0
      IF (NPOS2 .LT. 1 .OR. NPOS2 .GT. MM) THEN
        ISTAT = -1
        RETURN
      END IF
               CALL DELROW
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK3,INDIC3,NONZR3,
     *            NPOS2,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      DO 10 I=1,NONZR2
        IF (COLPNT(INDIC2(I)) .EQ. COLPNT(INDIC2(I)+1)) THEN
          IND = COLPNT(INDIC2(I))+1
        ELSE
          DO 20 J=COLPNT(INDIC2(I))+1,COLPNT(INDIC2(I)+1)
            IF (ROWNOS(J) .GT. NPOS2) THEN
              IND = J
              GOTO 21
            ENDIF
   20     CONTINUE
          IND = COLPNT(INDIC2(I)+1)+1
   21     CONTINUE
          IF (IND .EQ. 0) THEN
            WRITE(*,'(A)')' Program error... (1)'
            RETURN
          ENDIF
        ENDIF
        DO 30 J=COLPNT(NN+1),IND,-1
           ABYCOL(J+1) = ABYCOL(J)
           ROWNOS(J+1) = ROWNOS(J)
   30   CONTINUE
        ROWNOS(IND) = NPOS2
        ABYCOL(IND) = WORK2(I)
        DO 40 J=INDIC2(I),NN
   40   COLPNT(J+1)=COLPNT(J+1)+1
  10  CONTINUE
      END
C=====================================================================
      SUBROUTINE APPCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK1,INDIC1,NONZR1,
     *      ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Appends matrix with a new column,
C              the column to be appended is given by WORK1,INDIC1,NONZR
C
C     Return codes:
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK1(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC1(MAXVC)
      INTEGER NN,NONZR1,ISTAT,I

      DO 10 I=COLPNT(NN+1)+1,COLPNT(NN+1)+NONZR1
         ABYCOL(I)=WORK1(I-COLPNT(NN+1))
         ROWNOS(I)=INDIC1(I-COLPNT(NN+1))
   10 CONTINUE
      NN=NN+1
      COLPNT(NN+1)=COLPNT(NN)+NONZR1
      ISTAT=0
      END
C=====================================================================
      SUBROUTINE APPZC(COLPNT,NN,ISTAT,MAXCOP)
C
C     Purpose: Appends matrix with a zero column.
C
C     Return codes:
C
      INTEGER MAXCOP
      INTEGER COLPNT(MAXCOP)
      INTEGER NN,ISTAT
C
      NN=NN+1
      COLPNT(NN+1)=COLPNT(NN)
      ISTAT=0
      END
C=====================================================================
      SUBROUTINE EXCROW
     *      (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *       WORK1,INDIC1,NONZR1,
     *       WORK2,INDIC2,NONZR2,
     *       WORK3,INDIC3,NONZR3,
     *       NPOS1,NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Exchange for row index NPOS1 and NPOS2 the two rows
C              in the matrix.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false row index;
C                    1  -- row does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),
     *        INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC)
      INTEGER NONZR1,NONZR2,NONZR3,NPOS1,NPOS2
      INTEGER MM,NN,ISTAT

      ISTAT = 0
      IF (NPOS1 .LT. 1 .OR. NPOS1 .GT. MM) THEN
        ISTAT = -1
        RETURN
      END IF
      IF (NPOS2 .LT. 1 .OR. NPOS2 .GT. MM) THEN
        ISTAT = -1
        RETURN
      END IF
      CALL GETROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK1,INDIC1,NONZR1,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      CALL GETROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK2,INDIC2,NONZR2,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      CALL INSROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK1,INDIC1,NONZR1,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      CALL INSROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK2,INDIC2,NONZR2,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      END
C=====================================================================
      SUBROUTINE EXCCOL
     *      (ABYCOL,ROWNOS,COLPNT,NN,
     *       WORK1,INDIC1,NONZR1,
     *       WORK2,INDIC2,NONZR2,
     *       NPOS1,NPOS2,ISTAT,
     *       MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Exchange for column index NPOS1 and NPOS2 the
C              two columns in the matrix.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK1(MAXVC),WORK2(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC1(MAXVC),INDIC2(MAXVC)
      INTEGER NN,NONZR1,NONZR2,NPOS1,NPOS2,ISTAT

      ISTAT = 0
      IF (NPOS1 .LT. 1 .OR. NPOS1 .GT. NN) THEN
        ISTAT = -1
        RETURN
      END IF
      IF (NPOS2 .LT. 1 .OR. NPOS2 .GT. NN) THEN
        ISTAT = -1
        RETURN
      END IF
      CALL GETCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK1,INDIC1,NONZR1,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      CALL GETCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK2,INDIC2,NONZR2,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      CALL INSCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK1,INDIC1,NONZR1,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      CALL INSCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK2,INDIC2,NONZR2,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      END
C=====================================================================
      SUBROUTINE ADMCOL
     *       (ABYCOL,ROWNOS,COLPNT,NN,
     *        WORK1,INDIC1,NONZR1,
     *        WORK2,INDIC2,NONZR2,
     *        WORK3,INDIC3,NONZR3,
     *        NPOS1,NPOS2,MULTIP,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Adds the multiple of the column at index pos. NPOS2
C              to the column at pos. NPOS1.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),
     *        INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC)
      INTEGER NONZR1,NONZR2,NONZR3,NN,NPOS1,NPOS2,ISTAT
      REAL*8  MULTIP
      INTEGER IDX1,IDX2,IDX,I

      ISTAT = 0
      IF (NPOS1 .LT. 1 .OR. NPOS1 .GT. NN) THEN
        ISTAT = -1
        RETURN
      END IF
      IF (NPOS2 .LT. 1 .OR. NPOS2 .GT. NN) THEN
        ISTAT = -1
        RETURN
      END IF
      CALL GETCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      WORK3(NONZR3+1)=9999
      CALL GETCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK2,INDIC2,NONZR2,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      WORK2(NONZR2+1)=9999
C Do nothing for second vector beeing zero
      IF (NONZR2 .GT. 0 .AND. ABS(MULTIP) .GT. 1.E-5) THEN
C Otherwise generate multiple of second vector into itself.
        DO 90  I=1,NONZR2
          WORK2(I)  = MULTIP * WORK2(I)
 90     CONTINUE
C For first vector beeing zero, copy multiplied second into first;
C otherwise add WORK3 with second into first.
        IF (NONZR3 .EQ. 0) THEN
          NONZR1 = NONZR2
          DO 95  I=1,NONZR2
            WORK1(I)  = WORK2(I)
            INDIC1(I) = INDIC2(I)
  95      CONTINUE
        ELSE
          IDX  = 0
          IDX1 = 1
          IDX2 = 1
 92       CONTINUE
          IDX = IDX + 1
          IF (INDIC3(IDX1) .EQ. INDIC2(IDX2)) THEN
            WORK1(IDX)  = WORK3(IDX1)+WORK2(IDX2)
            INDIC1(IDX) = INDIC3(IDX1)
            IDX1 = IDX1+1
            IDX2 = IDX2+1
          ELSE IF (INDIC3(IDX1) .LT. INDIC2(IDX2)) THEN
            WORK1(IDX)  = WORK3(IDX1)
            INDIC1(IDX) = INDIC3(IDX1)
            IDX1 = IDX1+1
          ELSE IF (INDIC3(IDX1) .GT. INDIC2(IDX2)) THEN
            WORK1(IDX)  = WORK2(IDX2)
            INDIC1(IDX) = INDIC2(IDX2)
            IDX2 = IDX2+1
          ENDIF
          IF (IDX1 .LE. NONZR3 .OR. IDX2 .LE. NONZR2 ) GOTO 92
          NONZR1 = IDX
        ENDIF
        CALL INSCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK1,INDIC1,NONZR1,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      ENDIF
      END
C=====================================================================
      SUBROUTINE ADMROW
     *       (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *        WORK1,INDIC1,NONZR1,
     *        WORK2,INDIC2,NONZR2,
     *        WORK3,INDIC3,NONZR3,
     *        NPOS1,NPOS2,MULTIP,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Multiplies the row with index NPOS2 by MULTIP and
C              adds the result to the row with index NPOS1.
C
C     Method : Retrieves the row with index NPOS1 from the data
C              structure into HROW,HCOLP,NONZHR. Retrieves the row
C              with index NPOS2 from the data structure into
C              VVROW,VVCOLP,NONZVV. For MULTIP not equal zero and the
C              row with index NPOS2 not the zero vector multiply into
C              VVROW,VVCOLP,NONZVV and add into VROW,VCOLP,NONZV.
C              Then insert VROW,VCOLP,NONZV at position NPOS1.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),
     *        INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC)
      INTEGER NONZR1,NONZR2,NONZR3,MM,NN,NPOS1,NPOS2,ISTAT
      REAL*8  MULTIP
      INTEGER IDX1,IDX2,IDX,I

      ISTAT = 0
      IF (NPOS1 .LT. 1 .OR. NPOS1 .GT. MM) THEN
        ISTAT = -1
        RETURN
      END IF
      IF (NPOS2 .LT. 1 .OR. NPOS2 .GT. MM) THEN
        ISTAT = -1
        RETURN
      END IF
      CALL GETROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      CALL GETROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK2,INDIC2,NONZR2,
     *      NPOS2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      INDIC3(NONZR3+1)=9999
      INDIC2(NONZR2+1)=9999
C Do nothing for second vector being zero
      IF (NONZR2 .GT. 0 .AND. ABS(MULTIP) .GT. 1.E-5) THEN
C Otherwise generate multiple of second vector into itself.
        DO 90  I=1,NONZR2
          WORK2(I)  = MULTIP * WORK2(I)
 90     CONTINUE
C For first vector beeing zero, copy multiplied second into first;
C otherwise add temp. H... plus second into first.
        IF (NONZR3 .EQ. 0) THEN
          NONZR1 = NONZR2
          DO 95  I=1,NONZR2
            WORK1(I)  = WORK2(I)
            INDIC1(I) = INDIC2(I)
  95      CONTINUE
        ELSE
          IDX  = 0
          IDX1 = 1
          IDX2 = 1
 92       CONTINUE
          IDX = IDX + 1
          IF (INDIC3(IDX1) .EQ. INDIC2(IDX2)) THEN
            WORK1(IDX)  = WORK3(IDX1)+WORK2(IDX2)
            INDIC1(IDX) = INDIC3(IDX1)
            IDX1 = IDX1+1
            IDX2 = IDX2+1
          ELSE IF (INDIC3(IDX1) .LT. INDIC2(IDX2)) THEN
            WORK1(IDX)  = WORK3(IDX1)
            INDIC1(IDX) = INDIC3(IDX1)
            IDX1 = IDX1+1
          ELSE IF (INDIC3(IDX1) .GT. INDIC2(IDX2)) THEN
            WORK1(IDX)  = WORK2(IDX2)
            INDIC1(IDX) = INDIC2(IDX2)
            IDX2 = IDX2+1
          ENDIF
          IF (IDX1 .LE. NONZR3 .OR. IDX2 .LE. NONZR2 ) GOTO 92
          NONZR1 = IDX
        ENDIF
      CALL INSROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK1,INDIC1,NONZR1,
     *      WORK3,INDIC3,NONZR3,
     *      NPOS1,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      ENDIF
      END
C======================================================================
      SUBROUTINE INITM (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *                  MAXNZ,MAXCOP,MAXVC)
      INTEGER  MAXNZ,MAXCOP,MAXVC
      DOUBLE PRECISION  ABYCOL(MAXNZ)
      INTEGER           ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN,I
      DO 10 I=1,NN+1
   10 COLPNT(I)=0
      END
C======================================================================
      SUBROUTINE INPMATR
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,NONUL,IUNIT,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Loading a matrix from a file with unit number IUNIT.
C              The matrix is assumed to be stored in the following
C              format:
C              Record 1.       (3I5)
C                         MM    :  number of rows;
C                         NN    :  number of columns;
C                         NONUL :  number of nonzeros.
C              Records next.   (4F20.9)
C                         ABYCOL
C              Records next.   (20I4)
C                         ROWNOS
C              Records next.   (20I4)
C                         COLPNT.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- file not existent;
C                    1  -- data format error.
C
      INTEGER MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),J
      INTEGER MM,NN,NONUL,IUNIT,ISTAT
      READ (IUNIT,'(3I5)') MM,NN,NONUL
      READ (IUNIT,'(4F20.9)') (ABYCOL(J),J=1,NONUL)
      READ (IUNIT,'(20I4)') (ROWNOS(J),J=1,NONUL)
      READ (IUNIT,'(20I4)') (COLPNT(J),J=1,NN+1)
      ISTAT = 0
      END
C=====================================================================
      SUBROUTINE OUTMAX
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,NONUL,OUNIT,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Saving a matrix in a file with unit number OUNIT.
C              The data format corresponds to input of SINGXMP.
C              The matrix will be stored in the following format.
C              Record 1.       (3I5)
C                         MM    :  number of rows;
C                         NN    :  number of columns;
C                         NONUL :  number of nonzeros.
C              Records next. (for each column)
C                   Record 1.            (I5)
C                         NONULC:  number of nonzeros in column.
C                   Records next.        (I5,F20.9)
C                         Row index
C                         Entry
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8 ABYCOL(MAXNZ)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP)
      INTEGER MM,NN,NONUL,OUNIT,ISTAT
      INTEGER NONULC,INDIC,I,J
      WRITE (OUNIT,'(3I5)') MM,NN,NONUL
      DO 10 I=1,NN
         NONULC=COLPNT(I+1)-COLPNT(I)
         WRITE(OUNIT,'(I5)') NONULC
         DO 10  J=1,NONULC
            INDIC=COLPNT(I)+J
            WRITE(OUNIT,'(I5,F20.9)') ROWNOS(INDIC),ABYCOL(INDIC)
   10 CONTINUE
      ISTAT = 0
      END
C=====================================================================
      SUBROUTINE OUTMAT
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,NONUL,OUNIT,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Saving a matrix in a file with unit number OUNIT.
C              The matrix will be stored in the following format.
C              Record 1.       (3I5)
C                         MM    :  number of rows;
C                         NN    :  number of columns;
C                         NONUL :  number of nonzeros.
C              Records next.   (4F20.9)
C                         ABYCOL
C              Records next.   (20I4)
C                         ROWNOS
C              Records next.   (20I4)
C                         COLPNT.
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C
      INTEGER  MAXNZ,MAXCOP,MAXVC
      REAL*8 ABYCOL(MAXNZ)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP)
      INTEGER MM,NN,NONUL,OUNIT,ISTAT
      INTEGER J
      WRITE (OUNIT,'(3I5)') MM,NN,NONUL
      WRITE (OUNIT,'(4F20.9)') (ABYCOL(J),J=1,NONUL)
      WRITE (OUNIT,'(20I4)') (ROWNOS(J),J=1,NONUL)
      WRITE (OUNIT,'(20I4)') (COLPNT(J),J=1,NN+1)
      ISTAT = 0
      END
C=====================================================================
      SUBROUTINE SUMROW
     *       (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *        WORK2,INDIC2,NONZR2,
     *        FULLRW,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Adds the elements of a row.
C
C     Method :
C
C     Return codes:
C
C           ISTAT =  0  -- normal completion;
C                   -1  -- false column index;
C                    1  -- column does not exist in matrix.
C
      INTEGER MAXNZ,MAXCOP,MAXVC
      REAL*8  ABYCOL(MAXNZ),WORK2(MAXVC)
      INTEGER ROWNOS(MAXNZ),COLPNT(MAXCOP),INDIC2(MAXVC)
      INTEGER NONZR2,MM,NN,ISTAT,I,J
      REAL*8  FULLRW(MAXVC)

      DO 10 I=1,MM
         FULLRW(I)=0.0D0
         CALL GETROW
     *     (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *      WORK2,INDIC2,NONZR2,
     *      I,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
         IF (NONZR2 .EQ. 0) GOTO 10
         DO 20 J=1,NONZR2
            FULLRW(I)=FULLRW(I)+WORK2(J)
   20    CONTINUE
   10 CONTINUE
      END
C=======================================================================

