coloring.f90

牛顿优化算法源fortran代码

    IF (k == 0) EXIT
    INDEX(i) = k
    i = i + 1
    k = next(k)
  END DO
END DO
RETURN

!     Last card of subroutine numsrt.

END SUBROUTINE numsrt



SUBROUTINE sdpt(n, nghbrp, npntrp, nghbrs, npntrs, ngrp, maxgrp)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:49:17

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)  :: ngrp(:)
INTEGER, INTENT(OUT)  :: maxgrp

!  **********

!  subroutine sdpt

!  Given a loopless graph G = (V,E), this subroutine determines
!  a symmetric coloring of G by the Powell-Toint direct method.

!  The Powell-Toint method assigns the k-th color by examining the un-colored
!  vertices U(k) in order of non-increasing degree and assigning color k to
!  vertex v if there are no paths of length 1 or 2 (in the graph induced
!  by U(k)) between v and some k-colored vertex.

!  The subroutine statement is

!    subroutine sdpt(n, nghbrp, npntrp, nghbrs, npntrs, ngrp, maxgrp,
!                    iwa1, iwa2)

!  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.

!    ngrp is an integer output array of length n which specifies the
!      symmetric coloring of G. Vertex j is colored with color ngrp(j).

!    maxgrp is an integer output variable which specifies the
!      number of colors in the symmetric coloring of G.

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

!  Subprograms called

!    FORTRAN-supplied ... max

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

!  **********

INTEGER :: j, jp, k, kp, l, maxdeg, numdeg, numv, iwa1(0:n-1), iwa2(n)

!     Initialization block. Numv is the current number of un-colored
!     vertices, maxdeg is the maximum induced degree of these
!     vertices, and maxgrp is the current group number (color).

numv = n
maxdeg = 0
DO  j = 1, n
  ngrp(j) = (npntrp(j) - npntrp(j+1) + 1) + (npntrs(j) - npntrs(j+1) + 1)
  maxdeg = MAX(maxdeg,-ngrp(j))
  iwa2(j) = -j
END DO
maxgrp = 0

!     Beginning of iteration loop.

!        Sort the list of un-colored vertices so that their
!        induced degrees are in non-decreasing order.

20 iwa1(0:maxdeg) = 0
DO  l = 1, numv
  numdeg = -ngrp(-iwa2(l))
  iwa1(numdeg) = iwa1(numdeg) + 1
END DO
k = 1
DO  numdeg = maxdeg, 0, -1
  l = iwa1(numdeg)
  iwa1(numdeg) = k
  k = k + l
END DO
k = 1

60 j = iwa2(k)
IF (j > 0) THEN
  k = iwa1(-ngrp(j))
ELSE
  numdeg = -ngrp(-j)
  l = iwa1(numdeg)
  iwa2(k) = iwa2(l)
  iwa2(l) = -j
  iwa1(numdeg) = iwa1(numdeg) + 1
END IF
IF (k <= numv) GO TO 60
maxgrp = maxgrp + 1

!        Determine the vertices in group maxgrp.

loop150:   &
DO  l = 1, numv
  j = iwa2(l)
  
!          Examine each vertex k preceding vertex j and all the neighbors of
!          vertex k to determine if vertex j can be considered for group maxgrp.
  
  DO  jp = npntrp(j), npntrp(j+1)-1
    k = nghbrp(jp)
    IF (ngrp(k) == maxgrp) CYCLE loop150
    IF (ngrp(k) <= 0) THEN
      DO  kp = npntrp(k), npntrp(k+1)-1
        IF (ngrp(nghbrp(kp)) == maxgrp) CYCLE loop150
      END DO
      DO  kp = npntrs(k), npntrs(k+1)-1
        IF (ngrp(nghbrs(kp)) == maxgrp) CYCLE loop150
      END DO
    END IF
  END DO
  
