dsm Subroutine

public subroutine dsm(m, n, npairs, indrow, indcol, ngrp, maxgrp, mingrp, info, ipntr, jpntr)

The purpose of dsm is to determine an optimal or near- optimal consistent partition of the columns of a sparse m by n matrix a.

the sparsity pattern of the matrix a is specified by the arrays indrow and indcol. on input the indices for the non-zero elements of a are

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

the (indrow,indcol) pairs may be specified in any order. duplicate input pairs are permitted, but the subroutine eliminates them.

the subroutine partitions the columns of a into groups such that columns in the same group do not have a non-zero in the same row position. a partition of the columns of a with this property is consistent with the direct determination of a.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: m

number of rows of a (>0)
integer, intent(in) :: n

number of columns of a (>0)
integer, intent(in) :: npairs

number of (indrow,indcol) pairs used to describe the sparsity pattern of a (>0)
integer, intent(inout), dimension(npairs) :: indrow

an integer array of length npairs. 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 are in non-decreasing order. the column indices can be recovered from the array jpntr.
integer, intent(inout), dimension(npairs) :: indcol

an integer array of length npairs. on input indcol must contain the column indices of the non-zero elements 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.
integer, intent(out), dimension(n) :: ngrp

specifies the partition of the columns of a. column jcol belongs to group ngrp(jcol).
integer, intent(out) :: maxgrp

the number of groups in the partition of the columns of a.
integer, intent(out) :: mingrp

a lower bound for the number of groups in any consistent partition of the columns of a.
integer, intent(out) :: info

for normal termination info = 1. if m, n, or npairs is not positive, then info = 0. if the k-th element of indrow is not an integer between 1 and m or the k-th element of indcol is not an integer between 1 and n, then info = -k.
integer, intent(out), dimension(m+1) :: ipntr

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(m+1)-1 is then the number of non-zero elements of the matrix a.
integer, intent(out), dimension(n+1) :: jpntr

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 of the matrix a.

Calls

proc~~dsm~~CallsGraph proc~dsm dsm_module::dsm proc~degr dsm_module::degr proc~dsm->proc~degr proc~ido dsm_module::ido proc~dsm->proc~ido proc~numsrt dsm_module::numsrt proc~dsm->proc~numsrt proc~seq dsm_module::seq proc~dsm->proc~seq proc~setr dsm_module::setr proc~dsm->proc~setr proc~slo dsm_module::slo proc~dsm->proc~slo proc~srtdat dsm_module::srtdat proc~dsm->proc~srtdat proc~ido->proc~numsrt

Called by

proc~~dsm~~CalledByGraph proc~dsm dsm_module::dsm proc~dsm_wrapper numerical_differentiation_module::sparsity_pattern%dsm_wrapper proc~dsm_wrapper->proc~dsm proc~compute_sparsity_random numerical_differentiation_module::compute_sparsity_random proc~compute_sparsity_random->proc~dsm_wrapper proc~compute_sparsity_random_2 numerical_differentiation_module::compute_sparsity_random_2 proc~compute_sparsity_random_2->proc~dsm_wrapper proc~set_sparsity_pattern numerical_differentiation_module::numdiff_type%set_sparsity_pattern proc~set_sparsity_pattern->proc~dsm_wrapper

