coloring.f90

牛顿优化算法源fortran代码

MODULE coloring
! A module from MINPACK-2, used by DTRON

IMPLICIT NONE

INTEGER, PARAMETER, PRIVATE :: dp = SELECTED_REAL_KIND(14, 60)


CONTAINS


SUBROUTINE dssm(n, npairs, indrow, indcol, method, listp, ngrp,  &
                maxgrp, mingrp, info, ipntr, jpntr)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:40:40
! Latest revision - 21 November 1999

INTEGER, INTENT(IN)      :: n
INTEGER, INTENT(IN)      :: npairs
INTEGER, INTENT(IN OUT)  :: indrow(:)    ! indrow(npairs)
INTEGER, INTENT(IN OUT)  :: indcol(:)    ! indcol(npairs)
INTEGER, INTENT(IN)      :: method
INTEGER, INTENT(OUT)     :: listp(:)
INTEGER, INTENT(OUT)     :: ngrp(:)
INTEGER, INTENT(OUT)     :: maxgrp
INTEGER, INTENT(OUT)     :: mingrp
INTEGER, INTENT(OUT)     :: info
INTEGER, INTENT(IN OUT)  :: ipntr(:)    ! ipntr(n+1)
INTEGER, INTENT(IN OUT)  :: jpntr(:)    ! jpntr(n+1)

!  **********

!  subroutine dssm

!  Given the sparsity pattern of a symmetric matrix A of order n,
!  this subroutine determines a symmetric permutation of A and a
!  partition of the columns of A consistent with the determination
!  of A by a lower triangular substitution method.

!  The sparsity pattern of the matrix A is specified by the arrays indrow and
!  indcol.  On input the indices for the non-zero elements in the lower
!  triangular part of A are

!        (indrow(k),indcol(k)), k = 1,2,...,npairs.

!  The (indrow(k),indcol(k)) pairs may be specified in any order.
!  Duplicate input pairs are permitted, but the subroutine eliminates them.
!  The subroutine requires that all the diagonal elements be part of the
!  sparsity pattern and replaces any pair (indrow(k),indcol(k)) where indrow(k)
!  is less than indcol(k) by the pair (indcol(k),indrow(k)).

!  The direct method (method = 1) first determines a partition of the columns
!  of A such that two columns in a group have a non-zero element in row k only
!  if column k is in an earlier group.  Using this partition, the subroutine
!  then computes a symmetric permutation of A consistent with the determination
!  of A by a lower triangular substitution method.

!  The indirect method first computes a symmetric permutation of A which
!  minimizes the maximum number of non-zero elements in any row of L, where L
!  is the lower triangular part of the permuted matrix.  The subroutine then
!  partitions the columns of L into groups such that columns of L in a group
!  do not have a non-zero in the same row position.

!  The subroutine statement is

!    subroutine dssm(n, npairs, indrow, indcol, method, listp, ngrp,
!                    maxgrp, mingrp, info, ipntr, jpntr)

!  where

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

!    npairs is a positive integer input variable set to the number of
!      (indrow,indcol) pairs used to describe the sparsity pattern of A.

!    indrow is an integer array of length npairs.  On input indrow must contain
!      the row indices of the non-zero elements in the lower triangular part
!      of A.  On output indrow is permuted so that the corresponding column
!      indices are in non-decreasing order.  The column indices can be
!      recovered from the array jpntr.

!    indcol is an integer array of length npairs.  On input indcol must contain
!      the column indices of the non-zero elements in the lower triangular part
!      of A.  On output indcol is permuted so that the corresponding row
!      indices are in non-decreasing order.  The row indices can be recovered
!      from the array ipntr.

!    method is an integer input variable.  If method = 1, the direct method
!      is used to determine the partition and symmetric permutation.
!      Otherwise, the indirect method is used to determine the symmetric
!      permutation and partition.

!    listp is an integer output array of length n which specifies
!      the symmetric permutation of the matrix A.  Element (i,j)
!      of A is the (listp(i),listp(j)) element of the permuted matrix.

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

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

