* File: IDSM.for     Graph manipulation subroutines
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1992-04-07 CLL Changing BWA() from LOGICAL to INTEGER.
c>> 1991-05-21 CLL @ JPL Changing (1) to (*) in declarations.
c>> 1990-03-20 CLL @ JPL
c>> 1990-02-21 CLL @ JPL
*       IDSM         Graph manipulation subroutines from
*       Argonne National Laboratory (MINPACK Project)
*       used by the separable versions of the David Gay & Linda
*       Kaufman package for nonlinear least-squares.
*       Contains IDSM and 7 subrs needed by IDSM:
*       I7RTDT, IS7ETR, ID7EGR, M7SLO, M7SEQ, I7DO, N7MSRT
*       IDSM is called by the driver subroutines that need to compute
*       finite-difference approximations to partial derivatives.
*       These are [D/S]SFU and [D/S]SFB.
*       (The names IDSM, I7RTDT, IS7ETR, & ID7EGR were formerly
*                   DSM, S7RTDT,  S7ETR, &  D7EGR.
*       These subrs do not use any real or double precision variables.
c 1992-04-07 CLL Changing BWA() from LOGICAL to INTEGER because the
c    subrs that call IDSM pass it an integer array for BWA().
*
      SUBROUTINE IDSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP,
     *               INFO,IPNTR,JPNTR,IWA,LIWA,BWA)
      INTEGER M,N,NPAIRS,MAXGRP,MINGRP,INFO,LIWA
      INTEGER INDROW(NPAIRS),INDCOL(NPAIRS),NGRP(N)
      INTEGER IPNTR(*),JPNTR(*),IWA(LIWA)
      INTEGER BWA(N)
C     **********
C
C     SUBROUTINE IDSM
C
C     THE PURPOSE OF IDSM IS TO DETERMINE AN OPTIMAL OR NEAR-
C     OPTIMAL CONSISTENT PARTITION OF THE COLUMNS OF A SPARSE
C     M BY N MATRIX A.
C
C     THE SPARSITY PATTERN OF THE MATRIX A IS SPECIFIED BY
C     THE ARRAYS INDROW AND INDCOL. ON INPUT THE INDICES
C     FOR THE NON-ZERO ELEMENTS OF A ARE
C
C           INDROW(K),INDCOL(K), K = 1,2,...,NPAIRS.
C
C     THE (INDROW,INDCOL) PAIRS MAY BE SPECIFIED IN ANY ORDER.
C     DUPLICATE INPUT PAIRS ARE PERMITTED, BUT THE SUBROUTINE
C     ELIMINATES THEM.
C
C     THE SUBROUTINE PARTITIONS THE COLUMNS OF A INTO GROUPS
C     SUCH THAT COLUMNS IN THE SAME GROUP DO NOT HAVE A
C     NON-ZERO IN THE SAME ROW POSITION. A PARTITION OF THE
C     COLUMNS OF A WITH THIS PROPERTY IS CONSISTENT WITH THE
C     DIRECT DETERMINATION OF A.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE IDSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP,
C                      INFO,IPNTR,JPNTR,IWA,LIWA,BWA)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF A.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       NPAIRS IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE
C         NUMBER OF (INDROW,INDCOL) PAIRS USED TO DESCRIBE THE
C         SPARSITY PATTERN OF A.
C
C       INDROW IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDROW
C         MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS OF A.
C         ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING
C         COLUMN INDICES ARE IN NON-DECREASING ORDER. THE COLUMN
C         INDICES CAN BE RECOVERED FROM THE ARRAY JPNTR.
C
C       INDCOL IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDCOL
C         MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS OF
C         A. ON OUTPUT INDCOL IS PERMUTED SO THAT THE CORRESPONDING
C         ROW INDICES ARE IN NON-DECREASING ORDER. THE ROW INDICES
C         CAN BE RECOVERED FROM THE ARRAY IPNTR.
C
C       NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE PARTITION OF THE COLUMNS OF A. COLUMN JCOL BELONGS
C         TO GROUP NGRP(JCOL).
C
C       MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE
C         NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A.
C
C       MINGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES A LOWER
C         BOUND FOR THE NUMBER OF GROUPS IN ANY CONSISTENT PARTITION
C         OF THE COLUMNS OF A.
C
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS. FOR
C         NORMAL TERMINATION INFO = 1. IF M, N, OR NPAIRS IS NOT
C         POSITIVE OR LIWA IS LESS THAN MAX(M,6*N), THEN INFO = 0.
C         IF THE K-TH ELEMENT OF INDROW IS NOT AN INTEGER BETWEEN
C         1 AND M OR THE K-TH ELEMENT OF INDCOL IS NOT AN INTEGER
C         BETWEEN 1 AND N, THEN INFO = -K.
C
C       IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH LIWA.
C
C       LIWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
C         MAX(M,6*N).
C
C       BWA() is an INTEGER WORK ARRAY OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       MINPACK-SUPPLIED ...ID7EGR,I7DO,N7MSRT,M7SEQ,IS7ETR,M7SLO,I7RTDT
C
C       FORTRAN-SUPPLIED ... max
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER I,IR,J,JP,JPL,JPU,K,MAXCLQ,NNZ,NUMGRP
C
C     CHECK THE INPUT DATA.
C
      INFO = 0
      IF (M .LT. 1 .OR. N .LT. 1 .OR. NPAIRS .LT. 1 .OR.
     *    LIWA .LT. max(M,6*N)) GO TO 130
      DO 10 K = 1, NPAIRS
         INFO = -K
         IF (INDROW(K) .LT. 1 .OR. INDROW(K) .GT. M .OR.
     *       INDCOL(K) .LT. 1 .OR. INDCOL(K) .GT. N) GO TO 130
   10    CONTINUE
      INFO = 1
