* ============================================================================ *
      SUBROUTINE DSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP,
     *               INFO,IPNTR,JPNTR,IWA,LIWA)
      INTEGER M,N,NPAIRS,MAXGRP,MINGRP,INFO,LIWA
      INTEGER INDROW(NPAIRS),INDCOL(NPAIRS),NGRP(N),
     *        IPNTR(M+1),JPNTR(N+1),IWA(LIWA)
C     **********
C
C     SUBROUTINE DSM
C
C     THE PURPOSE OF DSM 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 DSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP,
C                      INFO,IPNTR,JPNTR,IWA,LIWA)
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     SUBPROGRAMS CALLED
C
C       MINPACK-SUPPLIED ... DEGR,IDO,NUMSRT,SEQ,SETR,SLO,SRTDAT
C
C       FORTRAN-SUPPLIED ... MAX
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER I,IR,J,JP,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)) RETURN
      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) RETURN
   10    CONTINUE
      INFO = 1
C
C     SORT THE DATA STRUCTURE BY COLUMNS.
C
      CALL SRTDAT(N,NPAIRS,INDROW,INDCOL,JPNTR,IWA)
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 = 1
      DO 40 J = 1, N
         K = NNZ
         DO 30 JP = JPNTR(J), JPNTR(J+1)-1
            IR = INDROW(JP)
            IF (IWA(IR) .NE. J) THEN
               INDROW(NNZ) = IR
               NNZ = NNZ + 1
               IWA(IR) = J
               END IF
   30       CONTINUE
         JPNTR(J) = K
   40    CONTINUE
      JPNTR(N+1) = NNZ
C
C     EXTEND THE DATA STRUCTURE TO ROWS.
C
      CALL SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA)
C
C     DETERMINE A LOWER BOUND FOR THE NUMBER OF GROUPS.
C
      MINGRP = 0
      DO 50 I = 1, M
         MINGRP = MAX(MINGRP,IPNTR(I+1)-IPNTR(I))
   50    CONTINUE
C
C     DETERMINE THE DEGREE SEQUENCE FOR THE INTERSECTION
C     GRAPH OF THE COLUMNS OF A.
C
      CALL DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(N+1))
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE SMALLEST-LAST (SL) ORDERING.
C
      CALL SLO(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))
      CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),NGRP,MAXGRP,
     *         IWA(N+1))
      MINGRP = MAX(MINGRP,MAXCLQ)
C
C     EXIT IF THE SMALLEST-LAST ORDERING IS OPTIMAL.
C
      IF (MAXGRP .EQ. MINGRP) RETURN
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE INCIDENCE-DEGREE (ID) ORDERING.
C
      CALL IDO(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))
      CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1),NUMGRP,
     *         IWA(N+1))
      MINGRP = MAX(MINGRP,MAXCLQ)
C
C     RETAIN THE BETTER OF THE TWO ORDERINGS SO FAR.
C
      IF (NUMGRP .LT. MAXGRP) THEN
         MAXGRP = NUMGRP
         DO 60 J = 1, N
            NGRP(J) = IWA(J)
   60       CONTINUE
C
C        EXIT IF THE INCIDENCE-DEGREE ORDERING IS OPTIMAL.
C
         IF (MAXGRP .EQ. MINGRP) RETURN
         END IF
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE LARGEST-FIRST (LF) ORDERING.
C
      CALL NUMSRT(N,N-1,IWA(5*N+1),-1,IWA(4*N+1),IWA(2*N+1),IWA(N+1))
      CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1),NUMGRP,
     *         IWA(N+1))
C
C     RETAIN THE BEST OF THE THREE ORDERINGS AND EXIT.
C
      IF (NUMGRP .LT. MAXGRP) THEN
         MAXGRP = NUMGRP
         DO 70 J = 1, N
            NGRP(J) = IWA(J)
   70       CONTINUE
         END IF
      RETURN