!    mingrp is an integer output variable which specifies a lower bound for
!      the number of groups in any partition of the columns of A consistent
!      with the determination of A by a lower triangular substitution method.

!    info is an integer output variable set as follows. For
!      normal termination info = 1. If n or npairs is not
!      positive or liwa is less than 6*n, then info = 0. If the
!      k-th element of indrow or the k-th element of indcol is
!      not an integer between 1 and n, or if the k-th diagonal
!      element is not in the sparsity pattern, then info = -k.

!    ipntr is an integer output array of length n + 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(n+1)-1 is then the number of non-zero
!      elements in the lower triangular part of the matrix A.

!    jpntr is an integer output 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 in the lower triangular part of the matrix A.

!  Subprograms called

!    MINPACK-supplied ... degr,ido,idog,numsrt,sdpt,seq,setr,slo,slog,srtdat

!    FORTRAN-supplied ... max,min

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

!  **********

INTEGER :: i, ir, j, jp, k, maxid, maxvd, maxclq, nnz, numgrp
INTEGER :: iwa(3*n)

!     Check the input data.

info = 0
IF (n < 1 .OR. npairs < 1) RETURN
iwa(1:n) = 0
DO  k = 1, npairs
  info = -k
  IF (indrow(k) < 1 .OR. indrow(k) > n .OR.  &
      indcol(k) < 1 .OR. indcol(k) > n) RETURN
  IF (indrow(k) == indcol(k)) iwa(indrow(k)) = 1
END DO
DO  k = 1, n
  info = -k
  IF (iwa(k) /= 1) RETURN
END DO
info = 1

!     Generate the sparsity pattern for the lower triangular part of A.

DO  k = 1, npairs
  i = indrow(k)
  j = indcol(k)
  indrow(k) = MAX(i,j)
  indcol(k) = MIN(i,j)
END DO

!     Sort the data structure by columns.

CALL srtdat(n, npairs, indrow, indcol, jpntr)

!     Compress the data and determine the number of non-zero elements
!     in the lower triangular part of A.

iwa(1:n) = 0
nnz = 0
DO  j = 1, n
  k = nnz
  DO  jp = jpntr(j), jpntr(j+1)-1
    ir = indrow(jp)
    IF (iwa(ir) /= j) THEN
      nnz = nnz + 1
      indrow(nnz) = ir
      iwa(ir) = j
    END IF
  END DO
  jpntr(j) = k + 1
END DO
jpntr(n+1) = nnz + 1

!     Extend the data structure to rows.

CALL setr(n, n, indrow, jpntr, indcol, ipntr)

!     Determine the smallest-last ordering of the vertices of the adjacency
!     graph of A, and from it determine a lower bound for the number of groups.

CALL slog(n, indrow, jpntr, indcol, ipntr, iwa, maxclq, maxvd)
mingrp =  1 + maxvd

!     Use the selected method.

IF (method == 1) THEN
  
!        Direct method. Determine a partition of the columns
!        of A by the Powell-Toint method.
  
  CALL sdpt(n, indrow, jpntr, indcol, ipntr, ngrp, maxgrp)
  
!        Define a symmetric permutation of A according to the
!        ordering of the column group numbers in the partition.
  
  CALL numsrt(n, maxgrp, ngrp, 1, iwa, iwa(2*n+1:), iwa(n+1:))
  DO  i = 1, n
    listp(iwa(i)) = i
  END DO
ELSE
  
!        Indirect method. Determine the incidence degree ordering of the
!        vertices of the adjacency graph of A and, together with the
!        smallest-last ordering, define a symmetric permutation of A.
  
  CALL idog(n, indrow, jpntr, indcol, ipntr, listp, maxclq, maxid)
  IF (maxid > maxvd) listp(1:n) = iwa(1:n)
  
!        Generate the sparsity pattern for the lower
!        triangular part L of the permuted matrix.
  
  DO  j = 1, n
    DO  jp = jpntr(j), jpntr(j+1)-1
      i = indrow(jp)
      indrow(jp) = MAX(listp(i),listp(j))
      indcol(jp) = MIN(listp(i),listp(j))
    END DO
  END DO
  