C
C     SORT THE DATA STRUCTURE BY COLUMNS.
C
      CALL I7RTDT(N,NPAIRS,INDROW,INDCOL,JPNTR,IWA(1))
C
C     COMPRESS THE DATA AND DETERMINE THE NUMBER OF
C     NON-ZERO ELEMENTS OF A.
C
      DO 20 I = 1, M
         IWA(I) = 0
   20    CONTINUE
      NNZ = 0
      DO 70 J = 1, N
         JPL = JPNTR(J)
         JPU = JPNTR(J+1) - 1
         JPNTR(J) = NNZ + 1
         IF (JPU .LT. JPL) GO TO 60
         DO 40 JP = JPL, JPU
            IR = INDROW(JP)
            IF (IWA(IR) .NE. 0) GO TO 30
            NNZ = NNZ + 1
            INDROW(NNZ) = IR
            IWA(IR) = 1
   30       CONTINUE
   40       CONTINUE
         JPL = JPNTR(J)
         DO 50 JP = JPL, NNZ
            IR = INDROW(JP)
            IWA(IR) = 0
   50       CONTINUE
   60    CONTINUE
   70    CONTINUE
      JPNTR(N+1) = NNZ + 1
C
C     EXTEND THE DATA STRUCTURE TO ROWS.
C
      CALL IS7ETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA(1))
C
C     DETERMINE A LOWER BOUND FOR THE NUMBER OF GROUPS.
C
      MINGRP = 0
      DO 80 I = 1, M
         MINGRP = max(MINGRP,IPNTR(I+1)-IPNTR(I))
   80    CONTINUE
C
C     DETERMINE THE DEGREE SEQUENCE FOR THE INTERSECTION
C     GRAPH OF THE COLUMNS OF A.
C
      CALL ID7EGR(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(N+1),BWA)
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE SMALLEST-LAST (SL) ORDERING.
C
      CALL M7SLO(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1),
     *         MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1),BWA)
      CALL M7SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),NGRP,MAXGRP,
     *         IWA(N+1),BWA)
      MINGRP = max(MINGRP,MAXCLQ)
      IF (MAXGRP .EQ. MINGRP) GO TO 130
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE INCIDENCE-DEGREE (ID) ORDERING.
C
      CALL I7DO(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1),
     *         MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1),BWA)
      CALL M7SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1),NUMGRP,
     *         IWA(N+1),BWA)
      MINGRP = max(MINGRP,MAXCLQ)
      IF (NUMGRP .GE. MAXGRP) GO TO 100
      MAXGRP = NUMGRP
      DO 90 J = 1, N
         NGRP(J) = IWA(J)
   90    CONTINUE
      IF (MAXGRP .EQ. MINGRP) GO TO 130
  100 CONTINUE
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE LARGEST-FIRST (LF) ORDERING.
C
      CALL N7MSRT(N,N-1,IWA(5*N+1),-1,IWA(4*N+1),IWA(2*N+1),IWA(N+1))
      CALL M7SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1),NUMGRP,
     *         IWA(N+1),BWA)
      IF (NUMGRP .GE. MAXGRP) GO TO 120
      MAXGRP = NUMGRP
      DO 110 J = 1, N
         NGRP(J) = IWA(J)
  110    CONTINUE
  120 CONTINUE
C
C     EXIT FROM PROGRAM.
C
  130 CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE IDSM.
C
      END
      SUBROUTINE I7RTDT(N,NNZ,INDROW,INDCOL,JPNTR,IWA)
      INTEGER N,NNZ
      INTEGER INDROW(NNZ),INDCOL(NNZ),JPNTR(*),IWA(N)
C     **********
C
C     SUBROUTINE I7RTDT
C
C     GIVEN THE NON-ZERO ELEMENTS OF AN M BY N MATRIX A IN
C     ARBITRARY ORDER AS SPECIFIED BY THEIR ROW AND COLUMN
C     INDICES, THIS SUBROUTINE PERMUTES THESE ELEMENTS SO
C     THAT THEIR COLUMN INDICES ARE IN NON-DECREASING ORDER.
C
C     ON INPUT IT IS ASSUMED THAT THE ELEMENTS ARE SPECIFIED IN
C
C           INDROW(K),INDCOL(K), K = 1,...,NNZ.
C
C     ON OUTPUT THE ELEMENTS ARE PERMUTED SO THAT INDCOL IS
C     IN NON-DECREASING ORDER. IN ADDITION, THE ARRAY JPNTR
C     IS SET SO THAT THE ROW INDICES FOR COLUMN J ARE
C
C           INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY I7RTDT AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE I7RTDT(N,NNZ,INDROW,INDCOL,JPNTR,IWA)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       NNZ IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF NON-ZERO ELEMENTS OF A.
C
C       INDROW IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDROW
C         MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS OF A.
C         ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING
C         COLUMN INDICES OF INDCOL ARE IN NON-DECREASING ORDER.
C
C       INDCOL IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDCOL
C         MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS
C         OF A. ON OUTPUT INDCOL IS PERMUTED SO THAT THESE INDICES
C         ARE IN NON-DECREASING ORDER.
C
C       JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN THE OUTPUT
C         INDROW. THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(1) IS SET TO 1 AND THAT JPNTR(N+1)-1
C         IS THEN NNZ.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... MAX
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER I,J,K,L
C
C     DETERMINE THE NUMBER OF NON-ZEROES IN THE COLUMNS.
C
      DO 10 J = 1, N
         IWA(J) = 0
   10    CONTINUE
      DO 20 K = 1, NNZ
         J = INDCOL(K)
         IWA(J) = IWA(J) + 1
   20    CONTINUE