C
C     LAST CARD OF SUBROUTINE DSM.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA)
      INTEGER N
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N),IWA(N)
C     **********
C
C     SUBROUTINE DEGR
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 DEGR AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA)
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     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,JCOL,JP
C
C     INITIALIZATION BLOCK.
C
      DO 10 JP = 1, N
         NDEG(JP) = 0
         IWA(JP) = 0
   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
      DO 40 JCOL = 2, N
         IWA(JCOL) = N
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA MARKS COLUMNS WHICH HAVE CONTRIBUTED TO
C              THE DEGREE COUNT OF COLUMN JCOL. UPDATE THE DEGREE
C              COUNTS OF THESE COLUMNS AS WELL AS COLUMN JCOL.
C
               IF (IWA(IC) .LT. JCOL) THEN
                  IWA(IC) = JCOL
                  NDEG(IC) = NDEG(IC) + 1
                  NDEG(JCOL) = NDEG(JCOL) + 1
                  END IF
   20          CONTINUE
   30       CONTINUE
   40    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE DEGR.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
     *               MAXCLQ,IWA1,IWA2,IWA3,IWA4)
      INTEGER M,N,MAXCLQ
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),NDEG(N),
     *        LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N)
C     **********
C
C     SUBROUTINE IDO
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     THE INCIDENCE-DEGREE ORDERING IS DETERMINED RECURSIVELY BY
C     LETTING LIST(K), K = 1,...,N BE A COLUMN WITH MAXIMAL
C     INCIDENCE TO THE SUBGRAPH SPANNED BY THE ORDERED COLUMNS.
C     AMONG ALL THE COLUMNS OF MAXIMAL INCIDENCE, IDO CHOOSES A
C     COLUMN OF MAXIMAL DEGREE.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
C                      MAXCLQ,IWA1,IWA2,IWA3,IWA4)
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     SUBPROGRAMS CALLED
C
C       MINPACK-SUPPLIED ... NUMSRT
C
C       FORTRAN-SUPPLIED ... MAX
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,JCOL,JP,
     *        MAXINC,MAXLST,NCOMP,NUMINC,NUMLST,NUMORD,NUMWGT
C
C     SORT THE DEGREE SEQUENCE.
C
      CALL NUMSRT(N,N-1,NDEG,-1,IWA4,IWA2,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 IC IS IN A LIST (THE INCIDENCE LIST)
C     OF COLUMNS WITH THE SAME INCIDENCE.
C
C     IWA1(NUMINC) IS THE FIRST COLUMN IN THE NUMINC LIST
C     UNLESS IWA1(NUMINC) = 0. IN THIS CASE THERE ARE
C     NO COLUMNS IN THE NUMINC LIST.
C
C     IWA2(IC) IS THE COLUMN BEFORE IC IN THE INCIDENCE LIST
C     UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST
C     COLUMN IN THIS INCIDENCE LIST.
C
C     IWA3(IC) IS THE COLUMN AFTER IC IN THE INCIDENCE LIST
C     UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST
C     COLUMN IN THIS INCIDENCE LIST.
C
C     IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE
C     INCIDENCE OF IC TO THE GRAPH INDUCED BY THE ORDERED
C     COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL)
C     IS THE INCIDENCE-DEGREE ORDER OF COLUMN JCOL.
C
      MAXINC = 0
      DO 10 JP = N, 1, -1
         IC = IWA4(JP)
         IWA1(N-JP) = 0
         IWA2(IC) = 0
         IWA3(IC) = IWA1(0)
         IF (IWA1(0) .GT. 0) IWA2(IWA1(0)) = IC
         IWA1(0) = IC
         IWA4(JP) = 0
         LIST(JP) = 0
   10    CONTINUE
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 = 0
      NUMORD = 1
C
C     BEGINNING OF ITERATION LOOP.
C
   30 CONTINUE
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        CHOOSE A COLUMN JCOL OF MAXIMAL DEGREE AMONG THE
C        COLUMNS OF MAXIMAL INCIDENCE MAXINC.
C
   40    CONTINUE
            JP = IWA1(MAXINC)
            IF (JP .GT. 0) GO TO 50
            MAXINC = MAXINC - 1
            GO TO 40
   50    CONTINUE
         NUMWGT = -1
         DO 60 NUMLST = 1, MAXLST
            IF (NDEG(JP) .GT. NUMWGT) THEN
               NUMWGT = NDEG(JP)
               JCOL = JP
               END IF
            JP = IWA3(JP)
            IF (JP .LE. 0) GO TO 70
   60       CONTINUE
   70    CONTINUE
         LIST(JCOL) = NUMORD
         NUMORD = NUMORD + 1
