C     ------------------------------------------------------------
C
C*  Title
C
C     Supplementary package implementing numerical differentation
C     on a given sparse pattern - for use in connection with NLEQ1S.
C
C*  Written by        U. Nowak, L. Weimann 
C*  Purpose           Evaluation of Jacobian by numerical differentation 
C                     for use with NLEQ1S
C*  Category          F2a. - Systems of nonlinear equations
C*  Version           1.0
C*  Revision          August 2003
C*  Latest Change     August 2003
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Fortran 77 environment on PC's,
C                     workstations and hosts.
C*  Copyright     (c) Konrad-Zuse-Zentrum fuer
C                     Informationstechnik Berlin (ZIB)
C                     Takustrasse 7, D-14195 Berlin-Dahlem
C                     phone : + 49/30/84185-0
C                     fax   : + 49/30/84185-125
C*  Contact           Bodo Erdmann
C                     ZIB, Division Scientific Computing, 
C                          Department Numerical Analysis and Modelling
C                     phone : + 49/30/84185-185
C                     fax   : + 49/30/84185-107
C                     e-mail: erdmann@zib.de
C
C     ------------------------------------------------------------
C
C     Description:
C
C     This package supplies a replacement routine for the 
C     analytical Jacobian, named NDJAC, which is required by NLEQ1S. 
C     This routine either provides the Jacobian by numerical 
C     differentation on a sparse pattern defined by row- and 
C     column-indices which may be either provided by a user-written
C     routine named JACPAT, or automatically determined by full mode
C     or band mode numerical differentation. In connection 
C     with the last mentioned option, the user must supply just
C     a dummy routine JACPAT.
C
C     Usage:
C
C     To use this package, follow these steps:
C
C     1) Search this file for a comment line starting with "c-->"
C        and adapt the constants nmax and nzmax to fit the
C        requirements of your problem function and its Jacobian.
C        
C     2) Include the following lines into your main program, 
C        which calls NLEQ1S:
C
C     program main
C     ...
C     integer inistate, iJacPat, ml, mu
C     common / jacinf / ml, mu, iJacPat, inistate
C     external ndjac
C     ...
C     inistate = 0
Cc    insert the appropriate numbers below for the lower and
Cc    upper bandwidth of the Jacobian (or set both to n for 
Cc    a full mode Jacobian).
C     ml = <lower bandwidth of Jacobian>
C     mu = <upper bandwidth of Jacobian>
Cc    note that the problem function must be named 'fcn'.
Cc
Cc    If the sparse pattern should be automatically determined, insert
C     iJacPat = 0
Cc    or, if the sparse pattern should be supplied by the routine JacPat,
Cc    then insert
C     iJacPat = 1
C     call nleq1s(n,nzmax,fcn,ndjac,x,xscal,rtol,iopt,ierr,liwk,iwk,
C    $            li2wk,i2wk,lrwk,rwk)
C     ...
C     (For a description of the nleq1s call above refer to the nleq1s
C     documentation.)
C
C     3) Provide a subroutine JACPAT according to these specifications:
C
C     subroutine jacpat(n, x, irow, icol, lenjac, info)
C
C     n       integer  in   Number of equations and unknowns
C     x       double   in   starting point of Newton iteration
C     irow(*) integer  out  Array of row indices ir of the pattern
C                           pairs (ir,ic)
C     icol(*) integer  out  Array of column indices ic of the pattern
C                           pairs (ir,ic)
C     lenjac  integer  in   Max_Non_Zeros_Jacobian
C                      out  Number of nonzero elements of the Jacobian,
C                           .i.e. number of pairs (ir,ic).
C     info    integer  out  set to 0 if jacpat evaluation is successful;
C                           set to a negative number to indicate an
C                           error return from routine jacpat
C                   
C     ------------------------------------------------------------
C
      subroutine NDJAC (n, x, a, ir, ic, lena, info)
c
      implicit none
c
      integer nzmax, nmax
c-->  adjust parameters in the code line below according to your needs
c     nmax  : Maximum number of Equations
c     nzmax : Maximum number of nonzero Jacobian elements
      parameter( nmax = 1000, nzmax = 10000 )
c
      integer n, lena, info, ir(*), ic(*)
      double precision x(n), a(*)
c
      integer i, info1, info2, lenas
c
      integer ml, mu, iJacPat, inistate
c
      common / jacinf / ml, mu, iJacPat, inistate
c
      save lenas
c
      integer lrwk, nrw, l1, l2, l3, l4, l5,
     2        liwk, niw, k1, k2, k3, k4, k5, k6, k7
