coloring.f90

牛顿优化算法源fortran代码

!    indcol is an integer input array which contains the
!      column indices for the non-zeroes in the matrix A.

!    ipntr is an integer input array of length m + 1 which
!      specifies the locations of the column indices in indcol.
!      The column indices for row i are

!            indcol(k), k = ipntr(i),...,ipntr(i+1)-1.

!      Note that ipntr(m+1)-1 is then the number of non-zero
!      elements of the matrix A.

!    ndeg is an integer input array of length n which specifies
!      the degree sequence. The degree of the j-th column
!      of A is ndeg(j).

!    list is an integer output array of length n which specifies
!      the smallest-last ordering of the columns of A.  The j-th
!      column in this order is list(j).

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

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

!  Subprograms called

!    FORTRAN-supplied ... min

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

!  **********

INTEGER :: ic, ip, ir, jcol, jp, mindeg, numdeg, numord,  &
           iwa1(0:n-1), iwa2(n), iwa3(n), iwa4(n)

!     Initialization block.

mindeg = n
DO  jp = 1, n
  iwa1(jp-1) = 0
  iwa4(jp) = n
  list(jp) = ndeg(jp)
  mindeg = MIN(mindeg, ndeg(jp))
END DO

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

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

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

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

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

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

DO  jp = 1, n
  numdeg = ndeg(jp)
  iwa2(jp) = 0
  iwa3(jp) = iwa1(numdeg)
  IF (iwa1(numdeg) > 0) iwa2(iwa1(numdeg)) = jp
  iwa1(numdeg) = jp
END DO
maxclq = 0
numord = n

!     Beginning of iteration loop.

!        Choose a column jcol of minimal degree mindeg.

30 DO
  jcol = iwa1(mindeg)
  IF (jcol > 0) EXIT
  mindeg = mindeg + 1
END DO

list(jcol) = numord

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

IF (mindeg+1 == numord .AND. maxclq == 0) maxclq = numord

!        Termination test.

numord = numord - 1
IF (numord == 0) GO TO 80

!        Delete column jcol from the mindeg list.

iwa1(mindeg) = iwa3(jcol)
IF (iwa3(jcol) > 0) iwa2(iwa3(jcol)) = 0

!        Find all columns adjacent to column jcol.

iwa4(jcol) = 0

!        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 degree lists.
      
      numdeg = list(ic)
      list(ic) = list(ic) - 1
      mindeg = MIN(mindeg,list(ic))
      
!                 Delete column ic from the numdeg list.
      
      IF (iwa2(ic) == 0) THEN
        iwa1(numdeg) = iwa3(ic)
      ELSE
        iwa3(iwa2(ic)) = iwa3(ic)
      END IF
      IF (iwa3(ic) > 0) iwa2(iwa3(ic)) = iwa2(ic)
      
!                 Add column ic to the numdeg-1 list.
      
      iwa2(ic) = 0
      iwa3(ic) = iwa1(numdeg-1)
      IF (iwa1(numdeg-1) > 0) iwa2(iwa1(numdeg-1)) = ic
      iwa1(numdeg-1) = ic
    END IF
  END DO
END DO

!        End of iteration loop.

GO TO 30

!     Invert the array list.

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

!     Last card of subroutine slo.

END SUBROUTINE slo



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

INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(IN)   :: nghbrp(:)
INTEGER, INTENT(IN)   :: npntrp(:)
INTEGER, INTENT(IN)   :: nghbrs(:)
INTEGER, INTENT(IN)   :: npntrs(:)
INTEGER, INTENT(OUT)  :: listp(:)
INTEGER, INTENT(OUT)  :: maxclq
INTEGER, INTENT(OUT)  :: maxvd

!  **********

!  subroutine slog

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

!  The smallest-last ordering is determined recursively by
!  letting list(k), k = n,...,1 be a vertex with least degree
!  in the subgraph spanned by the un-ordered vertices.
!  This subroutine determines the inverse of the smallest-last
!  ordering, that is, an array listp such that listp(list(k)) = k
!  for k = 1,2,...,n.

!  The subroutine statement is