C
C        TERMINATION TEST.
C
         IF (NUMORD .GT. N) GO TO 100
C
C        DELETE COLUMN JCOL FROM THE MAXINC LIST.
C
         IF (IWA2(JCOL) .EQ. 0) THEN
            IWA1(MAXINC) = IWA3(JCOL)
         ELSE
            IWA3(IWA2(JCOL)) = IWA3(JCOL)
            END IF
         IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = IWA2(JCOL)
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
         IWA4(JCOL) = N
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 90 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 80 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO
C              COLUMN JCOL.
C
               IF (IWA4(IC) .LT. NUMORD) THEN
                  IWA4(IC) = NUMORD
C
C                 UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS.
C
                  NUMINC = LIST(IC)
                  LIST(IC) = LIST(IC) + 1
                  MAXINC = MAX(MAXINC,LIST(IC))
C
C                 DELETE COLUMN IC FROM THE NUMINC LIST.
C
                  IF (IWA2(IC) .EQ. 0) THEN
                     IWA1(NUMINC) = IWA3(IC)
                  ELSE
                     IWA3(IWA2(IC)) = IWA3(IC)
                     END IF
                  IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC)
C
C                 ADD COLUMN IC TO THE NUMINC+1 LIST.
C
                  IWA2(IC) = 0
                  IWA3(IC) = IWA1(NUMINC+1)
                  IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = IC
                  IWA1(NUMINC+1) = IC
                  END IF
   80          CONTINUE
   90       CONTINUE
C
C        END OF ITERATION LOOP.
C
         GO TO 30
  100 CONTINUE
C
C     INVERT THE ARRAY LIST.
C
      DO 110 JCOL = 1, N
         IWA2(LIST(JCOL)) = JCOL
  110    CONTINUE
      DO 120 JP = 1, N
         LIST(JP) = IWA2(JP)
  120    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE IDO.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE NUMSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT)
      INTEGER N,NMAX,MODE
      INTEGER NUM(N),INDEX(N),LAST(0:NMAX),NEXT(N)
C     **********.
C
C     SUBROUTINE NUMSRT
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)
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 NUMSRT(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)
C         FOR ANY L = 0,1,...,NMAX UNLESS LAST(L) = 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. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER I,J,JINC,JL,JU,K,L
C
C     DETERMINE THE ARRAYS NEXT AND LAST.
C
      DO 10 I = 0, NMAX
         LAST(I) = 0
   10    CONTINUE
      DO 20 K = 1, N
         L = NUM(K)
         NEXT(K) = LAST(L)
         LAST(L) = K
   20    CONTINUE
      IF (MODE .EQ. 0) RETURN
C
C     STORE THE POINTERS TO THE SORTED ARRAY IN INDEX.
C
      I = 1
      IF (MODE .GT. 0) THEN
         JL = 0
         JU = NMAX
         JINC = 1
      ELSE
         JL = NMAX
         JU = 0
         JINC = -1
         END IF
      DO 50 J = JL, JU, JINC
         K = LAST(J)
   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
      RETURN
C
C     LAST CARD OF SUBROUTINE NUMSRT.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP,
     *               IWA)
      INTEGER N,MAXGRP
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),LIST(N),
     *        NGRP(N),IWA(N)
C     **********
C
C     SUBROUTINE SEQ
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 SEQ AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP,
C                      IWA)
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     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,J,JCOL,JP
C
C     INITIALIZATION BLOCK.
C
      MAXGRP = 0
      DO 10 JP = 1, N
         NGRP(JP) = N
         IWA(JP) = 0
   10    CONTINUE
C
C     BEGINNING OF ITERATION LOOP.
C
      DO 60 J = 1, N
         JCOL = LIST(J)
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA MARKS THE GROUP NUMBERS OF THE
C              COLUMNS WHICH ARE ADJACENT TO COLUMN JCOL.
C
               IWA(NGRP(IC)) = J
   20          CONTINUE
   30       CONTINUE