c
      parameter( lrwk = 4*nmax,
     2           liwk = 2*nzmax+9*nmax+2 )
      double precision rwk(lrwk)
      integer iwk(liwk)
      integer icnt
      save icnt
c
c  generation of sparse pattern and numerical differentiation on it 
c
      if (n.gt.nmax) stop 'NDJAC: nmax too small'
c     nzmax = lena
c
      l1 = 1
      l2 = l1+n
      l3 = l2+n
      l4 = l3+n
      l5 = l4+n
      nrw = l5-1
c
      k1 = 1
      k2 = k1+nzmax
      k3 = k2+nzmax
      k4 = k3+n+1
      k5 = k4+n+1
      k6 = k5+n
      k7 = k6+6*n
      niw = k7-1
c
      info=0
      info1=0
      info2=0
      if (inistate.eq.0) icnt = 1
c
c  get the pattern of Jacobian 
c
      if( iJacPat .eq. 1 ) then
c  Get Jacobian sparse pattern
         if ( icnt .eq. 1 ) then ! evaluate at very first call only
            call JacPat(n, x, iwk(k1), iwk(k2), lenas, info1)
         endif
c
      else
c
c  full mode numerical differentation and dropping very small terms
         if(inistate.eq.0) then
            if((ml+mu+1) .lt. n) then
               call getjpatb 
     1              (n, x, a, iwk(k1), iwk(k2), lenas,
     2               rwk(l1), rwk(l2), rwk(l3), rwk(l4), info1)
            else
               call getjpat (n, x, a, iwk(k1), iwk(k2), lenas,
     2                       rwk(l1), rwk(l2), rwk(l3), info1)
            endif
            if(info1.lt.0) then
                info=info1
                return
            endif
c           info2=1
            inistate=1
         endif
c
      endif
c
c  call evaluation routine
c
      call ndjpat (n, x, a, iwk(k1), iwk(k2), lenas, 
     2             rwk(l1), rwk(l2), rwk(l3), rwk(l4),
     3             iwk(k3), iwk(k4), iwk(k5), iwk(k6), info2)
c
      if(info2.lt.0) then
         info=info2
         return
      endif
c
      if(lenas.gt.nzmax) then
         info=-99
         return
      endif
c
      lena=lenas
      do i=1, lena
         ir(i)= iwk(k1-1+i)
         ic(i)= iwk(k2-1+i)
      enddo
c
      icnt = icnt+1
c
      return
c
c----------------------------------------------------------------------
c  end subroutine NDJAC
c----------------------------------------------------------------------
c
      end
      subroutine ndjpat (n, x, a, ir, ic, lena, xh, dh, f0, f1, 
     2                   ipntr, jpntr, ngrp, iwa, info )
c
      implicit none
c
c----------------------------------------------------------------------
c  do numerical differentiation using a given sparse pattern
c----------------------------------------------------------------------
c
      integer n, lena, info, ir(*), ic(*)
      double precision x(n), a(*), xh(n), dh(n), f0(n), f1(n)
      integer ipntr(n+1), jpntr(n+1), ngrp(n), iwa(6*n)
c
      double precision xpert
c
c
      integer i, info1, irow, j, jcol, jp, liwa, lumon,
     2 maxgrp, mingrp, maxrow, minrow, m, ml, mu, numgrp, nnz
c
      logical linfoj
c
      parameter (lumon = 6)
c
      linfoj=.false.
      if(info.gt.0) linfoj=.true.
      info=0
c
      call fcn (n, x, f0, info)
c
      m=n
      nnz=lena 
      liwa=6*n 
c
      if(linfoj) write(lumon,*) ' '
      if(linfoj) write(lumon,*) ' Info on sparse pattern:'
c
      call dsm (m, n, nnz, ir, ic, ngrp, maxgrp, mingrp,
     *          info1, ipntr, jpntr, iwa, liwa)
c
      if(info1 .le. 0) then
         write(lumon,*) ' ERROR return from dsm, info =',info1
         if (info1.lt.0) then
           write(lumon,*) 'ir(',-info1,')=',ir(-info1)
           write(lumon,*) 'ic(',-info1,')=',ic(-info1)
         endif
         info=-1
         return
      endif
c
      if(linfoj) then
         write(lumon,*) '  nnz old and new:',lena,nnz
         maxrow = 0
         minrow = n
         do i = 1, n
            maxrow = max(maxrow,ipntr(i+1)-ipntr(i))
            minrow = min(minrow,ipntr(i+1)-ipntr(i))
         enddo
         write(lumon,*) '  row nnz min and max:',minrow,maxrow
         write(lumon,*) '  col groups min and act:',mingrp,maxgrp
         ml=0
         mu=0
         do jcol=1,n
            do jp = jpntr(jcol), jpntr(jcol+1)-1
                  irow = ir(jp)
                  ml=max(ml,irow-jcol)
                  mu=max(mu,jcol-irow)
            enddo 
         enddo
         write(lumon,*) '  max ml and mu:', ml, mu 
      endif
