coloring.f90

牛顿优化算法源fortran代码

!    FORTRAN-supplied ... max

!  Argonne National Laboratory. MINPACK Project. August 1984.
!  Thomas F. Coleman, Burton S. Garbow, Jorge J. More'

!  **********

INTEGER :: ic, ip, ir, jcol, jp, maxinc, maxlst, ncomp, numinc, numlst,  &
           numord, numwgt, iwa1(0:n-1), iwa2(n), iwa3(n), iwa4(n)

!     Sort the degree sequence.

CALL numsrt(n, n-1, ndeg, -1, iwa4, iwa2, iwa3)

!     Initialization block.

!     Create a doubly-linked list to access the incidences of the
!     columns.  The pointers for the linked list are as follows.

!     Each un-ordered column ic is in a list (the incidence list)
!     of columns with the same incidence.

!     iwa1(numinc) is the first column in the numinc list
!     unless iwa1(numinc) = 0.  In this case there are
!     no columns in the numinc list.

!     iwa2(ic) is the column before ic in the incidence list unless
!     iwa2(ic) = 0.
!     In this case ic is the first column in this incidence list.

!     iwa3(ic) is the column after ic in the incidence list unless iwa3(ic) = 0.
!     In this case ic is the last column in this incidence list.

!     If ic is an un-ordered column, then list(ic) is the incidence of ic to
!     the graph induced by the ordered columns.  If jcol is an ordered column,
!     then list(jcol) is the incidence-degree order of column jcol.

maxinc = 0
DO  jp = n, 1, -1
  ic = iwa4(jp)
  iwa1(n-jp) = 0
  iwa2(ic) = 0
  iwa3(ic) = iwa1(0)
  IF (iwa1(0) > 0) iwa2(iwa1(0)) = ic
  iwa1(0) = ic
  iwa4(jp) = 0
  list(jp) = 0
END DO

!     Determine the maximal search length for the list
!     of columns of maximal incidence.

maxlst = 0
DO  ir = 1, m
  maxlst = maxlst + (ipntr(ir+1) - ipntr(ir))**2
END DO
maxlst = maxlst/n
maxclq = 0
numord = 1

!     Beginning of iteration loop.

!        Choose a column jcol of maximal degree among the
!        columns of maximal incidence maxinc.

30 DO
  jp = iwa1(maxinc)
  IF (jp > 0) EXIT
  maxinc = maxinc - 1
END DO
numwgt = -1
DO  numlst = 1, maxlst
  IF (ndeg(jp) > numwgt) THEN
    numwgt = ndeg(jp)
    jcol = jp
  END IF
  jp = iwa3(jp)
  IF (jp <= 0) EXIT
END DO
list(jcol) = numord

!        Update the size of the largest clique found during the ordering.

IF (maxinc == 0) ncomp = 0
ncomp = ncomp + 1
IF (maxinc+1 == ncomp) maxclq = MAX(maxclq,ncomp)

!        Termination test.

numord = numord + 1
IF (numord > n) GO TO 100

!        Delete column jcol from the maxinc list.

IF (iwa2(jcol) == 0) THEN
  iwa1(maxinc) = iwa3(jcol)
ELSE
  iwa3(iwa2(jcol)) = iwa3(jcol)
END IF
IF (iwa3(jcol) > 0) iwa2(iwa3(jcol)) = iwa2(jcol)

!        Find all columns adjacent to column jcol.

iwa4(jcol) = n

!        Determine all positions (ir,jcol) which correspond
!        to non-zeroes in the matrix.

DO  jp = jpntr(jcol), jpntr(jcol+1)-1
  ir = indrow(jp)
  
!           For each row ir, determine all positions (ir,ic)
!           which correspond to non-zeroes in the matrix.
  
  DO  ip = ipntr(ir), ipntr(ir+1)-1
    ic = indcol(ip)
    
!              Array iwa4 marks columns which are adjacent to column jcol.
    
    IF (iwa4(ic) < numord) THEN
      iwa4(ic) = numord
      
!                 Update the pointers to the current incidence lists.
      
      numinc = list(ic)
      list(ic) = list(ic) + 1
      maxinc = MAX(maxinc,list(ic))
      