!        Sort the data structure by columns.
  
  CALL srtdat(n, nnz, indrow, indcol, jpntr)
  
!        Extend the data structure to rows.
  
  CALL setr(n, n, indrow, jpntr, indcol, ipntr)
  
!        Determine the degree sequence for the intersection
!        graph of the columns of L.
  
  CALL degr(n, indrow, jpntr, indcol, ipntr, iwa(2*n+1:))
  
!        Color the intersection graph of the columns of L
!        with the smallest-last (SL) ordering.
  
  CALL slo(n, indrow, jpntr, indcol, ipntr, iwa(2*n+1:), iwa(n+1:), maxclq)
  CALL seq(n, indrow, jpntr, indcol, ipntr, iwa(n+1:), iwa, maxgrp)
  DO  j = 1, n
    ngrp(j) = iwa(listp(j))
  END DO
  
!        Exit if the smallest-last ordering is optimal.
  
  IF (maxgrp /= maxclq) THEN
  
!        Color the intersection graph of the columns of L
!        with the incidence degree (ID) ordering.

    CALL ido(n, n, indrow, jpntr, indcol, ipntr, iwa(2*n+1:), iwa(n+1:), maxclq)
    CALL seq(n, indrow, jpntr, indcol, ipntr, iwa(n+1:), iwa, numgrp)

!        Retain the better of the two orderings.
  
    IF (numgrp < maxgrp) THEN
      maxgrp = numgrp
      DO  j = 1, n
        ngrp(j) = iwa(listp(j))
      END DO
    END IF
  END IF
  
!        Generate the sparsity pattern for the lower
!        triangular part of the original matrix.
  
  DO  j = 1, n
    iwa(listp(j)) = j
  END DO
  DO  j = 1, n
    DO  jp = jpntr(j), jpntr(j+1)-1
      i = indrow(jp)
      indrow(jp) = MAX(iwa(i),iwa(j))
      indcol(jp) = MIN(iwa(i),iwa(j))
    END DO
  END DO
  
!        Sort the data structure by columns.
  
  CALL srtdat(n, nnz, indrow, indcol, jpntr)
  
!        Extend the data structure to rows.
  
  CALL setr(n, n, indrow, jpntr, indcol, ipntr)
END IF
RETURN

!     Last card of subroutine dssm.

END SUBROUTINE dssm



SUBROUTINE degr(n, indrow, jpntr, indcol, ipntr, ndeg)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-30  Time: 12:48:53

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

!  **********

!  subroutine degr

!  Given the sparsity pattern of an m by n matrix A, this subroutine determines
!  the degree sequence for the intersection graph of the columns of A.

!  In graph-theory terminology, the intersection graph of the columns of A is
!  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.

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

!  The subroutine statement is

!    subroutine degr(n, indrow, jpntr, indcol, ipntr, ndeg, 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.

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

!    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, jcol, jp, iwa(n)

!     Initialization block.

DO  jp = 1, n
  ndeg(jp) = 0
  iwa(jp) = 0
END DO

!     Compute the degree sequence by determining the contributions
!     to the degrees from the current(jcol) column and further
!     columns which have not yet been considered.

DO  jcol = 2, n
  iwa(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 iwa marks columns which have contributed to
!              the degree count of column jcol. Update the degree
!              counts of these columns as well as column jcol.
      
      IF (iwa(ic) < jcol) THEN
        iwa(ic) = jcol
        ndeg(ic) = ndeg(ic) + 1
        ndeg(jcol) = ndeg(jcol) + 1
      END IF
    END DO
  END DO
END DO
RETURN

!     Last card of subroutine degr.

END SUBROUTINE degr



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

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

!  **********

!  subroutine ido

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

!  The incidence-degree 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 incidence-degree ordering is determined recursively by letting
!  list(k), k = 1,...,n be a column with maximal incidence to the subgraph
!  spanned by the ordered columns.  Among all the columns of maximal incidence,
!  ido chooses a column of maximal degree.

!  The subroutine statement is

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

!  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 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 incidence-degree 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

!    MINPACK-supplied ... numsrt