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


C=======================================================================
C                          G E N M A T
C=======================================================================
C     The procedures here produce various kinds of arrays with
C     SUBROUTINE GENMAT controlling the generation
C=======================================================================
      SUBROUTINE GENMAT
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *           WORK1,INDIC1,NONZR1,
     *           WORK2,INDIC2,NONZR2,
     *           WORK3,INDIC3,NONZR3,
     *           SIGDIG,MATTYP,DENSI,MAXVAL,INTVAL,NONZRS,
     *           PATTER,DETERM,ISTAT,
     *           MAXNZ,MAXCO,MAXCOP,MAXRO,MAXVC,TINDIC,
     *           INDIC4,JOINT,NZROW)
      INTEGER   MAXNZ,MAXCO,MAXCOP,MAXRO,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN
      REAL*8    WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC)
      INTEGER   INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC),INDIC4(MAXVC)
      INTEGER   NONZR1,NONZR2,NONZR3,MATTYP,JOINT,NZROW
      REAL*8    MAXVAL,DENSI,DETERM
      INTEGER   INTVAL,INDSTA,ISTAT,SIGDIG,PATTER,TINDIC
      INTEGER   NONZRS,NONZRQ,J
      REAL*8    ZERO
      IF (MM .EQ. 0 .OR. NN .EQ. 0) RETURN
      ZERO=1.0D-8
      NONZRS=0

      IF (ISTAT .EQ. -1) GOTO 60

      ISTAT=0

      NONZRQ=IDINT( MM * NN * 1.0D-2 * DENSI )

      IF (NONZRQ .LT. MM*NN) NONZRQ=NONZRQ+1

      IF (PATTER .NE. 1 .AND. TINDIC .EQ. 0) GOTO 70

      CALL INITM(ABYCOL,ROWNOS,COLPNT,MM,NN,
     *           MAXNZ,MAXCOP,MAXVC)
      IF (MATTYP .NE. 0 .AND. MATTYP .NE. 1) THEN
        IF (NN .LT. MM+1) THEN
           MATTYP=1
                          ELSE
           GOTO 40
        ENDIF
      ENDIF
C     General sparse matrix (nonsingular/any), MM <= NN

      IF (MM .GT. NN) GOTO 20
      IF ( MATTYP .EQ. 0 ) THEN
        INDSTA = 1
        NONZRS=0
        GOTO 10
      ENDIF
      CALL NONSING
     *           (ABYCOL,ROWNOS,COLPNT,MM,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            SIGDIG,DENSI,MAXVAL,INTVAL,
     *            ZERO,NONZRS,DETERM,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      IF (MM .EQ. NN) GOTO 100
      INDSTA = MM + 1
   10 CONTINUE
      CALL FILLCO
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC,
     *            INDIC4,JOINT,NZROW)
      GOTO 100

C     General sparse matrix (nonsingular/any), MM > NN

   20 CONTINUE
      IF ( MATTYP .EQ. 0 ) THEN
        INDSTA = 1
        NONZRS=0
        GOTO 30
      ENDIF
      CALL NONSING
     *           (ABYCOL,ROWNOS,COLPNT,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            SIGDIG,DENSI,MAXVAL,INTVAL,
     *            ZERO,NONZRS,DETERM,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      IF (MM .EQ. NN) GOTO 100
      INDSTA = NN + 1
   30 CONTINUE
      CALL FILLRO
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC,
     *            INDIC4,JOINT,NZROW)
      GOTO 100
   40 CONTINUE
      IF (MATTYP .NE. 4) GOTO 50

C     Complete recourse matrix

      CALL INITM(ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            MAXNZ,MAXCOP,MAXVC)
      CALL NONSING
     *           (ABYCOL,ROWNOS,COLPNT,MM,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            SIGDIG,DENSI,MAXVAL,INTVAL,
     *            ZERO,NONZRS,DETERM,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      IF (MM .EQ. NN) GOTO 100
      INDSTA = MM + 1
      CALL COMPRC
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      GOTO 100
   50 CONTINUE

C     Simple recourse matrix

      IF (MATTYP .EQ. 3) GOTO 100

C     Invalid MATTYP code

      ISTAT=-1
      RETURN
   60 CONTINUE

C     Right-hand-side to given matrix, nonnegative solution should exist,
C     and objective to ensure boundedness

      CALL COLSUM
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            SIGDIG,MAXVAL,INTVAL,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      GOTO 100

C     Matrix pattern specially selected

   70 CONTINUE
      IF (PATTER .EQ. -1 .OR. PATTER .EQ. -4 ) THEN

C     Same pattern as given in ROWNOS, COLPNT

      CALL EQPATT
     *           (ABYCOL,COLPNT,NN,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      ENDIF
      IF (PATTER .EQ. 2 .OR. PATTER .EQ. -2 ) THEN

C     Subset of pattern given in ROWNOS, COLPNT

      CALL SUBPATT
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      ENDIF
      IF (PATTER .EQ. 3 .OR. PATTER .EQ. -3 ) THEN

C     Superset of pattern given in ROWNOS, COLPNT

      CALL SUPPATT
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)

      ENDIF
  100 CONTINUE
      IF ( JOINT .EQ. 1 ) THEN
         DO 80 J=1,NN
            IF ( COLPNT(J) .EQ. COLPNT(J+1) ) THEN
               INDIC4(J) = 0
                                              ELSE
               INDIC4(J) = 1
            ENDIF
   80    CONTINUE
      ENDIF
      END
C=======================================================================
      SUBROUTINE FILLCO
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC,
     *            INDIC4,JOINT,NZROW)
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN
      REAL*8    WORK1(MAXVC),WORK3(MAXVC),MAXVAL
      INTEGER   INDIC1(MAXVC),INDIC3(MAXVC),NONZR1,NONZR3,I
      INTEGER   INTVAL,NONZRS,INDSTA,ISTAT,SIGDIG
      INTEGER   INDIC4(MAXVC),JOINT,NZROW
      REAL*8    ZERO,NEWVAL
      INTEGER   NPOS1,NONZRQ,NNN,INDX1,INDX2
      REAL*8    XENTRY

      NNN = INDSTA-1
      NONZR1=1
      IF ( NZROW .GT. 0 ) THEN
         DO 5 I=1,MM
    5    INDIC4(I) = 0
      ENDIF
      DO 10 I=INDSTA,NN
         IF ( JOINT .EQ. 0 .OR.
     *        ( JOINT .EQ. 2 .AND. INDIC4(I) .EQ. 0 ) ) THEN