c
c  evaluate jacobian
c
      do numgrp=1,maxgrp
         do j=1,n
            dh(j) = 0.0
            if (ngrp(j) .eq. numgrp) then
               xpert=max(abs(x(j))*1.d-7,1.d-10)! TRSH evtl. adaptieren!
               dh(j)=xpert
            endif
            xh(j) = x(j) + dh(j)
         enddo    
c
         call fcn (n, xh, f1, info)
c
         do jcol=1,n
            if (ngrp(jcol) .eq. numgrp) then
               do jp = jpntr(jcol), jpntr(jcol+1)-1
                  irow = ir(jp)
                  a(jp) = (f1(irow)-f0(irow))/dh(jcol)
                  ic(jp)=jcol
               enddo    
            endif
         enddo   
      enddo   
      lena=nnz
c
c
      return
c
c----------------------------------------------------------------------
c  end subroutine ndjpat
c----------------------------------------------------------------------
c
      end
      subroutine getjpat (n, xin, a, ir, ic, lena, x, f0, f1, info)
c
      implicit none
c
c----------------------------------------------------------------------
c  create sparse pattern by full mode numerical differentation and
c  dropping very small terms
c----------------------------------------------------------------------
c
      integer n, lena, info, ir(*), ic(*)
      double precision xin(n), a(*)
c
      double precision x(n), f0(n), f1(n)
c
      double precision xsave, xpert, aval
c
      integer i, j, k, l
c
c
      info=0
c
      do i=1,n
         x(i)=xin(i)
      enddo
c     
      do i=1,n
         if(abs(x(i)).lt. 1.d-31) x(i)=dble(i)*1.d-4      
c        x(i)= x(i) + max( x(i), 1.d-3 )*dble(i)*1.d-3  
      enddo
c
      call fcn (n, x, f0, info)
      if (info.ne.0) return
c
      l=0
c
      do k=1,n
c
         xsave=x(k)
         xpert=max(abs(x(k))*1.d-4,1.d-10)
         x(k)=x(k)+xpert
c
         call fcn (n, x, f1, info)
         if (info.ne.0) return
c
         do j=1,n
            aval=(f1(j)-f0(j))/xpert
            if(abs(aval).gt.1.d-50) then
               l=l+1
               a(l)=aval
               ir(l)=j
               ic(l)=k
            endif
         enddo
c
         x(k)=xsave
c
      enddo
c
      lena=l
c
      return
c
c----------------------------------------------------------------------
c  end subroutine getjpat
c----------------------------------------------------------------------
c
      end
      subroutine getjpatb (n, xin, a, ir, ic, lena, x, f0, f1, xh,
     2                     info)
c
      implicit none
c
c----------------------------------------------------------------------
c  create sparse pattern by band mode numerical differentation and
c  dropping very small terms
c----------------------------------------------------------------------
c
      integer n, lena, info, ir(*), ic(*)
      double precision xin(n), a(*)
      double precision x(n), f0(n), f1(n), xh(n)
c
      double precision xpert, aval
c
      integer i, iupper, ilower, j, k, l, mb
c
      integer ml, mu, iJacPat, inistate
c
      common / jacinf / ml, mu, iJacPat, inistate
c
c

      info=0
c
c
      do i=1,n
         x(i)=xin(i)
      enddo
c
      do i=1,n
         if(abs(x(i)).lt. 1.d-31) x(i)=dble(i)*1.d-4      
      enddo
c
      call fcn (n, x, f0, info)
      if (info.ne.0) return
c
      l=0
c
      mb = ml+mu+1
c
      do k=1,mb
c
         do j=1,n
           xh(j)=x(j)
         enddo   
c
         do j=k,n,mb
            xpert=max(abs(x(j))*1.d-4,1.d-10)
            xh(j)=xh(j)+xpert
         enddo
c
         call fcn (n, xh, f1, info)
         if (info.ne.0) return
c
         do j=k,n,mb
c
            iupper= max(1, j-mu)
            ilower= min(n, j+ml)
c
            do i= iupper, ilower
               aval=( f1(i)-f0(i) ) / ( xh(j)-x(j) )
               if(abs(aval).gt.1.d-25) then
                  l=l+1
                  a(l)=aval
                  ir(l)=i
                  ic(l)=j
               endif
            enddo
         enddo
c
      enddo
c
      lena=l
c
c
      return
c
c  end subroutine getjpatb
c
      end
c
      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