!          Examine each vertex k succeeding vertex j and all the neighbors of
!          vertex k to determine if vertex j can be added to group maxgrp.
  
  DO  jp = npntrs(j), npntrs(j+1)-1
    k = nghbrs(jp)
    IF (ngrp(k) == maxgrp) CYCLE loop150
    IF (ngrp(k) <= 0) THEN
      DO  kp = npntrp(k), npntrp(k+1)-1
        IF (ngrp(nghbrp(kp)) == maxgrp) CYCLE loop150
      END DO
      DO  kp = npntrs(k), npntrs(k+1)-1
        IF (ngrp(nghbrs(kp)) == maxgrp) CYCLE loop150
      END DO
    END IF
  END DO
  
!           Add vertex j to group maxgrp and remove vertex j
!           from the list of un-colored vertices.
  
  ngrp(j) = maxgrp
  iwa2(l) = 0
  
!           Update the degrees of the neighbors of vertex j.
  
  DO  jp = npntrp(j), npntrp(j+1)-1
    k = nghbrp(jp)
    IF (ngrp(k) < 0) ngrp(k) = ngrp(k) + 1
  END DO
  DO  jp = npntrs(j), npntrs(j+1)-1
    k = nghbrs(jp)
    IF (ngrp(k) < 0) ngrp(k) = ngrp(k) + 1
  END DO

END DO loop150

!        Compress the updated list of un-colored vertices.
!        Reset numv and recompute maxdeg.

k = 0
maxdeg = 0
DO  l = 1, numv
  IF (iwa2(l) /= 0) THEN
    k = k + 1
    iwa2(k) = -iwa2(l)
    maxdeg = MAX(maxdeg,-ngrp(iwa2(l)))
  END IF
END DO
numv = k

!        End of iteration loop.

IF (numv > 0) GO TO 20
RETURN

!     Last card of subroutine sdpt.

END SUBROUTINE sdpt



SUBROUTINE seq(n, indrow, jpntr, indcol, ipntr, list, ngrp, maxgrp)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:49:22

INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(IN)   :: indrow(:)
INTEGER, INTENT(IN)   :: jpntr(n+1)
INTEGER, INTENT(IN)   :: indcol(:)
INTEGER, INTENT(IN)   :: ipntr(:)
INTEGER, INTENT(IN)   :: list(:)
INTEGER, INTENT(OUT)  :: ngrp(:)
INTEGER, INTENT(OUT)  :: maxgrp

!  **********

!  subroutine seq

!  Given the sparsity pattern of an m by n matrix A, this
!  subroutine determines a consistent partition of the
!  columns of A by a sequential algorithm.

!  A consistent partition is defined in terms of the loopless
!  graph G with vertices a(j), j = 1,2,...,n where a(j) is the
!  j-th column of A and with edge (a(i),a(j)) if and only if
!  columns i and j have a non-zero in the same row position.

!  A partition of the columns of A into groups is consistent
!  if the columns in any group are not adjacent in the graph G.
!  In graph-theory terminology, a consistent partition of the
!  columns of A corresponds to a coloring of the graph G.

!  The subroutine examines the columns in the order specified
!  by the array list, and assigns the current column to the
!  group with the smallest possible number.

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

!  The subroutine statement is

!    subroutine seq(n, indrow, jpntr, indcol, ipntr, list, ngrp, maxgrp, iwa)

!  where

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

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

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

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

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

!    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.

!    list is an integer input array of length n which specifies
!      the order to be used by the sequential algorithm.
!      The j-th column in this order is list(j).

!    ngrp is an integer output array of length n which specifies the
!      partition of the columns of A.  Column jcol belongs to group ngrp(jcol).

!    maxgrp is an integer output variable which specifies the
!      number of groups in the partition of the columns of A.

!    iwa is an integer work array of length n.

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

!  **********

INTEGER :: ic, ip, ir, j, jcol, jp, iwa(n)

!     Initialization block.

maxgrp = 0
DO  jp = 1, n
  ngrp(jp) = n
  iwa(jp) = 0