C
C        ASSIGN THE SMALLEST UN-MARKED GROUP NUMBER TO JCOL.
C
         DO 40 JP = 1, MAXGRP
            IF (IWA(JP) .NE. J) GO TO 50
   40       CONTINUE
         MAXGRP = MAXGRP + 1
   50    CONTINUE
         NGRP(JCOL) = JP
   60    CONTINUE
C
C        END OF ITERATION LOOP.
C
      RETURN
C
C     LAST CARD OF SUBROUTINE SEQ.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA)
      INTEGER M,N
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),IWA(M)
C     **********
C
C     SUBROUTINE SETR
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 SETR(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. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IR,JCOL,JP
C
C     STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE ROWS.
C
      DO 10 IR = 1, M
         IWA(IR) = 0
   10    CONTINUE
      DO 20 JP = 1, JPNTR(N+1)-1
         IWA(INDROW(JP)) = IWA(INDROW(JP)) + 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 50 JCOL = 1, N
         DO 40 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
            INDCOL(IWA(IR)) = JCOL
            IWA(IR) = IWA(IR) + 1
   40       CONTINUE
   50    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE SETR.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
     *               MAXCLQ,IWA1,IWA2,IWA3,IWA4)
      INTEGER N,MAXCLQ
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N),
     *        LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N)
C     **********
C
C     SUBROUTINE SLO
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 SLO AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
C                      MAXCLQ,IWA1,IWA2,IWA3,IWA4)
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     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... MIN
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,JCOL,JP,MINDEG,NUMDEG,NUMORD
C
C     INITIALIZATION BLOCK.
C
      MINDEG = N
      DO 10 JP = 1, N
         IWA1(JP-1) = 0
         IWA4(JP) = N
         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 IC IS IN A LIST (THE DEGREE LIST)
C     OF COLUMNS WITH THE SAME DEGREE.
C
C     IWA1(NUMDEG) IS THE FIRST COLUMN IN THE NUMDEG LIST
C     UNLESS IWA1(NUMDEG) = 0. IN THIS CASE THERE ARE
C     NO COLUMNS IN THE NUMDEG LIST.
C
C     IWA2(IC) IS THE COLUMN BEFORE IC IN THE DEGREE LIST
C     UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST
C     COLUMN IN THIS DEGREE LIST.
C
C     IWA3(IC) IS THE COLUMN AFTER IC IN THE DEGREE LIST
C     UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST
C     COLUMN IN THIS DEGREE LIST.
C
C     IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE
C     DEGREE OF IC 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)
         IWA2(JP) = 0
         IWA3(JP) = IWA1(NUMDEG)
         IF (IWA1(NUMDEG) .GT. 0) IWA2(IWA1(NUMDEG)) = JP
         IWA1(NUMDEG) = 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)
            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 80
C
C        DELETE COLUMN JCOL FROM THE MINDEG LIST.
C
         IWA1(MINDEG) = IWA3(JCOL)
         IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = 0
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
         IWA4(JCOL) = 0
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 70 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 60 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO
C              COLUMN JCOL.
C
               IF (IWA4(IC) .GT. NUMORD) THEN
                  IWA4(IC) = NUMORD
C
C                 UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS.
C
                  NUMDEG = LIST(IC)
                  LIST(IC) = LIST(IC) - 1
                  MINDEG = MIN(MINDEG,LIST(IC))
C
C                 DELETE COLUMN IC FROM THE NUMDEG LIST.
C
                  IF (IWA2(IC) .EQ. 0) THEN
                     IWA1(NUMDEG) = IWA3(IC)
                  ELSE
                     IWA3(IWA2(IC)) = IWA3(IC)
                     END IF
                  IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC)
C
C                 ADD COLUMN IC TO THE NUMDEG-1 LIST.
C
                  IWA2(IC) = 0
                  IWA3(IC) = IWA1(NUMDEG-1)
                  IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = IC
                  IWA1(NUMDEG-1) = IC
                  END IF
   60          CONTINUE
   70       CONTINUE
C
C        END OF ITERATION LOOP.
C
         GO TO 30
   80 CONTINUE
C
C     INVERT THE ARRAY LIST.
C
      DO 90 JCOL = 1, N
         IWA2(LIST(JCOL)) = JCOL
   90    CONTINUE
      DO 100 JP = 1, N
         LIST(JP) = IWA2(JP)
  100    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE SLO.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA)
      INTEGER N,NNZ
      INTEGER INDROW(NNZ),INDCOL(NNZ),JPNTR(N+1),IWA(N)