C
C     SET POINTERS TO THE START OF THE COLUMNS IN INDROW.
C
      JPNTR(1) = 1
      DO 30 J = 1, N
         JPNTR(J+1) = JPNTR(J) + IWA(J)
         IWA(J) = JPNTR(J)
   30    CONTINUE
      K = 1
C
C     BEGIN IN-PLACE SORT.
C
   40 CONTINUE
         J = INDCOL(K)
         IF (K .LT. JPNTR(J) .OR. K .GE. JPNTR(J+1)) GO TO 50
C
C           CURRENT ELEMENT IS IN POSITION. NOW EXAMINE THE
C           NEXT ELEMENT OR THE FIRST UN-SORTED ELEMENT IN
C           THE J-TH GROUP.
C
            K = max(K+1,IWA(J))
            GO TO 60
   50    CONTINUE
C
C           CURRENT ELEMENT IS NOT IN POSITION. PLACE ELEMENT
C           IN POSITION AND MAKE THE DISPLACED ELEMENT THE
C           CURRENT ELEMENT.
C
            L = IWA(J)
            IWA(J) = IWA(J) + 1
            I = INDROW(K)
            INDROW(K) = INDROW(L)
            INDCOL(K) = INDCOL(L)
            INDROW(L) = I
            INDCOL(L) = J
   60    CONTINUE
         IF (K .LE. NNZ) GO TO 40
      RETURN
C
C     LAST CARD OF SUBROUTINE I7RTDT.
C
      END
      SUBROUTINE IS7ETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA)
      INTEGER M,N
      INTEGER INDROW(*),JPNTR(*),INDCOL(*),IPNTR(*),IWA(M)
C     **********
C
C     SUBROUTINE IS7ETR
C
C     GIVEN A COLUMN-ORIENTED DEFINITION OF THE SPARSITY PATTERN
C     OF AN M BY N MATRIX A, THIS SUBROUTINE DETERMINES A
C     ROW-ORIENTED DEFINITION OF THE SPARSITY PATTERN OF A.
C
C     ON INPUT THE COLUMN-ORIENTED DEFINITION IS SPECIFIED BY
C     THE ARRAYS INDROW AND JPNTR. ON OUTPUT THE ROW-ORIENTED
C     DEFINITION IS SPECIFIED BY THE ARRAYS INDCOL AND IPNTR.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE IS7ETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF A.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER OUTPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(1) IS SET TO 1 AND THAT IPNTR(M+1)-1 IS
C         THEN THE NUMBER OF NON-ZERO ELEMENTS OF THE MATRIX A.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH M.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER IR,JCOL,JP,JPL,JPU,L,NNZ
C
C     DETERMINE THE NUMBER OF NON-ZEROES IN THE ROWS.
C
      DO 10 IR = 1, M
         IWA(IR) = 0
   10    CONTINUE
      NNZ = JPNTR(N+1) - 1
      DO 20 JP = 1, NNZ
         IR = INDROW(JP)
         IWA(IR) = IWA(IR) + 1
   20    CONTINUE
C
C     SET POINTERS TO THE START OF THE ROWS IN INDCOL.
C
      IPNTR(1) = 1
      DO 30 IR = 1, M
         IPNTR(IR+1) = IPNTR(IR) + IWA(IR)
         IWA(IR) = IPNTR(IR)
   30    CONTINUE
C
C     FILL INDCOL.
C
      DO 60 JCOL = 1, N
         JPL = JPNTR(JCOL)
         JPU = JPNTR(JCOL+1) - 1
         IF (JPU .LT. JPL) GO TO 50
         DO 40 JP = JPL, JPU
            IR = INDROW(JP)
            L = IWA(IR)
            INDCOL(L) = JCOL
            IWA(IR) = IWA(IR) + 1
   40       CONTINUE
   50    CONTINUE
   60    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE IS7ETR.
C
      END
      SUBROUTINE ID7EGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA,BWA)
      INTEGER N
      INTEGER INDROW(*),JPNTR(*),INDCOL(*),IPNTR(*),NDEG(N),IWA(N)
      INTEGER BWA(N), ITRUE, IFALSE
      parameter(ITRUE = 1, IFALSE = 0)
C     **********
C
C     SUBROUTINE ID7EGR
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A,
C     THIS SUBROUTINE DETERMINES THE DEGREE SEQUENCE FOR
C     THE INTERSECTION GRAPH OF THE COLUMNS OF A.
C
C     IN GRAPH-THEORY TERMINOLOGY, THE INTERSECTION GRAPH OF
C     THE COLUMNS OF A IS THE LOOPLESS GRAPH G WITH VERTICES
C     A(J), J = 1,2,...,N WHERE A(J) IS THE J-TH COLUMN OF A
C     AND WITH EDGE (A(I),A(J)) IF AND ONLY IF COLUMNS I AND J
C     HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY ID7EGR AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE ID7EGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA,BWA)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       NDEG IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH
C         SPECIFIES THE DEGREE SEQUENCE. THE DEGREE OF THE
C         J-TH COLUMN OF A IS NDEG(J).
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH N.
C
C       BWA() is an INTEGER WORK ARRAY OF LENGTH N.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER DEG,IC,IP,IPL,IPU,IR,JCOL,JP,JPL,JPU
C
C     INITIALIZATION BLOCK.
C
      DO 10 JP = 1, N
         NDEG(JP) = 0
         BWA(JP) = IFALSE
   10    CONTINUE