!                 Delete column ic from the numinc list.
      
      IF (iwa2(ic) == 0) THEN
        iwa1(numinc) = iwa3(ic)
      ELSE
        iwa3(iwa2(ic)) = iwa3(ic)
      END IF
      IF (iwa3(ic) > 0) iwa2(iwa3(ic)) = iwa2(ic)
      
!                 Add column ic to the numinc+1 list.
      
      iwa2(ic) = 0
      iwa3(ic) = iwa1(numinc+1)
      IF (iwa1(numinc+1) > 0) iwa2(iwa1(numinc+1)) = ic
      iwa1(numinc+1) = ic
    END IF
  END DO
END DO

!        End of iteration loop.

GO TO 30

!     Invert the array list.

100 DO  jcol = 1, n
  iwa2(list(jcol)) = jcol
END DO
list(1:n) = iwa2(1:n)
RETURN

!     Last card of subroutine ido.

END SUBROUTINE ido



SUBROUTINE idog(n, nghbrp, npntrp, nghbrs, npntrs, listp, maxclq, maxid)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:49:04

INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(IN)   :: nghbrp(:)
INTEGER, INTENT(IN)   :: npntrp(:)    ! npntrp(n+1)
INTEGER, INTENT(IN)   :: nghbrs(:)
INTEGER, INTENT(IN)   :: npntrs(:)    ! npntrs(n+1)
INTEGER, INTENT(OUT)  :: listp(:)
INTEGER, INTENT(OUT)  :: maxclq
INTEGER, INTENT(OUT)  :: maxid

!  **********

!  subroutine idog

!  Given a loopless graph G = (V,E), this subroutine determines
!  the incidence degree ordering of the vertices of G.

!  The incidence degree ordering is determined recursively by
!  letting list(k), k = 1,...,n be a vertex with maximal
!  incidence to the subgraph spanned by the ordered vertices.
!  Among all the vertices of maximal incidence, a vertex of
!  maximal degree is chosen. This subroutine determines the
!  inverse of the incidence degree ordering, that is, an array
!  listp such that listp(list(k)) = k for k = 1,2,...,n.

!  The subroutine statement is

!    subroutine idog(n, nghbrp, npntrp, nghbrs, npntrs, listp,
!                    maxclq, maxid, iwa1, iwa2, iwa3)

!  where

!    n is a positive integer input variable set to the number of vertices of G.

!    nghbrp is an integer input array which contains the
!      predecessor adjacency lists for the graph G.

!    npntrp is an integer input array of length n + 1 which specifies the
!      locations of the predecessor adjacency lists in nghbrp.
!      The vertices preceding and adjacent to vertex j are

!            nghbrp(k), k = npntrp(j),...,npntrp(j+1)-1.

!      Note that npntrp(n+1)-1 is then the number of vertices
!      plus edges of the graph G.

!    nghbrs is an integer input array which contains the
!      successor adjacency lists for the graph G.

!    npntrs is an integer input array of length n + 1 which specifies the
!      locations of the successor adjacency lists in nghbrs.
!      The vertices succeeding and adjacent to vertex j are

!            nghbrs(k), k = npntrs(j),...,npntrs(j+1)-1.

!      Note that npntrs(n+1)-1 is then the number of vertices
!      plus edges of the graph G.

!    listp is an integer output array of length n which specifies
!      the inverse of the incidence degree ordering of the
!      vertices.  Vertex j is in position listp(j) of this ordering.

!    maxclq is an integer output variable set to the size
!      of the largest clique found during the ordering.

!    maxid is an integer output variable set to the maximum
!      incidence degree found during the ordering.

!    iwa1, iwa2, and iwa3 are integer work arrays of length n.

!  Subprograms called

!    MINPACK-supplied ... numsrt

!    FORTRAN-supplied ... max

!  Argonne National Laboratory. MINPACK Project. December 1984.
!  Thomas F. Coleman, Burton S. Garbow, Jorge J. More'

!  **********