Source Code

    subroutine dsm(m,n,Npairs,Indrow,Indcol,Ngrp,Maxgrp,Mingrp,Info,Ipntr,Jpntr)

    implicit none

      integer,intent(in)  :: m       !! number of rows of `a` (>0)
      integer,intent(in)  :: n       !! number of columns of `a` (>0)
      integer,intent(in)  :: npairs  !! number of (`indrow`,`indcol`) pairs used
                                     !! to describe the sparsity pattern of `a` (>0)
      integer,intent(out) :: maxgrp  !! the number of groups in the partition
                                     !! of the columns of `a`.
      integer,intent(out) :: mingrp  !! a lower bound for the number of groups
                                     !! in any consistent partition of the
                                     !! columns of `a`.
      integer,intent(out) :: info    !! for normal termination `info = 1`.
                                     !! if `m`, `n`, or `npairs` is not positive,
                                     !! then `info = 0`. if the k-th element of
                                     !! `indrow` is not an integer between
                                     !! 1 and m or the k-th element of `indcol`
                                     !! is not an integer between 1 and n,
                                     !! then `info = -k`.
      integer,dimension(npairs),intent(inout) :: indrow  !! an integer array of length `npairs`. 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 are in non-decreasing order. the column
                                                         !! indices can be recovered from the array `jpntr`.
      integer,dimension(npairs),intent(inout) :: indcol  !! an integer array of length `npairs`. on input `indcol`
                                                         !! must contain the column indices of the non-zero elements 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`.
      integer,dimension(n),intent(out) :: ngrp   !! specifies the partition of the columns of `a`.
                                                 !! column `jcol` belongs to group `ngrp(jcol)`.
      integer,dimension(m+1),intent(out) :: ipntr   !! 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(m+1)-1` is then the number of non-zero
                                                    !! elements of the matrix `a`.
      integer,dimension(n+1),intent(out) :: jpntr   !! 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 of the matrix `a`.

    integer,dimension(max(m,6*n)) :: iwa     !! an integer work array
    integer :: i , ir , j , jp , k , maxclq , nnz , numgrp

    !  check the input data.

    Info = 0
    if ( m<1 .or. n<1 .or. Npairs<1 ) return
    do k = 1 , Npairs
        Info = -k
        if ( Indrow(k)<1 .or. Indrow(k)>m .or. Indcol(k)<1 .or. Indcol(k)>n ) return
    enddo
    Info = 1

    !  sort the data structure by columns.

    call srtdat(n,Npairs,Indrow,Indcol,Jpntr,Iwa)

    !  compress the data and determine the number of
    !  non-zero elements of a.

    do i = 1 , m
        Iwa(i) = 0
    enddo
    nnz = 1
    do j = 1 , n
        k = nnz
        do jp = Jpntr(j) , Jpntr(j+1) - 1
            ir = Indrow(jp)
            if ( Iwa(ir)/=j ) then
                Indrow(nnz) = ir
                nnz = nnz + 1
                Iwa(ir) = j
            endif
        enddo
        Jpntr(j) = k
    enddo
    Jpntr(n+1) = nnz

    !  extend the data structure to rows.

    call setr(m,n,Indrow,Jpntr,Indcol,Ipntr,Iwa)

    !  determine a lower bound for the number of groups.

    Mingrp = 0
    do i = 1 , m
        Mingrp = max(Mingrp,Ipntr(i+1)-Ipntr(i))
    enddo

    !  determine the degree sequence for the intersection
    !  graph of the columns of a.

    call degr(n,Indrow,Jpntr,Indcol,Ipntr,Iwa(5*n+1),Iwa(n+1))

    !  color the intersection graph of the columns of a
    !  with the smallest-last (sl) ordering.

    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)

    !  exit if the smallest-last ordering is optimal.

    if ( Maxgrp==Mingrp ) return

    !  color the intersection graph of the columns of a
    !  with the incidence-degree (id) ordering.

    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)

    !  retain the better of the two orderings so far.

    if ( numgrp<Maxgrp ) then
        Maxgrp = numgrp
        do j = 1 , n
            Ngrp(j) = Iwa(j)
        enddo

        !  exit if the incidence-degree ordering is optimal.

        if ( Maxgrp==Mingrp ) return
    endif

    !  color the intersection graph of the columns of a
    !  with the largest-first (lf) ordering.

    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))

    !  retain the best of the three orderings and exit.

    if ( numgrp<Maxgrp ) then
        Maxgrp = numgrp
        do j = 1 , n
            Ngrp(j) = Iwa(j)
        enddo
    endif

    end subroutine dsm