C
C     COMPUTE THE DEGREE SEQUENCE BY DETERMINING THE CONTRIBUTIONS
C     TO THE DEGREES FROM THE CURRENT(JCOL) COLUMN AND FURTHER
C     COLUMNS WHICH HAVE NOT YET BEEN CONSIDERED.
C
      IF (N .LT. 2) GO TO 90
      DO 80 JCOL = 2, N
         BWA(JCOL) = ITRUE
         DEG = 0
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         JPL = JPNTR(JCOL)
         JPU = JPNTR(JCOL+1) - 1
         IF (JPU .LT. JPL) GO TO 50
         DO 40 JP = JPL, JPU
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            IPL = IPNTR(IR)
            IPU = IPNTR(IR+1) - 1
            DO 30 IP = IPL, IPU
               IC = INDCOL(IP)
C
C            Array BWA() MARKS COLUMNS WHICH HAVE CONTRIBUTED TO
C              THE DEGREE COUNT OF COLUMN JCOL. UPDATE THE DEGREE
C              COUNTS OF THESE COLUMNS. ARRAY IWA RECORDS THE
C              MARKED COLUMNS.
C
               IF (BWA(IC) .eq. ITRUE) GO TO 20
               BWA(IC) = ITRUE
               NDEG(IC) = NDEG(IC) + 1
               DEG = DEG + 1
               IWA(DEG) = IC
   20          CONTINUE
   30          CONTINUE
   40       CONTINUE
   50    CONTINUE
C
C        UN-MARK THE COLUMNS RECORDED BY IWA AND FINALIZE THE
C        DEGREE COUNT OF COLUMN JCOL.
C
         IF (DEG .LT. 1) GO TO 70
         DO 60 JP = 1, DEG
            IC = IWA(JP)
            BWA(IC) = IFALSE
   60       CONTINUE
         NDEG(JCOL) = NDEG(JCOL) + DEG
   70    CONTINUE
   80    CONTINUE
   90 CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE ID7EGR.
C
      END
      SUBROUTINE M7SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
     *               MAXCLQ,IWA1,IWA2,IWA3,IWA4,BWA)
      INTEGER N,MAXCLQ
      INTEGER INDROW(*),JPNTR(*),INDCOL(*),IPNTR(*),NDEG(N)
      INTEGER LIST(N),IWA1(N),IWA2(N),IWA3(N),IWA4(N)
      INTEGER BWA(N), ITRUE, IFALSE
      parameter(ITRUE = 1, IFALSE = 0)
C     **********
C
C     SUBROUTINE M7SLO
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS
C     SUBROUTINE DETERMINES THE SMALLEST-LAST ORDERING OF THE
C     COLUMNS OF A.
C
C     THE SMALLEST-LAST ORDERING IS DEFINED FOR THE LOOPLESS
C     GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE
C     J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF
C     COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     THE SMALLEST-LAST ORDERING IS DETERMINED RECURSIVELY BY
C     LETTING LIST(K), K = N,...,1 BE A COLUMN WITH LEAST DEGREE
C     IN THE SUBGRAPH SPANNED BY THE UN-ORDERED COLUMNS.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY M7SLO AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE M7SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
C                      MAXCLQ,IWA1,IWA2,IWA3,IWA4,BWA)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN
C         OF A IS NDEG(J).
C
C       LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE SMALLEST-LAST ORDERING OF THE COLUMNS OF A. THE J-TH
C         COLUMN IN THIS ORDER IS LIST(J).
C
C       MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE
C         OF THE LARGEST CLIQUE FOUND DURING THE ORDERING.
C
C       IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N.
C
C       BWA() is an INTEGER WORK ARRAY OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... min
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER DEG,HEAD,IC,IP,IPL,IPU,IR,JCOL,JP,JPL,JPU,
     *        L,MINDEG,NUMDEG,NUMORD
C
C     INITIALIZATION BLOCK.
C
      MINDEG = N
      DO 10 JP = 1, N
         IWA1(JP) = 0
         BWA(JP) = IFALSE
         LIST(JP) = NDEG(JP)
         MINDEG = min(MINDEG,NDEG(JP))
   10    CONTINUE
C
C     CREATE A DOUBLY-LINKED LIST TO ACCESS THE DEGREES OF THE
C     COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS.
C
C     EACH UN-ORDERED COLUMN JCOL IS IN A LIST (THE DEGREE
C     LIST) OF COLUMNS WITH THE SAME DEGREE.
C
C     IWA1(NUMDEG+1) IS THE FIRST COLUMN IN THE NUMDEG LIST
C     UNLESS IWA1(NUMDEG+1) = 0. IN THIS CASE THERE ARE
C     NO COLUMNS IN THE NUMDEG LIST.
C
C     IWA2(JCOL) IS THE COLUMN BEFORE JCOL IN THE DEGREE LIST
C     UNLESS IWA2(JCOL) = 0. IN THIS CASE JCOL IS THE FIRST
C     COLUMN IN THIS DEGREE LIST.
C
C     IWA3(JCOL) IS THE COLUMN AFTER JCOL IN THE DEGREE LIST
C     UNLESS IWA3(JCOL) = 0. IN THIS CASE JCOL IS THE LAST
C     COLUMN IN THIS DEGREE LIST.
C
C     IF JCOL IS AN UN-ORDERED COLUMN, THEN LIST(JCOL) IS THE
C     DEGREE OF JCOL IN THE GRAPH INDUCED BY THE UN-ORDERED
C     COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL)
C     IS THE SMALLEST-LAST ORDER OF COLUMN JCOL.
C
      DO 20 JP = 1, N
         NUMDEG = NDEG(JP)
         HEAD = IWA1(NUMDEG+1)
         IWA1(NUMDEG+1) = JP
         IWA2(JP) = 0
         IWA3(JP) = HEAD
         IF (HEAD .GT. 0) IWA2(HEAD) = JP
   20    CONTINUE
      MAXCLQ = 0
      NUMORD = N