END DO

!     Beginning of iteration loop.

DO  j = 1, n
  jcol = list(j)
  
!        Find all columns adjacent to column jcol.
  
!        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 iwa marks the group numbers of the
!              columns which are adjacent to column jcol.
      
      iwa(ngrp(ic)) = j
    END DO
  END DO
  
!        Assign the smallest un-marked group number to jcol.
  
  DO  jp = 1, maxgrp
    IF (iwa(jp) /= j) GO TO 50
  END DO
  maxgrp = maxgrp + 1
  50 ngrp(jcol) = jp
END DO

!        End of iteration loop.

RETURN

!     Last card of subroutine seq.

END SUBROUTINE seq



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

INTEGER, INTENT(IN)   :: m
INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(IN)   :: indrow(:)
INTEGER, INTENT(IN)   :: jpntr(:)
INTEGER, INTENT(OUT)  :: indcol(:)
INTEGER, INTENT(OUT)  :: ipntr(:)

!  **********

!  subroutine setr

!  Given a column-oriented definition of the sparsity pattern
!  of an m by n matrix A, this subroutine determines a
!  row-oriented definition of the sparsity pattern of A.

!  On input the column-oriented definition is specified by
!  the arrays indrow and jpntr. On output the row-oriented
!  definition is specified by the arrays indcol and ipntr.

!  The subroutine statement is

!    subroutine setr(m, n, indrow, jpntr, indcol, ipntr, iwa)

!  where

!    m is a positive integer input variable set to the number of rows of A.

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

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

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

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

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

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

!    ipntr is an integer output 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(1) is set to 1 and that ipntr(m+1)-1 is
!      then the number of non-zero elements of the matrix A.

!    iwa is an integer work array of length m.

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

!  **********

INTEGER :: ir, jcol, jp, iwa(n)

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

iwa(1:m) = 0
DO  jp = 1, jpntr(n+1)-1
  iwa(indrow(jp)) = iwa(indrow(jp)) + 1
END DO

!     Set pointers to the start of the rows in indcol.

ipntr(1) = 1
DO  ir = 1, m
  ipntr(ir+1) = ipntr(ir) + iwa(ir)
  iwa(ir) = ipntr(ir)
END DO

!     Fill indcol.

DO  jcol = 1, n
  DO  jp = jpntr(jcol), jpntr(jcol+1)-1
    ir = indrow(jp)
    indcol(iwa(ir)) = jcol
    iwa(ir) = iwa(ir) + 1
  END DO
END DO
RETURN

!     Last card of subroutine setr.

END SUBROUTINE setr



SUBROUTINE slo(n, indrow, jpntr, indcol, ipntr, ndeg, list, maxclq)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:49:41

INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(IN)   :: indrow(:)
INTEGER, INTENT(IN)   :: jpntr(:)
INTEGER, INTENT(IN)   :: indcol(:)
INTEGER, INTENT(IN)   :: ipntr(:)
INTEGER, INTENT(IN)   :: ndeg(:)
INTEGER, INTENT(OUT)  :: list(:)
INTEGER, INTENT(OUT)  :: maxclq

!  **********

!  subroutine slo

!  Given the sparsity pattern of an m by n matrix A, this subroutine
!  determines the smallest-last ordering of the columns of A.

!  The smallest-last ordering is defined for the loopless
!  graph G with vertices a(j), j = 1,2,...,n where a(j) is the
!  j-th column of A and with edge (a(i),a(j)) if and only if
!  columns i and j have a non-zero in the same row position.

!  The smallest-last ordering is determined recursively by
!  letting list(k), k = n,...,1 be a column with least degree
!  in the subgraph spanned by the un-ordered columns.

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

!  The subroutine statement is

!    subroutine slo(n, indrow, jpntr, indcol, ipntr, ndeg, list,
!                   maxclq, iwa1, iwa2, iwa3, iwa4)

!  where

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

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

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

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

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