C     **********
C
C     SUBROUTINE SRTDAT
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 SRTDAT AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE SRTDAT(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. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER I,J,K,L
C
C     STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE COLUMNS.
C
      DO 10 J = 1, N
         IWA(J) = 0
   10    CONTINUE
      DO 20 K = 1, NNZ
         IWA(INDCOL(K)) = IWA(INDCOL(K)) + 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 .GE. JPNTR(J)) THEN
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))
         ELSE
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
            END IF
         IF (K .LE. NNZ) GO TO 40
      RETURN
C
C     LAST CARD OF SUBROUTINE SRTDAT.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE FDJS(M,N,COL,IND,NPNTR,NGRP,NUMGRP,d,FJACD,FJAC)
      INTEGER M,N,NUMGRP
      INTEGER IND(*),NPNTR(*),NGRP(N)
      REAL D(N),FJACD(M),FJAC(*)
      LOGICAL COL
C     **********
C
C     SUBROUTINE FDJS
C
C     GIVEN A CONSISTENT PARTITION OF THE COLUMNS OF AN M BY N
C     JACOBIAN MATRIX INTO GROUPS, THIS SUBROUTINE COMPUTES
C     APPROXIMATIONS TO THOSE COLUMNS IN A GIVEN GROUP.  THE
C     APPROXIMATIONS ARE STORED INTO EITHER A COLUMN-ORIENTED
C     OR A ROW-ORIENTED PATTERN.
C
C     A PARTITION IS CONSISTENT IF THE COLUMNS IN ANY GROUP
C     DO NOT HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     APPROXIMATIONS TO THE COLUMNS OF THE JACOBIAN MATRIX IN A
C     GIVEN GROUP CAN BE OBTAINED BY SPECIFYING A DIFFERENCE
C     PARAMETER ARRAY D WITH D(JCOL) NON-ZERO IF AND ONLY IF
C     JCOL IS A COLUMN IN THE GROUP, AND AN APPROXIMATION TO
C     JAC*D WHERE JAC DENOTES THE JACOBIAN MATRIX OF A MAPPING F.
C
C     D CAN BE DEFINED WITH THE FOLLOWING SEGMENT OF CODE.
C
C           DO 10 JCOL = 1, N
C              D(JCOL) = 0.0
C              IF (NGRP(JCOL) .EQ. NUMGRP) D(JCOL) = ETA(JCOL)
C        10    CONTINUE
C
C     IN THE ABOVE CODE NUMGRP IS THE GIVEN GROUP NUMBER,
C     NGRP(JCOL) IS THE GROUP NUMBER OF COLUMN JCOL, AND
C     ETA(JCOL) IS THE DIFFERENCE PARAMETER USED TO
C     APPROXIMATE COLUMN JCOL OF THE JACOBIAN MATRIX.
C     SUITABLE VALUES FOR THE ARRAY ETA MUST BE PROVIDED.
C
C     AS MENTIONED ABOVE, AN APPROXIMATION TO JAC*D MUST
C     ALSO BE PROVIDED. FOR EXAMPLE, THE APPROXIMATION
C
C           F(X+D) - F(X)
C
C     CORRESPONDS TO THE FORWARD DIFFERENCE FORMULA AT X.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE FDJS(M,N,COL,IND,NPNTR,NGRP,NUMGRP,D,FJACD,FJAC)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF THE JACOBIAN MATRIX.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF THE JACOBIAN MATRIX.
C
C       COL IS A LOGICAL INPUT VARIABLE. IF COL IS SET TRUE, THEN THE
C         JACOBIAN APPROXIMATIONS ARE STORED INTO A COLUMN-ORIENTED
C         PATTERN. IF COL IS SET FALSE, THEN THE JACOBIAN
C         APPROXIMATIONS ARE STORED INTO A ROW-ORIENTED PATTERN.
C
C       IND IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE JACOBIAN MATRIX
C         IF COL IS TRUE, AND CONTAINS THE COLUMN INDICES FOR
C         THE NON-ZEROES IN THE JACOBIAN MATRIX IF COL IS FALSE.
C
C       NPNTR IS AN INTEGER INPUT ARRAY WHICH SPECIFIES THE
C         LOCATIONS OF THE ROW INDICES IN IND IF COL IS TRUE, AND
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN IND IF
C         COL IS FALSE. IF COL IS TRUE, THE INDICES FOR COLUMN J ARE
C
C               IND(K), K = NPNTR(J),...,NPNTR(J+1)-1.
C
C         IF COL IS FALSE, THE INDICES FOR ROW I ARE
C
C               IND(K), K = NPNTR(I),...,NPNTR(I+1)-1.
C
C         NOTE THAT NPNTR(N+1)-1 IF COL IS TRUE, OR NPNTR(M+1)-1
C         IF COL IS FALSE, IS THEN THE NUMBER OF NON-ZERO ELEMENTS
C         OF THE JACOBIAN MATRIX.
C
C       NGRP IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE PARTITION OF THE COLUMNS OF THE JACOBIAN MATRIX.
C         COLUMN JCOL BELONGS TO GROUP NGRP(JCOL).
C
C       NUMGRP IS A POSITIVE INTEGER INPUT VARIABLE SET TO A GROUP
C         NUMBER IN THE PARTITION. THE COLUMNS OF THE JACOBIAN
C         MATRIX IN THIS GROUP ARE TO BE ESTIMATED ON THIS CALL.
C
C       D IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE
C         DIFFERENCE PARAMETER VECTOR FOR THE ESTIMATE OF
C         THE JACOBIAN MATRIX COLUMNS IN GROUP NUMGRP.
C
C       FJACD IS AN INPUT ARRAY OF LENGTH M WHICH CONTAINS
C         AN APPROXIMATION TO THE DIFFERENCE VECTOR JAC*D,
C         WHERE JAC DENOTES THE JACOBIAN MATRIX.
C
C       FJAC IS AN OUTPUT ARRAY OF LENGTH NNZ, WHERE NNZ IS THE
C         NUMBER OF ITS NON-ZERO ELEMENTS. AT EACH CALL OF FDJS,
C         FJAC IS UPDATED TO INCLUDE THE NON-ZERO ELEMENTS OF THE
C         JACOBIAN MATRIX FOR THOSE COLUMNS IN GROUP NUMGRP. FJAC
C         SHOULD NOT BE ALTERED BETWEEN SUCCESSIVE CALLS TO FDJS.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IP,IROW,JCOL,JP
C
C     COMPUTE ESTIMATES OF JACOBIAN MATRIX COLUMNS IN GROUP
C     NUMGRP. THE ARRAY FJACD MUST CONTAIN AN APPROXIMATION
C     TO JAC*D, WHERE JAC DENOTES THE JACOBIAN MATRIX AND D
C     IS A DIFFERENCE PARAMETER VECTOR WITH D(JCOL) NON-ZERO
C     IF AND ONLY IF JCOL IS A COLUMN IN GROUP NUMGRP.
c  M,N,COL,IND,NPNTR,NGRP,NUMGRP,D,FJACD,FJAC
C
      IF (COL) THEN
C
C        COLUMN ORIENTATION.
C
         DO 20 JCOL = 1, N
            IF (NGRP(JCOL) .EQ. NUMGRP) THEN
               DO 10 JP = NPNTR(JCOL), NPNTR(JCOL+1)-1
                  IROW = IND(JP)
                  FJAC(JP) = FJACD(IROW)/D(JCOL)
   10             CONTINUE
               END IF
   20       CONTINUE
      ELSE
C
C        ROW ORIENTATION.
C
         DO 50 IROW = 1, M
            DO 30 IP = NPNTR(IROW), NPNTR(IROW+1)-1
               JCOL = IND(IP)
               IF (NGRP(JCOL) .EQ. NUMGRP) THEN
c	          write (81,'(3i5, es16.8)') irow, ip, jcol, d(jcol)
                  FJAC(IP) = FJACD(IROW)/D(JCOL)
                  GO TO 40
                  END IF
   30          CONTINUE
   40       CONTINUE
   50       CONTINUE
         END IF
      RETURN
C
C     LAST CARD OF SUBROUTINE FDJS.
C
      END
* ============================================================================ *