C
C     BEGINNING OF ITERATION LOOP.
C
   30 CONTINUE
C
C        MARK THE SIZE OF THE LARGEST CLIQUE
C        FOUND DURING THE ORDERING.
C
         IF (MINDEG + 1 .EQ. NUMORD .AND. MAXCLQ .EQ. 0)
     *       MAXCLQ = NUMORD
C
C        CHOOSE A COLUMN JCOL OF MINIMAL DEGREE MINDEG.
C
   40    CONTINUE
            JCOL = IWA1(MINDEG+1)
            IF (JCOL .GT. 0) GO TO 50
            MINDEG = MINDEG + 1
            GO TO 40
   50    CONTINUE
         LIST(JCOL) = NUMORD
         NUMORD = NUMORD - 1
C
C        TERMINATION TEST.
C
         IF (NUMORD .EQ. 0) GO TO 120
C
C        DELETE COLUMN JCOL FROM THE MINDEG LIST.
C
         L = IWA3(JCOL)
         IWA1(MINDEG+1) = L
         IF (L .GT. 0) IWA2(L) = 0
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
         BWA(JCOL) = ITRUE
         DEG = 0
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         JPL = JPNTR(JCOL)
         JPU = JPNTR(JCOL+1) - 1
         IF (JPU .LT. JPL) GO TO 90
         DO 80 JP = JPL, JPU
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            IPL = IPNTR(IR)
            IPU = IPNTR(IR+1) - 1
            DO 70 IP = IPL, IPU
               IC = INDCOL(IP)
C
C            Array BWA() MARKS COLUMNS WHICH ARE ADJACENT TO
C              COLUMN JCOL. ARRAY IWA4 RECORDS THE MARKED COLUMNS.
C
               IF (BWA(IC) .eq. ITRUE) GO TO 60
               BWA(IC) = ITRUE
               DEG = DEG + 1
               IWA4(DEG) = IC
   60          CONTINUE
   70          CONTINUE
   80       CONTINUE
   90    CONTINUE
C
C        UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS.
C
         IF (DEG .LT. 1) GO TO 110
         DO 100 JP = 1, DEG
            IC = IWA4(JP)
            NUMDEG = LIST(IC)
            LIST(IC) = LIST(IC) - 1
            MINDEG = min(MINDEG,LIST(IC))
C
C           DELETE COLUMN IC FROM THE NUMDEG LIST.
C
            L = IWA2(IC)
            IF (L .EQ. 0) IWA1(NUMDEG+1) = IWA3(IC)
            IF (L .GT. 0) IWA3(L) = IWA3(IC)
            L = IWA3(IC)
            IF (L .GT. 0) IWA2(L) = IWA2(IC)
C
C           ADD COLUMN IC TO THE NUMDEG-1 LIST.
C
            HEAD = IWA1(NUMDEG)
            IWA1(NUMDEG) = IC
            IWA2(IC) = 0
            IWA3(IC) = HEAD
            IF (HEAD .GT. 0) IWA2(HEAD) = IC
C
C           UN-MARK COLUMN IC IN THE ARRAY BWA.
C
            BWA(IC) = IFALSE
  100       CONTINUE
  110    CONTINUE
C
C        END OF ITERATION LOOP.
C
         GO TO 30
  120 CONTINUE
C
C     INVERT THE ARRAY LIST.
C
      DO 130 JCOL = 1, N
         NUMORD = LIST(JCOL)
         IWA1(NUMORD) = JCOL
  130    CONTINUE
      DO 140 JP = 1, N
         LIST(JP) = IWA1(JP)
  140    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE M7SLO.
C
      END
      SUBROUTINE M7SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP,
     *               IWA,BWA)
      INTEGER N,MAXGRP
      INTEGER INDROW(*),JPNTR(*),INDCOL(*),IPNTR(*),LIST(N),NGRP(N)
      INTEGER IWA(N)
      INTEGER BWA(N), ITRUE, IFALSE
      parameter(ITRUE = 1, IFALSE = 0)