INTEGER :: i, j, k, maxinc, maxdeg, maxlst, ncomp, numdeg, numinc, numord,  &
           iwa1(0:n-1), iwa2(n), iwa3(n)

!     Initialization block.

DO  j = 1, n
  listp(j) = (npntrp(j+1) - npntrp(j) - 1) + (npntrs(j+1) - npntrs(j) - 1)
END DO
maxlst = (npntrp(n+1) + npntrs(n+1))/n

!     Sort the degree sequence.

CALL numsrt(n, n-1, listp, 1, iwa1, iwa2, iwa3)

!     Create a doubly-linked list to access the incidences of the
!     vertices. The pointers for the linked list are as follows.

!     Each un-ordered vertex i is in a list (the incidence list)
!     of vertices with the same incidence.

!     iwa1(numinc) is the first vertex in the numinc list unless
!     iwa1(numinc) = 0.  In this case there are no vertices in the numinc list.

!     iwa2(i) is the vertex before i in the incidence list unless iwa2(i) = 0.
!     In this case i is the first vertex in this incidence list.

!     iwa3(i) is the vertex after i in the incidence list unless iwa3(i) = 0.
!     In this case i is the last vertex in this incidence list.

!     If i is an un-ordered vertex, then -listp(i) is the incidence of i to the
!     graph induced by the ordered vertices.  If j is an ordered vertex, then
!     listp(j) is the incidence degree order of vertex j.

maxinc = 0
DO  j = 1, n
  i = iwa1(j-1)
  iwa1(j-1) = 0
  iwa2(i) = 0
  iwa3(i) = iwa1(0)
  IF (iwa1(0) > 0) iwa2(iwa1(0)) = i
  iwa1(0) = i
  listp(j) = 0
END DO
maxclq = 0
maxid = 0
numord = 1

!     Beginning of iteration loop.

!        Choose a vertex j of maximal degree among the
!        vertices of maximal incidence maxinc.

30 DO
  k = iwa1(maxinc)
  IF (k > 0) EXIT
  maxinc = maxinc - 1
END DO

maxdeg = -1
DO  i = 1, maxlst
  numdeg = (npntrp(k+1) - npntrp(k) - 1) + (npntrs(k+1) - npntrs(k) - 1)
  IF (numdeg > maxdeg) THEN
    maxdeg = numdeg
    j = k
  END IF
  k = iwa3(k)
  IF (k <= 0) EXIT
END DO

listp(j) = numord
maxid = MAX(maxid, maxinc)

!        Update the size of the largest clique found during the ordering.

IF (maxinc == 0) ncomp = 0
ncomp = ncomp + 1
IF (maxinc+1 == ncomp) maxclq = MAX(maxclq, ncomp)

!        Termination test.

numord = numord + 1
IF (numord > n) GO TO 100

!        Delete vertex j from the maxinc list.

IF (iwa2(j) == 0) THEN
  iwa1(maxinc) = iwa3(j)
ELSE
  iwa3(iwa2(j)) = iwa3(j)
END IF
IF (iwa3(j) > 0) iwa2(iwa3(j)) = iwa2(j)

!        Determine all the neighbors of vertex j which precede j
!        in the subgraph spanned by the un-ordered vertices.

DO  k = npntrp(j), npntrp(j+1)-1
  i = nghbrp(k)
  
!           Update the pointers to the current incidence lists.
  
  numinc = -listp(i)
  IF (numinc >= 0) THEN
    listp(i) = listp(i) - 1
    maxinc = MAX(maxinc,-listp(i))
    
!              Delete vertex i from the numinc list.
    
    IF (iwa2(i) == 0) THEN
      iwa1(numinc) = iwa3(i)
    ELSE
      iwa3(iwa2(i)) = iwa3(i)
    END IF
    IF (iwa3(i) > 0) iwa2(iwa3(i)) = iwa2(i)
    
!              Add vertex i to the numinc+1 list.
    
    iwa2(i) = 0
    iwa3(i) = iwa1(numinc+1)
    IF (iwa1(numinc+1) > 0) iwa2(iwa1(numinc+1)) = i
    iwa1(numinc+1) = i
  END IF
END DO