!    subroutine slog(n, nghbrp, npntrp, nghbrs, npntrs, listp,
!                    maxclq, maxvd, 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 smallest-last 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.

!    maxvd is an integer output variable set to the maximum
!      vertex degree found during the ordering.

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

!  Subprograms called

!    FORTRAN-supplied ... max,min

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

!  **********

INTEGER :: i, j, k, mindeg, numdeg, numord, iwa1(0:n-1), iwa2(n), iwa3(n)

!     Initialization block.

mindeg = n
DO  j = 1, n
  iwa1(j-1) = 0
  listp(j) = (npntrp(j) - npntrp(j+1) + 1) + (npntrs(j) - npntrs(j+1) + 1)
  mindeg = MIN(mindeg,-listp(j))
END DO

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

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

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

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

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

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

DO  j = 1, n
  numdeg = -listp(j)
  iwa2(j) = 0
  iwa3(j) = iwa1(numdeg)
  IF (iwa1(numdeg) > 0) iwa2(iwa1(numdeg)) = j
  iwa1(numdeg) = j
END DO
maxclq = 0
maxvd = 0
numord = n

!     Beginning of iteration loop.

!        Choose a vertex j of minimal degree mindeg.

30 DO
  j = iwa1(mindeg)
  IF (j > 0) EXIT
  mindeg = mindeg + 1
END DO

listp(j) = numord
maxvd = MAX(maxvd, mindeg)

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

IF (mindeg+1 == numord .AND. maxclq == 0) maxclq = numord

!        Termination test.

numord = numord - 1
IF (numord == 0) GO TO 80

!        Delete vertex j from the mindeg list.

iwa1(mindeg) = iwa3(j)
IF (iwa3(j) > 0) iwa2(iwa3(j)) = 0

!        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 degree lists.
  
  numdeg = -listp(i)
  IF (numdeg >= 0) THEN
    listp(i) = listp(i) + 1
    mindeg = MIN(mindeg,-listp(i))
    
!              Delete vertex i from the numdeg list.
    
    IF (iwa2(i) == 0) THEN
      iwa1(numdeg) = iwa3(i)
    ELSE
      iwa3(iwa2(i)) = iwa3(i)
    END IF
    IF (iwa3(i) > 0) iwa2(iwa3(i)) = iwa2(i)
    
!              Add vertex i to the numdeg-1 list.
    
    iwa2(i) = 0
    iwa3(i) = iwa1(numdeg-1)
    IF (iwa1(numdeg-1) > 0) iwa2(iwa1(numdeg-1)) = i
    iwa1(numdeg-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 degree lists.
  
  numdeg = -listp(i)
  IF (numdeg >= 0) THEN
    listp(i) = listp(i) + 1
    mindeg = MIN(mindeg,-listp(i))
    
!              Delete vertex i from the numdeg list.
    
    IF (iwa2(i) == 0) THEN
      iwa1(numdeg) = iwa3(i)
    ELSE
      iwa3(iwa2(i)) = iwa3(i)
    END IF
    IF (iwa3(i) > 0) iwa2(iwa3(i)) = iwa2(i)
    
!              Add vertex i to the numdeg-1 list.
    
    iwa2(i) = 0
    iwa3(i) = iwa1(numdeg-1)
    IF (iwa1(numdeg-1) > 0) iwa2(iwa1(numdeg-1)) = i
    iwa1(numdeg-1) = i
  END IF
END DO

!        End of iteration loop.

GO TO 30
80 RETURN

!     Last card of subroutine slog.

END SUBROUTINE slog



SUBROUTINE srtdat(n, nnz, indrow, indcol, jpntr)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:49:52

INTEGER, INTENT(IN)      :: n
INTEGER, INTENT(IN)      :: nnz
INTEGER, INTENT(IN OUT)  :: indrow(:)
INTEGER, INTENT(IN OUT)  :: indcol(:)
INTEGER, INTENT(OUT)     :: jpntr(:)    ! jpntr(n+1)

!  **********

!  subroutine srtdat

!  Given the non-zero elements of an m by n matrix A in
!  arbitrary order as specified by their row and column
!  indices, this subroutine permutes these elements so
!  that their column indices are in non-decreasing order.

!  On input it is assumed that the elements are specified in

!        indrow(k),indcol(k), k = 1,...,nnz.

!  On output the elements are permuted so that indcol is
!  in non-decreasing order. In addition, the array jpntr
!  is set so that the row indices for column j are

!        indrow(k), k = jpntr(j),...,jpntr(j+1)-1.

!  Note that the value of m is not needed by srtdat and is
!  therefore not present in the subroutine statement.

!  The subroutine statement is

!    subroutine srtdat(n,nnz,indrow,indcol,jpntr,iwa)

!  where

!    n is a positive integer input variable set to the number of columns of A.

!    nnz is a positive integer input variable set to the number
!      of non-zero elements of A.

!    indrow is an integer array of length nnz. On input indrow
!      must contain the row indices of the non-zero elements of A.
!      On output indrow is permuted so that the corresponding
!      column indices of indcol are in non-decreasing order.

!    indcol is an integer array of length nnz. On input indcol
!      must contain the column indices of the non-zero elements
!      of A. On output indcol is permuted so that these indices
!      are in non-decreasing order.

!    jpntr is an integer output array of length n + 1 which
!      specifies the locations of the row indices in the output
!      indrow. The row indices for column j are

!            indrow(k), k = jpntr(j),...,jpntr(j+1)-1.

!      Note that jpntr(1) is set to 1 and that jpntr(n+1)-1 is then nnz.

!    iwa is an integer work array of length n.

!  Subprograms called

!    FORTRAN-supplied ... max

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

!  **********

INTEGER :: i, j, k, l, iwa(n)

!     Store in array iwa the counts of non-zeroes in the columns.

iwa(1:n) = 0
DO  k = 1, nnz
  iwa(indcol(k)) = iwa(indcol(k)) + 1
END DO

!     Set pointers to the start of the columns in indrow.

jpntr(1) = 1
DO  j = 1, n
  jpntr(j+1) = jpntr(j) + iwa(j)
  iwa(j) = jpntr(j)
END DO
k = 1

!     Begin in-place sort.

40 j = indcol(k)
IF (k >= jpntr(j)) THEN
  
!           Current element is in position.  Now examine the next element
!           or the first un-sorted element in the j-th group.
  
  k = MAX(k+1,iwa(j))
ELSE
  
!           Current element is not in position.  Place element in position
!           and make the displaced element the current element.
  
  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 <= nnz) GO TO 40
RETURN

!     Last card of subroutine srtdat.

END SUBROUTINE srtdat

END MODULE coloring