C     **********
C
C     SUBROUTINE M7SEQ
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS
C     SUBROUTINE DETERMINES A CONSISTENT PARTITION OF THE
C     COLUMNS OF A BY A SEQUENTIAL ALGORITHM.
C
C     A CONSISTENT PARTITION IS DEFINED IN TERMS OF THE LOOPLESS
C     GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE
C     J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF
C     COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     A PARTITION OF THE COLUMNS OF A INTO GROUPS IS CONSISTENT
C     IF THE COLUMNS IN ANY GROUP ARE NOT ADJACENT IN THE GRAPH G.
C     IN GRAPH-THEORY TERMINOLOGY, A CONSISTENT PARTITION OF THE
C     COLUMNS OF A CORRESPONDS TO A COLORING OF THE GRAPH G.
C
C     THE SUBROUTINE EXAMINES THE COLUMNS IN THE ORDER SPECIFIED
C     BY THE ARRAY LIST, AND ASSIGNS THE CURRENT COLUMN TO THE
C     GROUP WITH THE SMALLEST POSSIBLE NUMBER.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY M7SEQ AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE M7SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP,
C                      IWA,BWA)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       LIST IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE ORDER TO BE USED BY THE SEQUENTIAL ALGORITHM.
C         THE J-TH COLUMN IN THIS ORDER IS LIST(J).
C
C       NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE PARTITION OF THE COLUMNS OF A. COLUMN JCOL BELONGS
C         TO GROUP NGRP(JCOL).
C
C       MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE
C         NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH N.
C
C       BWA() is an INTEGER WORK ARRAY OF LENGTH N.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER DEG,IC,IP,IPL,IPU,IR,J,JCOL,JP,JPL,JPU,L,NUMGRP
C
C     INITIALIZATION BLOCK.
C
      MAXGRP = 0
      DO 10 JP = 1, N
         NGRP(JP) = N
         BWA(JP) = IFALSE
   10    CONTINUE
      BWA(N) = ITRUE
C
C     BEGINNING OF ITERATION LOOP.
C
      DO 100 J = 1, N
         JCOL = LIST(J)
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
         DEG = 0
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         JPL = JPNTR(JCOL)
         JPU = JPNTR(JCOL+1) - 1
         IF (JPU .LT. JPL) GO TO 50
         DO 40 JP = JPL, JPU
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            IPL = IPNTR(IR)
            IPU = IPNTR(IR+1) - 1
            DO 30 IP = IPL, IPU
               IC = INDCOL(IP)
               L = NGRP(IC)
C
C            Array BWA() MARKS THE GROUP NUMBERS OF THE
C              COLUMNS WHICH ARE ADJACENT TO COLUMN JCOL.
C              ARRAY IWA RECORDS THE MARKED GROUP NUMBERS.
C
               IF (BWA(L) .eq. ITRUE) GO TO 20
               BWA(L) = ITRUE
               DEG = DEG + 1
               IWA(DEG) = L
   20          CONTINUE
   30          CONTINUE
   40       CONTINUE
   50    CONTINUE
C
C        ASSIGN THE SMALLEST UN-MARKED GROUP NUMBER TO JCOL.
C
         DO 60 JP = 1, N
            NUMGRP = JP
            IF (BWA(JP) .eq. IFALSE) GO TO 70
   60       CONTINUE
   70    CONTINUE
         NGRP(JCOL) = NUMGRP
         MAXGRP = max(MAXGRP,NUMGRP)
C
C        UN-MARK THE GROUP NUMBERS.
C
         IF (DEG .LT. 1) GO TO 90
         DO 80 JP = 1, DEG
            L = IWA(JP)
            BWA(L) = IFALSE
   80       CONTINUE
   90    CONTINUE
  100    CONTINUE
C
C        END OF ITERATION LOOP.
C
      RETURN
C
C     LAST CARD OF SUBROUTINE M7SEQ.
C
      END
      SUBROUTINE I7DO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
     *               MAXCLQ,IWA1,IWA2,IWA3,IWA4,BWA)
      INTEGER M,N,MAXCLQ
      INTEGER INDROW(*),JPNTR(*),INDCOL(*),IPNTR(*),NDEG(N),LIST(N)
      INTEGER IWA1(N),IWA2(N),IWA3(N),IWA4(N)
      INTEGER BWA(N), ITRUE, IFALSE
      parameter(ITRUE = 1, IFALSE = 0)
C     **********
C
C     SUBROUTINE I7DO
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS
C     SUBROUTINE DETERMINES AN INCIDENCE-DEGREE ORDERING OF THE
C     COLUMNS OF A.
C
C     THE INCIDENCE-DEGREE ORDERING IS DEFINED FOR THE LOOPLESS
C     GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE
C     J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF
C     COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     AT EACH STAGE OF I7DO, A COLUMN OF MAXIMAL INCIDENCE IS
C     CHOSEN AND ORDERED. IF JCOL IS AN UN-ORDERED COLUMN, THEN
C     THE INCIDENCE OF JCOL IS THE NUMBER OF ORDERED COLUMNS
C     ADJACENT TO JCOL IN THE GRAPH G. AMONG ALL THE COLUMNS OF
C     MAXIMAL INCIDENCE,I7DO CHOOSES A COLUMN OF MAXIMAL DEGREE.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE I7DO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
C                      MAXCLQ,IWA1,IWA2,IWA3,IWA4,BWA)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF A.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN
C         OF A IS NDEG(J).
C
C       LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE INCIDENCE-DEGREE ORDERING OF THE COLUMNS OF A. THE J-TH
C         COLUMN IN THIS ORDER IS LIST(J).
C
C       MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE
C         OF THE LARGEST CLIQUE FOUND DURING THE ORDERING.
C
C       IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N.
C
C       BWA() is an INTEGER WORK ARRAY OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       MINPACK-SUPPLIED ... N7MSRT
C
C       FORTRAN-SUPPLIED ... max
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER DEG,HEAD,IC,IP,IPL,IPU,IR,JCOL,JP,JPL,JPU,L,MAXINC,
     *        MAXLST,NCOMP,NUMINC,NUMLST,NUMORD,NUMWGT
C
C     SORT THE DEGREE SEQUENCE.
C
      CALL N7MSRT(N,N-1,NDEG,-1,IWA4,IWA1,IWA3)