C     Add nonzero column
            CALL INDEX(MM,NPOS1)
            WORK1(1)=XENTRY(MAXVAL,INTVAL,SIGDIG,ZERO)
            NONZRS=NONZRS+1
            IF ( MOD(NONZRS,100) .EQ. 0 )
     *      WRITE(*,'(A,I5)')' Counter = ',NONZRS
            INDIC1(1)=NPOS1
            CALL APPCOL
     *          (ABYCOL,ROWNOS,COLPNT,NNN,
     *           WORK1,INDIC1,NONZR1,
     *           ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
                                                        ELSE
C     Add zero column
            CALL APPZC(COLPNT,NNN,ISTAT,MAXCOP)
         ENDIF
   10 CONTINUE
      IF (NONZRQ .LE. NONZRS) THEN
         NONZRQ=NONZRS
         RETURN
      ENDIF
   20 CONTINUE
      NNN=NN-INDSTA+1
      CALL INDEX(NNN,INDX2)
      INDX2=INDX2+INDSTA-1
      CALL GETCOL
     *       (ABYCOL,ROWNOS,COLPNT,NN,
     *        WORK1,INDIC1,NONZR1,
     *        INDX2,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
      CALL XINDIC(INDIC1,NONZR1,INDX1,MM,MAXVC)
      IF (INDX1 .EQ. 0) GOTO 20
      IF ( NZROW .GT. 0 ) INDIC4( INDX1 ) = 1
      CALL INSIND(  WORK1,INDIC1,NONZR1,
     *              WORK3,INDIC3,NONZR3,NEWVAL,
     *              ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,NONZRS,ISTAT,MAXVC)
      CALL INSCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK3,INDIC3,NONZR3,
     *      INDX2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      IF (NONZRS .LT. NONZRQ) GOTO 20
      IF ( NZROW .GT. 0 ) THEN
          DO 30 I=1,MM
             IF ( INDIC4(I) .EQ. 0 ) THEN
                NNN=NN-INDSTA+1
                CALL INDEX(NNN,INDX2)
                INDX2=INDX2+INDSTA-1
                CALL GETCOL
     *                     (ABYCOL,ROWNOS,COLPNT,NN,
     *                      WORK1,INDIC1,NONZR1,
     *                      INDX2,ISTAT,
     *                      MAXNZ,MAXCOP,MAXVC)
                INDX1 = I
                CALL INSIND(  WORK1,INDIC1,NONZR1,
     *                        WORK3,INDIC3,NONZR3,NEWVAL,
     *                        ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,
     *                        NONZRS,ISTAT,MAXVC)
                CALL INSCOL
     *                     ( ABYCOL,ROWNOS,COLPNT,NN,
     *                       WORK3,INDIC3,NONZR3,
     *                       INDX2,ISTAT,
     *                       MAXNZ,MAXCOP,MAXVC)
             ENDIF
   30     CONTINUE
      ENDIF
      WRITE(*,'(A,I5)')' Counter = ',NONZRS
      END
C=======================================================================
      SUBROUTINE FILLRO
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC,
     *            INDIC4,JOINT,NZROW)
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN
      REAL*8    WORK1(MAXVC),WORK3(MAXVC),MAXVAL
      INTEGER   INDIC1(MAXVC),INDIC3(MAXVC),NONZR1,NONZR3,I,J
      INTEGER   INTVAL,NONZRS,INDSTA,ISTAT,SIGDIG
      INTEGER   INDIC4(MAXVC),JOINT,NZROW
      REAL*8    ZERO,NEWVAL
      INTEGER   NONZRQ,DIMMY,INDX1,INDX2
      INTEGER   INCOL,INREST,JJ
      LOGICAL   COIN
C
      IF ( NZROW .GT. 0 ) THEN
         DO 5 I=1,MM
    5    INDIC4(I) = 0
      ENDIF
      INCOL=IDINT( (NONZRQ-NONZRS) / NN * 1.0D0 )
      INREST= NONZRQ-NONZRS-INCOL*NN
      IF ( NONZRS .EQ. 0 .AND. INCOL .EQ. 0 ) THEN
           INCOL=1
           INREST=0
      ENDIF
      IF (INCOL .LT. 0) THEN
         WRITE(*,'(A)') 'Density too low...'
         RETURN
      ENDIF
      DO 10 I=1,NN
         CALL GETCOL
     *          (ABYCOL,ROWNOS,COLPNT,NN,
     *           WORK1,INDIC1,NONZR1,
     *           I,ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
         IF (INCOL .GT. 0) THEN
           DO 20 J=1,INCOL
              CALL YINDIC(INDIC1,NONZR1,INDX1,INDSTA,MM,MAXVC)
              IF ( NZROW .GT. 0 ) INDIC4( INDX1 ) = 1
              CALL INSIND(  WORK1,INDIC1,NONZR1,
     *                      WORK3,INDIC3,NONZR3,NEWVAL,
     *                      ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,
     *                      NONZRS,ISTAT,MAXVC)
              CALL COPYVC(  WORK1,INDIC1,NONZR1,
     *                      WORK3,INDIC3,NONZR3,MAXVC)
   20      CONTINUE
         ENDIF
         IF (I .EQ. NN .AND. INREST. GT. 0) THEN
           DO 30 J=1,INREST
              CALL YINDIC(INDIC1,NONZR1,INDX1,INDSTA,MM,MAXVC)
              IF ( NZROW .GT. 0 ) INDIC4( INDX1 ) = 1
              CALL INSIND(WORK1,INDIC1,NONZR1,
     *                    WORK3,INDIC3,NONZR3,NEWVAL,
     *                    ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,NONZRS,
     *                    ISTAT,MAXVC)
              CALL COPYVC(  WORK1,INDIC1,NONZR1,
     *                      WORK3,INDIC3,NONZR3,MAXVC)
   30      CONTINUE
         ENDIF
         IF (I. LT. NN .AND. INREST .GT. 0) THEN
            IF (COIN(DIMMY)) THEN
               INREST=INREST-1
               CALL YINDIC(INDIC1,NONZR1,INDX1,INDSTA,MM,MAXVC)
               CALL INSIND(WORK1,INDIC1,NONZR1,
     *                     WORK3,INDIC3,NONZR3,NEWVAL,
     *                     ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,NONZRS,
     *                     ISTAT,MAXVC)
               CALL COPYVC(  WORK1,INDIC1,NONZR1,
     *                       WORK3,INDIC3,NONZR3,MAXVC)
            ENDIF
         ENDIF
         CALL INSCOL
     *            (ABYCOL,ROWNOS,COLPNT,NN,
     *             WORK3,INDIC3,NONZR3,
     *             I,ISTAT,
     *             MAXNZ,MAXCOP,MAXVC)
   10 CONTINUE
      IF ( NZROW .GT. 0 ) THEN
          DO 40 I=1,MM
             IF ( INDIC4(I) .EQ. 0 ) THEN
                CALL INDEX(NN,INDX2)
                CALL GETCOL
     *                     (ABYCOL,ROWNOS,COLPNT,NN,
     *                      WORK1,INDIC1,NONZR1,
     *                      INDX2,ISTAT,
     *                      MAXNZ,MAXCOP,MAXVC)
                INDX1 = I
                CALL INSIND(  WORK1,INDIC1,NONZR1,
     *                        WORK3,INDIC3,NONZR3,NEWVAL,
     *                        ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,
     *                        NONZRS,ISTAT,MAXVC)
                CALL INSCOL
     *                     ( ABYCOL,ROWNOS,COLPNT,NN,
     *                       WORK3,INDIC3,NONZR3,
     *                       INDX2,ISTAT,
     *                       MAXNZ,MAXCOP,MAXVC)
             ENDIF
   40     CONTINUE
      ENDIF
      NONZRQ=NONZRS
      WRITE(*,'(A,I5)')' Counter = ',NONZRS
      END
C=======================================================================
      SUBROUTINE NONSING
     *           (ABYCOL,ROWNOS,COLPNT,MM,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            SIGDIG,DENSI,MAXVAL,INTVAL,
     *            ZERO,NONZRS,DETERM,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
C
C     Return codes:
C             ISTAT  = 1 SPARSITY CONTRADICTS SIZE
C
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM
      REAL*8    WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC),MAXVAL,
     *          DENSI,ZERO,DETERM
      INTEGER   INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC),
     *          NONZR1,NONZR2,NONZR3
      INTEGER   INTVAL,ISTAT,SIGDIG
      INTEGER   NPOS1,NPOS2,DIMMY
      INTEGER   NONZRS,NONZTR,NONZRQ
      LOGICAL   COIN
C
C     Initialize data.
C
      ISTAT=0
C
C     Compute number of nonzeros.
C
      NONZRQ=IDINT( MM * MM * 1.0D0 * 0.01 * DENSI )
      IF (NONZRQ .LT. MM*MM) NONZRQ=NONZRQ+1
C
C     Check data.
C
      IF (MM .GT. NONZRQ) THEN
        ISTAT=1
        NONZRQ=MM
      ENDIF
C
C     Compute triangular matrix.
C
      NONZTR=MM*(MM+1)/2
      NONZTR=MIN(NONZRQ,NONZTR)
      CALL TRIANG
     *     (ABYCOL,ROWNOS,COLPNT,MM,
     *      WORK1,INDIC1,NONZR1,
     *      WORK2,INDIC2,NONZR2,
     *      WORK3,INDIC3,NONZR3,
     *      MAXVAL,INTVAL,
     *      DETERM,ZERO,NONZRS,NONZTR,SIGDIG,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
C     Shake matrix.
C
      CALL ENTROP(MM,
     *      ABYCOL,ROWNOS,COLPNT,MM,MM,NPOS1,WORK1,INDIC1,NONZR1,
     *      NPOS2,WORK2,INDIC2,NONZR2,WORK3,INDIC3,NONZR3,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C     CHFREQ=MM/5
C     CHFREQ=MAX(CHFREQ,20)
C     CALL SHAKER
C    *     (CHFREQ,
C    *      ABYCOL,ROWNOS,COLPNT,MM,MM,NPOS1,WORK1,INDIC1,NONZR1,
C    *      NPOS2,WORK2,INDIC2,NONZR2,WORK3,INDIC3,NONZR3,ISTAT,
C    *      MAXNZ,MAXCOP,MAXVC)
C
C     If sufficient amount of nonzeros achieved then exit.
C
      IF (NONZRQ .LE. NONZRS) RETURN
C
C     Propagate nonzeros by addig a multiple of a row/column to another.
C
C     Iteration starts here.
C

   30 CONTINUE
C
C     Get indices.
C
      CALL INDICS(MM,NPOS1,NPOS2)
C
C     Add multiple of one row/column to another.
C
      IF (COIN(DIMMY))  GOTO 100
          CALL GETROW
     *         (ABYCOL,ROWNOS,COLPNT,MM,MM,
     *          WORK1,INDIC1,NONZR1,
     *          NPOS1,ISTAT,
     *          MAXNZ,MAXCOP,MAXVC)
          CALL GETROW
     *         (ABYCOL,ROWNOS,COLPNT,MM,MM,
     *          WORK2,INDIC2,NONZR2,
     *          NPOS2,ISTAT,
     *          MAXNZ,MAXCOP,MAXVC)
          CALL MIXARR
     *           (WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            MM,MAXVAL,INTVAL,NONZRS,ZERO,SIGDIG,
     *            MAXNZ,MAXCOP,MAXVC)
          CALL INSROW
     *         (ABYCOL,ROWNOS,COLPNT,MM,MM,
     *          WORK2,INDIC2,NONZR2,
     *          WORK1,INDIC1,NONZR1,
     *          NPOS2,ISTAT,
     *          MAXNZ,MAXCOP,MAXVC)
          GOTO 110
  100 CONTINUE
           CALL GETCOL
     *          (ABYCOL,ROWNOS,COLPNT,MM,
     *           WORK1,INDIC1,NONZR1,
     *           NPOS1,ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
           CALL GETCOL
     *          (ABYCOL,ROWNOS,COLPNT,MM,
     *           WORK2,INDIC2,NONZR2,
     *           NPOS2,ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
          CALL MIXARR
     *           (WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            MM,MAXVAL,INTVAL,NONZRS,ZERO,SIGDIG,
     *            MAXNZ,MAXCOP,MAXVC)
          CALL INSCOL
     *         (ABYCOL,ROWNOS,COLPNT,MM,
     *          WORK2,INDIC2,NONZR2,
     *          NPOS2,ISTAT,
     *          MAXNZ,MAXCOP,MAXVC)
  110  CONTINUE
C
C     Check if done.
C
      IF (NONZRS .LT. NONZRQ) GOTO 30
C
C     Shake matrix.
C
      CALL ENTROP(MM,
     *      ABYCOL,ROWNOS,COLPNT,MM,MM,NPOS1,WORK1,INDIC1,NONZR1,
     *      NPOS2,WORK2,INDIC2,NONZR2,WORK3,INDIC3,NONZR3,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      ISTAT = 0
      WRITE(*,'(A,I5)')' Counter = ',NONZRS
      END
C=======================================================================
      SUBROUTINE TRIANG
     *           (ABYCOL,ROWNOS,COLPNT,MM,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            UPPERB,BNDTYP,
     *            DETERM,ZERO,NONZRS,NONZTR,PRECIS,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
C
C     Purpose: Generate nonsingular triangular matrix with prescribed
C              sparsity, lower and upper bounds for elements.
C     Return codes:
C              none.
C
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM
      REAL*8    WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC),UPPERB
      REAL*8    DETERM,ZERO
      INTEGER   INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC),
     *          NONZR1,NONZR2,NONZR3
      INTEGER   BNDTYP,PRECIS,ISTAT
      REAL*8    XENTRY,NEWVAL
      INTEGER   NN,CHFREQ,NPOS1,NPOS2,INDX1,INDX2
      INTEGER   NONZRS,NONZTR,I
C
C     Initialize data.
C
      ISTAT=0
      CHFREQ=MAX(MM,20)
C
C     Compute diagonal matrix.
C
      NN=1
      DETERM=1.0D0
      NONZR1=1
      NPOS1=1
      DO 10 I=1,MM
        WORK1(1)=XENTRY(UPPERB,BNDTYP,PRECIS,ZERO)
        NONZRS=NONZRS+1
         IF ( MOD(NONZRS,100) .EQ. 0 )
     *  WRITE(*,'(A,I5)')' Counter = ',NONZRS
        INDIC1(1)=I
        IF (I .GT. 1) THEN
           CALL APPCOL
     *          (ABYCOL,ROWNOS,COLPNT,NN,
     *           WORK1,INDIC1,NONZR1,
     *           ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
        ELSE
            CALL INSCOL
     *           (ABYCOL,ROWNOS,COLPNT,NN,
     *            WORK1,INDIC1,NONZR1,
     *            NPOS1,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
        ENDIF
   10 CONTINUE

C
C     If MM nonzeros required then shake matrix and exit.
C
      IF (NONZTR .EQ. MM) THEN
        CALL SHAKER
     *       (CHFREQ,
     *        ABYCOL,ROWNOS,COLPNT,MM,NN,NPOS1,WORK1,INDIC1,NONZR1,
     *        NPOS2,WORK2,INDIC2,NONZR2,WORK3,INDIC3,NONZR3,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
        RETURN
      ENDIF
C
C     Compute upper off-diagonal part.
C
   20 CONTINUE
      CALL INDEX(MM,INDX2)
      IF (INDX2 .EQ. 1) GOTO 20
      CALL GETCOL
     *       (ABYCOL,ROWNOS,COLPNT,MM,
     *        WORK1,INDIC1,NONZR1,
     *        INDX2,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
      IF (NONZR1 .EQ. INDX2) GOTO 20
      CALL XINDIC(INDIC1,NONZR1,INDX1,INDX2,MAXVC)
      IF (INDX1 .EQ. 0) GOTO 20
      CALL INSIND(  WORK1,INDIC1,NONZR1,
     *              WORK3,INDIC3,NONZR3,NEWVAL,
     *              ZERO,PRECIS,UPPERB,BNDTYP,INDX1,NONZRS,ISTAT,MAXVC)
      CALL INSCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK3,INDIC3,NONZR3,
     *      INDX2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      IF (NONZRS .LT. NONZTR) GOTO 20
C
C     Shake matrix.
C
      CALL SHAKER
     *     (CHFREQ,
     *      ABYCOL,ROWNOS,COLPNT,MM,NN,NPOS1,WORK1,INDIC1,NONZR1,
     *      NPOS2,WORK2,INDIC2,NONZR2,WORK3,INDIC3,NONZR3,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
C
      ISTAT = 0
      END
C=======================================================================
      SUBROUTINE COMPRC
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN
      REAL*8    WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC),MAXVAL
      INTEGER   INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC)
      INTEGER   NONZR1,NONZR2,NONZR3,I,J
      INTEGER   INTVAL,NONZRS,INDSTA,ISTAT,SIGDIG
      REAL*8    ZERO,NEWVAL
      INTEGER   NPOS1,NONZRQ,NNN,INDX1,INDX2,INX
      REAL*8    XENTRY
C
C     Initialize data.
C
      ISTAT=0

C     Compute -B*1

      DO 10 I=1,MM
      INDIC2(I)=-1
   10 WORK2(I)=0.0D0
      DO 20 I=1,MM
         CALL GETCOL
     *            ( ABYCOL,ROWNOS,COLPNT,NN,
     *              WORK1,INDIC1,NONZR1,
     *              I,ISTAT,
     *              MAXNZ,MAXCOP,MAXVC)
         DO 30 J=1,NONZR1
            INX=INDIC1(J)
            WORK2(INX)=WORK2(INX)-WORK1(J)
   30    CONTINUE
   20 CONTINUE

C     Put one nonzero into each remaining column.

      NNN=MM
      NONZR1=1
      DO 40 I=INDSTA,NN
         CALL INDEX(MM,NPOS1)
         IF (I .GT. INDIC2(NPOS1)) INDIC2(NPOS1)=I
         WORK1(1)=XENTRY(MAXVAL,INTVAL,SIGDIG,ZERO)
         WORK2(NPOS1)=WORK2(NPOS1)-WORK1(1)
         NONZRS=NONZRS+1
         IF ( MOD(NONZRS,100) .EQ. 0 )
     *   WRITE(*,'(A,I5)')' Counter = ',NONZRS
         INDIC1(1)=NPOS1
         CALL APPCOL
     *          (ABYCOL,ROWNOS,COLPNT,NNN,
     *           WORK1,INDIC1,NONZR1,
     *           ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
   40 CONTINUE

C     Enforce that each row of N contains at least one nonzero.

      DO 50 I=1,MM
         IF (INDIC2(I) .LT. 0 ) THEN
            NNN=NN-INDSTA+1
            CALL INDEX(NNN,INDX2)
            INDX2=INDX2+INDSTA-1
            CALL GETCOL
     *                (ABYCOL,ROWNOS,COLPNT,NN,
     *                 WORK1,INDIC1,NONZR1,
     *                 INDX2,ISTAT,
     *                 MAXNZ,MAXCOP,MAXVC)
            INDX1=I
            IF (INDX2 .GT. INDIC2(I)) INDIC2(I)=INDX2
            CALL INSIND(  WORK1,INDIC1,NONZR1,
     *                    WORK3,INDIC3,NONZR3,NEWVAL,
     *                    ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,NONZRS,
     *                    ISTAT,MAXVC)
            WORK2(I)=WORK2(I)-NEWVAL
            CALL INSCOL
     *                 (ABYCOL,ROWNOS,COLPNT,NN,
     *                  WORK3,INDIC3,NONZR3,
     *                  INDX2,ISTAT,
     *                  MAXNZ,MAXCOP,MAXVC)
         ENDIF
   50 CONTINUE

C     Fill up according to sparsity.

      IF (NONZRS .GT. NONZRQ) THEN
         NONZRQ=NONZRS
         GOTO 70
      ENDIF
   60 CONTINUE
      NNN=NN-INDSTA+1
      CALL INDEX(NNN,INDX2)
      INDX2=INDX2+INDSTA-1
      CALL GETCOL
     *       (ABYCOL,ROWNOS,COLPNT,NN,
     *        WORK1,INDIC1,NONZR1,
     *        INDX2,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
      CALL XINDIC(INDIC1,NONZR1,INDX1,MM,MAXVC)
      IF (INDX1 .EQ. 0) GOTO 60
      IF (INDX2 .GT. INDIC2(INDX1)) INDIC2(INDX1)=INDX2
      CALL INSIND(  WORK1,INDIC1,NONZR1,
     *              WORK3,INDIC3,NONZR3,NEWVAL,
     *              ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,NONZRS,ISTAT,MAXVC)
      WORK2(INDX1)=WORK2(INDX1)-NEWVAL
      CALL INSCOL
     *     (ABYCOL,ROWNOS,COLPNT,NN,
     *      WORK3,INDIC3,NONZR3,
     *      INDX2,ISTAT,
     *      MAXNZ,MAXCOP,MAXVC)
      IF (NONZRS .LT. NONZRQ) GOTO 60
   70 CONTINUE

C     Correct last element in each row to ensure completeness.

      DO 80 I=1,MM
         INX=INDIC2(I)
         CALL GETCOL
     *            (ABYCOL,ROWNOS,COLPNT,NN,
     *             WORK1,INDIC1,NONZR1,
     *             INX,ISTAT,
     *             MAXNZ,MAXCOP,MAXVC)
         DO 90 J=1,NONZR1
            IF (INDIC1(J) .EQ. I) THEN
               WORK1(J)=WORK2(I)+WORK1(J)
               CALL INSCOL
     *                  (ABYCOL,ROWNOS,COLPNT,NN,
     *                   WORK1,INDIC1,NONZR1,
     *                   INX,ISTAT,
     *                   MAXNZ,MAXCOP,MAXVC)
               GOTO 80
            ENDIF
   90    CONTINUE
   80 CONTINUE
      WRITE(*,'(A,I5)')' Counter = ',NONZRS
      END
C=======================================================================
      SUBROUTINE INSIND(  WORK1,INDIC1,NONZR1,
     *                    WORK3,INDIC3,NONZR3,NEWVAL,
     *                    ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,NONZI,
     *                    ISTAT,MAXVC)
      INTEGER   MAXVC
      REAL*8    WORK1(MAXVC),WORK3(MAXVC),MAXVAL
      INTEGER   INDIC1(MAXVC),INDIC3(MAXVC),NONZR1,NONZR3,NONZI,I
      INTEGER   INTVAL,ISTAT,SIGDIG
      REAL*8    ZERO,XENTRY,NEWVAL
      INTEGER   INDX1
      NONZR3=0
      IF ( NONZR1 .EQ. 0 ) GOTO 100
      DO 30 I=1,NONZR1
         IF ( (I .EQ. 1 .AND. INDIC1(1)   .GT. INDX1)   .OR.
     *        (I .GT. 1 .AND. INDIC1(I-1) .LT. INDX1    .AND.
     *                        INDIC1(I)   .GT. INDX1) )  THEN
           NONZR3=NONZR3+1
           WORK3(NONZR3)=XENTRY(MAXVAL,INTVAL,SIGDIG,ZERO)
           NEWVAL=WORK3(NONZR3)
           INDIC3(NONZR3)=INDX1
         ENDIF
         NONZR3=NONZR3+1
         WORK3(NONZR3)=WORK1(I)
         INDIC3(NONZR3)=INDIC1(I)
   30 CONTINUE
C     IF ( INDIC1(NONZR1) .GE. INDX1 ) GOTO 200
      IF ( INDIC1(NONZR1) .GE. INDX1 ) GOTO 300
  100 CONTINUE
      NONZR3=NONZR3+1
      WORK3(NONZR3)=XENTRY(MAXVAL,INTVAL,SIGDIG,ZERO)
      NEWVAL=WORK3(NONZR3)
      INDIC3(NONZR3)=INDX1
  200 CONTINUE
      IF (NONZR3 .NE. NONZR1+1) THEN
         WRITE(*,'(A)')' System error, sorry...'
         STOP
      ENDIF
  300 CONTINUE
      NONZI=NONZI+1
      IF ( MOD(NONZI,100) .EQ. 0 )
     *WRITE(*,'(A,I5)')' Counter = ',NONZI
      END
C=======================================================================
      SUBROUTINE COPYVC(  WORK1,INDIC1,NONZR1,
     *                    WORK3,INDIC3,NONZR3,MAXVC)
      INTEGER   MAXVC
      REAL*8    WORK1(MAXVC),WORK3(MAXVC)
      INTEGER   INDIC1(MAXVC),INDIC3(MAXVC),NONZR1,NONZR3,I
      IF (NONZR3 .EQ. 0) RETURN
      DO 10 I=1,NONZR3
         WORK1(I)=WORK3(I)
         INDIC1(I)=INDIC3(I)
   10 CONTINUE
      NONZR1=NONZR3
      END
C=======================================================================
      SUBROUTINE COLSUM
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK2,INDIC2,NONZR2,
     *            WORK3,INDIC3,NONZR3,
     *            SIGDIG,MAXVAL,INTVAL,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN
      REAL*8    WORK1(MAXVC),WORK2(MAXVC),WORK3(MAXVC),MAXVAL
      INTEGER   INDIC1(MAXVC),INDIC2(MAXVC),INDIC3(MAXVC)
      INTEGER   NONZR1,NONZR2,NONZR3,I,J
      INTEGER   INTVAL,ISTAT,SIGDIG,SIGDOG,INX
      REAL*8    CUTENT,WMAX
C
C     For a given matrix ( in ABYCOL, ROWNOS, COLPNT; dimensions MM, NN,
C     a RHS and objective is computed ensuring that a strictly
C     positive optimal solution exists.
C
C     On input: INDIC2(I) -->  negative:   >,
C                                     0:   =,
C                              positive:   <.
C     On output:
C                WORK3 --> RHS (dense vector ),
C                WORK2 --> obj (dense vector).


      SIGDOG=SIGDIG/2

C     Generate RHS

      DO 10 I=1,MM
         WORK3(I)=0.0D0
   10 CONTINUE
      NONZR3=MM
      DO 20 I=1,NN
         CALL GETCOL
     *          (ABYCOL,ROWNOS,COLPNT,NN,
     *           WORK1,INDIC1,NONZR1,
     *           I,ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
         DO 30 J=1,NONZR1
            INX=INDIC1(J)
   30    WORK3(INX)=WORK3(INX)+WORK1(J)
   20 CONTINUE
      WMAX=0.0D0
      DO 40 I=1,MM
         IF (INDIC2(I) .LT. 0) WORK3(I)=WORK3(I)-MAXVAL*1.0D-1
         IF (INDIC2(I) .GT. 0) WORK3(I)=WORK3(I)+MAXVAL*1.0D-1
         IF( DABS(WORK3(I)) .GT. WMAX) WMAX=DABS(WORK3(I))
   40 CONTINUE
      WMAX=DABS(MAXVAL)/WMAX*0.8D0
      WMAX=CUTENT(WMAX,SIGDOG)
      DO 50 I=1,MM
         WORK3(I)=WORK3(I)*WMAX
         WORK3(I)=CUTENT(WORK3(I),SIGDIG)
   50 CONTINUE
C     For the general case the things below need improvement!
      NONZR2=NN
      DO 60 I=1,NN
         CALL GETCOL
     *          (ABYCOL,ROWNOS,COLPNT,NN,
     *           WORK1,INDIC1,NONZR1,
     *           I,ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
         WORK2(I)=0.0D0
         DO 70 J=1,NONZR1
   70    WORK2(I)=WORK2(I)+WORK1(J)
   60 CONTINUE
      WMAX=0.0D0
      DO 80 I=1,NN
   80 IF( DABS(WORK2(I)) .GT. WMAX) WMAX=DABS(WORK2(I))
      WMAX=DABS(MAXVAL)/WMAX*0.8D0
      WMAX=CUTENT(WMAX,SIGDOG)
      DO 90 I=1,NN
         WORK2(I)=WORK2(I)*WMAX+MAXVAL*1.0D-1
         WORK2(I)=CUTENT(WORK2(I),SIGDIG)
   90 CONTINUE
      END
C=======================================================================
      SUBROUTINE EQPATT
     *           (ABYCOL,COLPNT,NN,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   COLPNT(MAXCOP),NN
      INTEGER   SIGDIG,INTVAL,NONZRS,ISTAT,INX,I
      REAL*8    MAXVAL,XENTRY,ZERO

      NONZRS=0
      INX=COLPNT(NN+1)
      DO 10 I=1,INX
        ABYCOL(I)=XENTRY(MAXVAL,INTVAL,SIGDIG,ZERO)
        NONZRS=NONZRS+1
        IF ( MOD(NONZRS,100) .EQ. 0 )
     *  WRITE(*,'(A,I5)')' Counter = ',NONZRS
   10 CONTINUE
      WRITE(*,'(A,I5)')' Counter = ',NONZRS
      END
C=======================================================================
      SUBROUTINE SUPPATT
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN
      REAL*8    WORK1(MAXVC),WORK3(MAXVC),MAXVAL,ZERO
      INTEGER   INDIC1(MAXVC),INDIC3(MAXVC),NONZR1,NONZR3
      INTEGER   INTVAL,NONZRS,NONZRQ,INDSTA,ISTAT,SIGDIG
      REAL*8    NEWVAL
      INTEGER   NNN,INDX1,INDX2

      NNN = COLPNT(NN+1)
      IF ( NNN .LT. NONZRQ) THEN
         NONZRS=NNN
   10    CONTINUE
         CALL INDEX(NN,INDX2)
         CALL GETCOL
     *          (ABYCOL,ROWNOS,COLPNT,NN,
     *           WORK1,INDIC1,NONZR1,
     *           INDX2,ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
         CALL XINDIC(INDIC1,NONZR1,INDX1,MM,MAXVC)
         IF (INDX1 .EQ. 0) GOTO 10
         CALL INSIND(  WORK1,INDIC1,NONZR1,
     *                 WORK3,INDIC3,NONZR3,NEWVAL,
     *                 ZERO,SIGDIG,MAXVAL,INTVAL,INDX1,NONZRS,
     *                 ISTAT,MAXVC)
         CALL INSCOL
     *       (ABYCOL,ROWNOS,COLPNT,NN,
     *        WORK3,INDIC3,NONZR3,
     *        INDX2,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
         IF (NONZRS .LT. NONZRQ) GOTO 10
      ENDIF
      CALL EQPATT
     *           (ABYCOL,COLPNT,NN,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,NONZRS,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      WRITE(*,'(A,I5)')' Counter = ',NONZRS
      END
C=======================================================================
      SUBROUTINE SUBPATT
     *           (ABYCOL,ROWNOS,COLPNT,MM,NN,
     *            WORK1,INDIC1,NONZR1,
     *            WORK3,INDIC3,NONZR3,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,
     *            NONZRS,NONZRQ,INDSTA,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
      INTEGER   MAXNZ,MAXCOP,MAXVC
      REAL*8    ABYCOL(MAXNZ)
      INTEGER   ROWNOS(MAXNZ),COLPNT(MAXCOP),MM,NN
      REAL*8    WORK1(MAXVC),WORK3(MAXVC),MAXVAL,ZERO
      INTEGER   INDIC1(MAXVC),INDIC3(MAXVC),NONZR1,NONZR3,I
      INTEGER   INTVAL,NONZRS,NONZRQ,INDSTA,ISTAT,SIGDIG
      INTEGER   NNN,INDX1,INDX2

      NNN = COLPNT(NN+1)
      IF ( NNN .GT. NONZRQ) THEN
         NONZRS=NNN
   10    CONTINUE
         CALL INDEX(NN,INDX2)
         CALL GETCOL
     *          (ABYCOL,ROWNOS,COLPNT,NN,
     *           WORK1,INDIC1,NONZR1,
     *           INDX2,ISTAT,
     *           MAXNZ,MAXCOP,MAXVC)
         IF (NONZR1 .EQ. 0) GOTO 10
         CALL INDEX(NONZR1,INDX1)
         NONZR3=0
         DO 20 I=1,NONZR1
            IF (I .EQ. INDX1) GOTO 20
            NONZR3=NONZR3+1
            WORK3(NONZR3)=WORK1(I)
            INDIC3(NONZR3)=INDIC1(I)
   20    CONTINUE
         CALL INSCOL
     *       (ABYCOL,ROWNOS,COLPNT,NN,
     *        WORK3,INDIC3,NONZR3,
     *        INDX2,ISTAT,
     *        MAXNZ,MAXCOP,MAXVC)
         NONZRS=NONZRS-1
         IF (NONZRS .EQ. 0) GOTO 30
         IF (NONZRS .GT. NONZRQ) GOTO 10
      ENDIF
      CALL EQPATT
     *           (ABYCOL,COLPNT,NN,
     *            ZERO,SIGDIG,MAXVAL,INTVAL,NONZRS,ISTAT,
     *            MAXNZ,MAXCOP,MAXVC)
   30 CONTINUE
      WRITE(*,'(A,I5)')' Counter = ',NONZRS
      END
C=======================================================================