!        Determine all the neighbors of vertex j which succeed j
!        in the subgraph spanned by the un-ordered vertices.

DO  k = npntrs(j), npntrs(j+1)-1
  i = nghbrs(k)
  
!           Update the pointers to the current incidence lists.
  
  numinc = -listp(i)
  IF (numinc >= 0) THEN
    listp(i) = listp(i) - 1
    maxinc = MAX(maxinc,-listp(i))
    
!              Delete vertex i from the numinc list.
    
    IF (iwa2(i) == 0) THEN
      iwa1(numinc) = iwa3(i)
    ELSE
      iwa3(iwa2(i)) = iwa3(i)
    END IF
    IF (iwa3(i) > 0) iwa2(iwa3(i)) = iwa2(i)
    
!              Add vertex i to the numinc+1 list.
    
    iwa2(i) = 0
    iwa3(i) = iwa1(numinc+1)
    IF (iwa1(numinc+1) > 0) iwa2(iwa1(numinc+1)) = i
    iwa1(numinc+1) = i
  END IF
END DO

!        End of iteration loop.

GO TO 30
100 RETURN

!     Last card of subroutine idog.

END SUBROUTINE idog



SUBROUTINE numsrt(n, nmax, num, mode, INDEX, last, next)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:49:10

INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(IN)   :: nmax
INTEGER, INTENT(IN)   :: num(:)
INTEGER, INTENT(IN)   :: mode
INTEGER, INTENT(OUT)  :: INDEX(:)
INTEGER, INTENT(OUT)  :: last(0:)
INTEGER, INTENT(OUT)  :: next(:)

!  **********.

!  subroutine numsrt

!  Given a sequence of integers, this subroutine groups together those indices
!  with the same sequence value and, optionally, sorts the sequence into either
!  ascending or descending order.

!  The sequence of integers is defined by the array num, and it is assumed
!  that the integers are each from the set 0,1,...,nmax.  On output the indices
!  k such that num(k) = l for any l = 0,1,...,nmax can be obtained from the
!  arrays last and next as follows.

!        k = last(l)
!        while (k .ne. 0) k = next(k)

!  Optionally, the subroutine produces an array index so that
!  the sequence num(index(i)), i = 1,2,...,n is sorted.

!  The subroutine statement is

!    subroutine numsrt(n, nmax, num, mode, index, last, next)

!  where

!    n is a positive integer input variable.

!    nmax is a positive integer input variable.

!    num is an input array of length n which contains the sequence of
!      integers to be grouped and sorted.  It is assumed that the integers
!      are each from the set 0,1,...,nmax.

!    mode is an integer input variable.  The sequence num is sorted in
!      ascending order if mode is positive and in descending order if mode is
!      negative.  If mode is 0, no sorting is done.

!    index is an integer output array of length n set so that the sequence

!            num(index(i)), i = 1,2,...,n

!      is sorted according to the setting of mode.  If mode is 0,
!      index is not referenced.

!    last is an integer output array of length nmax + 1.  The index of num for
!      the last occurrence of l is last(l) for any l = 0,1,...,nmax unless
!      last(l) = 0.  In this case l does not appear in num.

!    next is an integer output array of length n.  If num(k) = l, then the
!      index of num for the previous occurrence of l is next(k) for any
!      l = 0,1,...,nmax unless next(k) = 0.
!      In this case there is no previous occurrence of l in num.

!  Argonne National Laboratory. MINPACK Project. July 1983.
!  Thomas F. Coleman, Burton S. Garbow, Jorge J. More'

!  **********

INTEGER :: i, j, jinc, jl, ju, k, l

!     Determine the arrays next and last.

last(0:nmax) = 0
DO  k = 1, n
  l = num(k)
  next(k) = last(l)
  last(l) = k
END DO
IF (mode == 0) RETURN

!     Store the pointers to the sorted array in index.

i = 1
IF (mode > 0) THEN
  jl = 0
  ju = nmax
  jinc = 1
ELSE
  jl = nmax
  ju = 0
  jinc = -1
END IF
DO  j = jl, ju, jinc
  k = last(j)
  DO