C
C     INITIALIZATION BLOCK.
C
C     CREATE A DOUBLY-LINKED LIST TO ACCESS THE INCIDENCES OF THE
C     COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS.
C
C     EACH UN-ORDERED COLUMN JCOL IS IN A LIST (THE INCIDENCE LIST)
C     OF COLUMNS WITH THE SAME INCIDENCE.
C
C     IWA1(NUMINC+1) IS THE FIRST COLUMN IN THE NUMINC LIST
C     UNLESS IWA1(NUMINC+1) = 0. IN THIS CASE THERE ARE
C     NO COLUMNS IN THE NUMINC LIST.
C
C     IWA2(JCOL) IS THE COLUMN BEFORE JCOL IN THE INCIDENCE LIST
C     UNLESS IWA2(JCOL) = 0. IN THIS CASE JCOL IS THE FIRST
C     COLUMN IN THIS INCIDENCE LIST.
C
C     IWA3(JCOL) IS THE COLUMN AFTER JCOL IN THE INCIDENCE LIST
C     UNLESS IWA3(JCOL) = 0. IN THIS CASE JCOL IS THE LAST
C     COLUMN IN THIS INCIDENCE LIST.
C
C     IF JCOL IS AN UN-ORDERED COLUMN, THEN LIST(JCOL) IS THE
C     INCIDENCE OF JCOL IN THE GRAPH. IF JCOL IS AN ORDERED COLUMN,
C     THEN LIST(JCOL) IS THE INCIDENCE-DEGREE ORDER OF COLUMN JCOL.
C
      MAXINC = 0
      DO 10 JP = 1, N
         LIST(JP) = 0
         BWA(JP) = IFALSE
         IWA1(JP) = 0
         L = IWA4(JP)
         IF (JP .NE. 1) IWA2(L) = IWA4(JP-1)
         IF (JP .NE. N) IWA3(L) = IWA4(JP+1)
   10    CONTINUE
      IWA1(1) = IWA4(1)
      L = IWA4(1)
      IWA2(L) = 0
      L = IWA4(N)
      IWA3(L) = 0
C
C     DETERMINE THE MAXIMAL SEARCH LENGTH FOR THE LIST
C     OF COLUMNS OF MAXIMAL INCIDENCE.
C
      MAXLST = 0
      DO 20 IR = 1, M
         MAXLST = MAXLST + (IPNTR(IR+1) - IPNTR(IR))**2
   20    CONTINUE
      MAXLST = MAXLST/N
      MAXCLQ = 1
C
C     BEGINNING OF ITERATION LOOP.
C
      DO 140 NUMORD = 1, N
C
C        CHOOSE A COLUMN JCOL OF MAXIMAL DEGREE AMONG THE
C        COLUMNS OF MAXIMAL INCIDENCE.
C
         JP = IWA1(MAXINC+1)
         NUMLST = 1
         NUMWGT = -1
   30    CONTINUE
            IF (NDEG(JP) .LE. NUMWGT) GO TO 40
            NUMWGT = NDEG(JP)
            JCOL = JP
   40       CONTINUE
            JP = IWA3(JP)
            NUMLST = NUMLST + 1
            IF (JP .GT. 0 .AND. NUMLST .LE. MAXLST) GO TO 30
         LIST(JCOL) = NUMORD
C
C        DELETE COLUMN JCOL FROM THE LIST OF COLUMNS OF
C        MAXIMAL INCIDENCE.
C
         L = IWA2(JCOL)
         IF (L .EQ. 0) IWA1(MAXINC+1) = IWA3(JCOL)
         IF (L .GT. 0) IWA3(L) = IWA3(JCOL)
         L = IWA3(JCOL)
         IF (L .GT. 0) IWA2(L) = IWA2(JCOL)
C
C        UPDATE THE SIZE OF THE LARGEST CLIQUE
C        FOUND DURING THE ORDERING.
C
         IF (MAXINC .EQ. 0) NCOMP = 0
         NCOMP = NCOMP + 1
         IF (MAXINC + 1 .EQ. NCOMP) MAXCLQ = max(MAXCLQ,NCOMP)
C
C        UPDATE THE MAXIMAL INCIDENCE COUNT.
C
   50    CONTINUE
            IF (IWA1(MAXINC+1) .GT. 0) GO TO 60
            MAXINC = MAXINC - 1
            IF (MAXINC .GE. 0) GO TO 50
   60    CONTINUE
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
         BWA(JCOL) = ITRUE
         DEG = 0
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         JPL = JPNTR(JCOL)
         JPU = JPNTR(JCOL+1) - 1
         IF (JPU .LT. JPL) GO TO 100
         DO 90 JP = JPL, JPU
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            IPL = IPNTR(IR)
            IPU = IPNTR(IR+1) - 1
            DO 80 IP = IPL, IPU
               IC = INDCOL(IP)
C
C            Array BWA() MARKS COLUMNS WHICH ARE ADJACENT TO
C              COLUMN JCOL. ARRAY IWA4 RECORDS THE MARKED COLUMNS.
C
               IF (BWA(IC) .eq. ITRUE) GO TO 70
               BWA(IC) = ITRUE
               DEG = DEG + 1
               IWA4(DEG) = IC
   70          CONTINUE
   80          CONTINUE
   90       CONTINUE
  100    CONTINUE
C
C        UPDATE THE POINTERS TO THE INCIDENCE LISTS.
C
         IF (DEG .LT. 1) GO TO 130
         DO 120 JP = 1, DEG
            IC = IWA4(JP)
            IF (LIST(IC) .GT. 0) GO TO 110
            NUMINC = -LIST(IC) + 1
            LIST(IC) = -NUMINC
            MAXINC = max(MAXINC,NUMINC)
C
C           DELETE COLUMN IC FROM THE NUMINC-1 LIST.
C
            L = IWA2(IC)
            IF (L .EQ. 0) IWA1(NUMINC) = IWA3(IC)
            IF (L .GT. 0) IWA3(L) = IWA3(IC)
            L = IWA3(IC)
            IF (L .GT. 0) IWA2(L) = IWA2(IC)
C
C           ADD COLUMN IC TO THE NUMINC LIST.
C
            HEAD = IWA1(NUMINC+1)
            IWA1(NUMINC+1) = IC
            IWA2(IC) = 0
            IWA3(IC) = HEAD
            IF (HEAD .GT. 0) IWA2(HEAD) = IC
  110       CONTINUE
C
C           UN-MARK COLUMN IC IN THE ARRAY BWA.
C
            BWA(IC) = IFALSE
  120       CONTINUE
  130    CONTINUE
         BWA(JCOL) = IFALSE
C
C        END OF ITERATION LOOP.
C
  140    CONTINUE
C
C     INVERT THE ARRAY LIST.
C
      DO 150 JCOL = 1, N
         NUMORD = LIST(JCOL)
         IWA1(NUMORD) = JCOL
  150    CONTINUE
      DO 160 JP = 1, N
         LIST(JP) = IWA1(JP)
  160    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE I7DO.
C
      END
      SUBROUTINE N7MSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT)
      INTEGER N,NMAX,MODE
      INTEGER NUM(N),INDEX(N),LAST(*),NEXT(N)
C     **********.
C
C     SUBROUTINE N7MSRT
C
C     GIVEN A SEQUENCE OF INTEGERS, THIS SUBROUTINE GROUPS
C     TOGETHER THOSE INDICES WITH THE SAME SEQUENCE VALUE
C     AND, OPTIONALLY, SORTS THE SEQUENCE INTO EITHER
C     ASCENDING OR DESCENDING ORDER.
C
C     THE SEQUENCE OF INTEGERS IS DEFINED BY THE ARRAY NUM,
C     AND IT IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET
C     0,1,...,NMAX. ON OUTPUT THE INDICES K SUCH THAT NUM(K) = L
C     FOR ANY L = 0,1,...,NMAX CAN BE OBTAINED FROM THE ARRAYS
C     LAST AND NEXT AS FOLLOWS.
C
C           K = LAST(L+1)
C           WHILE (K .NE. 0) K = NEXT(K)
C
C     OPTIONALLY, THE SUBROUTINE PRODUCES AN ARRAY INDEX SO THAT
C     THE SEQUENCE NUM(INDEX(I)), I = 1,2,...,N IS SORTED.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE N7MSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE.
C
C       NMAX IS A POSITIVE INTEGER INPUT VARIABLE.
C
C       NUM IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE
C         SEQUENCE OF INTEGERS TO BE GROUPED AND SORTED. IT
C         IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET
C         0,1,...,NMAX.
C
C       MODE IS AN INTEGER INPUT VARIABLE. THE SEQUENCE NUM IS
C         SORTED IN ASCENDING ORDER IF MODE IS POSITIVE AND IN
C         DESCENDING ORDER IF MODE IS NEGATIVE. IF MODE IS 0,
C         NO SORTING IS DONE.
C
C       INDEX IS AN INTEGER OUTPUT ARRAY OF LENGTH N SET SO
C         THAT THE SEQUENCE
C
C               NUM(INDEX(I)), I = 1,2,...,N
C
C         IS SORTED ACCORDING TO THE SETTING OF MODE. IF MODE
C         IS 0, INDEX IS NOT REFERENCED.
C
C       LAST IS AN INTEGER OUTPUT ARRAY OF LENGTH NMAX + 1. THE
C         INDEX OF NUM FOR THE LAST OCCURRENCE OF L IS LAST(L+1)
C         FOR ANY L = 0,1,...,NMAX UNLESS LAST(L+1) = 0. IN
C         THIS CASE L DOES NOT APPEAR IN NUM.
C
C       NEXT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IF
C         NUM(K) = L, THEN THE INDEX OF NUM FOR THE PREVIOUS
C         OCCURRENCE OF L IS NEXT(K) FOR ANY L = 0,1,...,NMAX
C         UNLESS NEXT(K) = 0. IN THIS CASE THERE IS NO PREVIOUS
C         OCCURRENCE OF L IN NUM.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1982.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE
C
C     **********
      INTEGER I,J,JP,K,L,NMAXP1,NMAXP2
C
C     DETERMINE THE ARRAYS NEXT AND LAST.
C
      NMAXP1 = NMAX + 1
      DO 10 I = 1, NMAXP1
         LAST(I) = 0
   10    CONTINUE
      DO 20 K = 1, N
         L = NUM(K)
         NEXT(K) = LAST(L+1)
         LAST(L+1) = K
   20    CONTINUE
      IF (MODE .EQ. 0) GO TO 60
C
C     STORE THE POINTERS TO THE SORTED ARRAY IN INDEX.
C
      I = 1
      NMAXP2 = NMAXP1 + 1
      DO 50 J = 1, NMAXP1
         JP = J
         IF (MODE .LT. 0) JP = NMAXP2 - J
         K = LAST(JP)
   30    CONTINUE
            IF (K .EQ. 0) GO TO 40
            INDEX(I) = K
            I = I + 1
            K = NEXT(K)
            GO TO 30
   40    CONTINUE
   50    CONTINUE
   60 CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE N7MSRT.
C
      END