Amesos Package Browser (Single Doxygen Collection) Development
amesos_btf_decl.h
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === BTF package ========================================================== */
00003 /* ========================================================================== */
00004 
00005 /* BTF_MAXTRANS:  find a column permutation Q to give A*Q a zero-free diagonal
00006  * BTF_STRONGCOMP:  find a symmetric permutation P to put P*A*P' into block
00007  *  upper triangular form.
00008  * BTF_ORDER: do both of the above (btf_maxtrans then btf_strongcomp).
00009  *
00010  * Copyright (c) 2004-2007.  Tim Davis, University of Florida,
00011  * with support from Sandia National Laboratories.  All Rights Reserved.
00012  */
00013 
00014 
00015 /* ========================================================================== */
00016 /* === BTF_MAXTRANS ========================================================= */
00017 /* ========================================================================== */
00018 
00019 /* BTF_MAXTRANS: finds a permutation of the columns of a matrix so that it has a
00020  * zero-free diagonal.  The input is an m-by-n sparse matrix in compressed
00021  * column form.  The array Ap of size n+1 gives the starting and ending
00022  * positions of the columns in the array Ai.  Ap[0] must be zero. The array Ai
00023  * contains the row indices of the nonzeros of the matrix A, and is of size
00024  * Ap[n].  The row indices of column j are located in Ai[Ap[j] ... Ap[j+1]-1].
00025  * Row indices must be in the range 0 to m-1.  Duplicate entries may be present
00026  * in any given column.  The input matrix  is not checked for validity (row
00027  * indices out of the range 0 to m-1 will lead to an undeterminate result -
00028  * possibly a core dump, for example).  Row indices in any given column need
00029  * not be in sorted order.  However, if they are sorted and the matrix already
00030  * has a zero-free diagonal, then the identity permutation is returned.
00031  *
00032  * The output of btf_maxtrans is an array Match of size n.  If row i is matched
00033  * with column j, then A(i,j) is nonzero, and then Match[i] = j.  If the matrix
00034  * is structurally nonsingular, all entries in the Match array are unique, and
00035  * Match can be viewed as a column permutation if A is square.  That is, column
00036  * k of the original matrix becomes column Match[k] of the permuted matrix.  In
00037  * MATLAB, this can be expressed as (for non-structurally singular matrices):
00038  *
00039  *  Match = maxtrans (A) ;
00040  *  B = A (:, Match) ;
00041  *
00042  * except of course here the A matrix and Match vector are all 0-based (rows
00043  * and columns in the range 0 to n-1), not 1-based (rows/cols in range 1 to n).
00044  * The MATLAB dmperm routine returns a row permutation.  See the maxtrans
00045  * mexFunction for more details.
00046  *
00047  * If row i is not matched to any column, then Match[i] is == -1.  The
00048  * btf_maxtrans routine returns the number of nonzeros on diagonal of the
00049  * permuted matrix.
00050  *
00051  * In the MATLAB mexFunction interface to btf_maxtrans, 1 is added to the Match
00052  * array to obtain a 1-based permutation.  Thus, in MATLAB where A is m-by-n:
00053  *
00054  *  q = maxtrans (A) ;  % has entries in the range 0:n
00055  *  q     % a column permutation (only if sprank(A)==n)
00056  *  B = A (:, q) ;    % permuted matrix (only if sprank(A)==n)
00057  *  sum (q > 0) ;   % same as "sprank (A)"
00058  *
00059  * This behaviour differs from p = dmperm (A) in MATLAB, which returns the
00060  * matching as p(j)=i if row i and column j are matched, and p(j)=0 if column j
00061  * is unmatched.
00062  *
00063  *  p = dmperm (A) ;  % has entries in the range 0:m
00064  *  p                       % a row permutation (only if sprank(A)==m)
00065  *  B = A (p, :) ;    % permuted matrix (only if sprank(A)==m)
00066  *  sum (p > 0) ;   % definition of sprank (A)
00067  *
00068  * This algorithm is based on the paper "On Algorithms for obtaining a maximum
00069  * transversal" by Iain Duff, ACM Trans. Mathematical Software, vol 7, no. 1,
00070  * pp. 315-330, and "Algorithm 575: Permutations for a zero-free diagonal",
00071  * same issue, pp. 387-390.  Algorithm 575 is MC21A in the Harwell Subroutine
00072  * Library.  This code is not merely a translation of the Fortran code into C.
00073  * It is a completely new implementation of the basic underlying method (depth
00074  * first search over a subgraph with nodes corresponding to columns matched so
00075  * far, and cheap matching).  This code was written with minimal observation of
00076  * the MC21A/B code itself.  See comments below for a comparison between the
00077  * maxtrans and MC21A/B codes.
00078  *
00079  * This routine operates on a column-form matrix and produces a column
00080  * permutation.  MC21A uses a row-form matrix and produces a row permutation.
00081  * The difference is merely one of convention in the comments and interpretation
00082  * of the inputs and outputs.  If you want a row permutation, simply pass a
00083  * compressed-row sparse matrix to this routine and you will get a row
00084  * permutation (just like MC21A).  Similarly, you can pass a column-oriented
00085  * matrix to MC21A and it will happily return a column permutation.
00086  */
00087 
00088 #ifndef AMESOS_BTF_DECL_H
00089 #define AMESOS_BTF_DECL_H
00090 
00091 /* make it easy for C++ programs to include BTF */
00092 #ifdef __cplusplus
00093 extern "C" {
00094 #endif
00095 
00096 #include "amesos_UFconfig.h"
00097 
00098 int amesos_btf_maxtrans    /* returns # of columns matched */
00099 (
00100     /* --- input, not modified: --- */
00101     int nrow,     /* A is nrow-by-ncol in compressed column form */
00102     int ncol,
00103     int Ap [ ],     /* size ncol+1 */
00104     int Ai [ ],     /* size nz = Ap [ncol] */
00105     double maxwork, /* maximum amount of work to do is maxwork*nnz(A); no limit
00106          * if <= 0 */
00107 
00108     /* --- output, not defined on input --- */
00109     double *work,   /* work = -1 if maxwork > 0 and the total work performed
00110          * reached the maximum of maxwork*nnz(A).
00111          * Otherwise, work = the total work performed. */
00112 
00113     int Match [ ],  /* size nrow.  Match [i] = j if column j matched to row i
00114          * (see above for the singular-matrix case) */
00115 
00116     /* --- workspace, not defined on input or output --- */
00117     int Work [ ]    /* size 5*ncol */
00118 ) ;
00119 
00120 /* long integer version (all "int" parameters become "UF_long") */
00121 UF_long amesos_btf_l_maxtrans (UF_long, UF_long, UF_long *, UF_long *, double,
00122     double *, UF_long *, UF_long *) ;
00123 
00124 
00125 /* ========================================================================== */
00126 /* === BTF_STRONGCOMP ======================================================= */
00127 /* ========================================================================== */
00128 
00129 /* BTF_STRONGCOMP finds the strongly connected components of a graph, returning
00130  * a symmetric permutation.  The matrix A must be square, and is provided on
00131  * input in compressed-column form (see BTF_MAXTRANS, above).  The diagonal of
00132  * the input matrix A (or A*Q if Q is provided on input) is ignored.
00133  *
00134  * If Q is not NULL on input, then the strongly connected components of A*Q are
00135  * found.  Q may be flagged on input, where Q[k] < 0 denotes a flagged column k.
00136  * The permutation is j = BTF_UNFLIP (Q [k]).  On output, Q is modified (the
00137  * flags are preserved) so that P*A*Q is in block upper triangular form.
00138  *
00139  * If Q is NULL, then the permutation P is returned so that P*A*P' is in upper
00140  * block triangular form.
00141  *
00142  * The vector R gives the block boundaries, where block b is in rows/columns
00143  * R[b] to R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the
00144  * number of strongly connected components found.
00145  */
00146 
00147 int amesos_btf_strongcomp  /* return # of strongly connected components */
00148 (
00149     /* input, not modified: */
00150     int n,      /* A is n-by-n in compressed column form */
00151     int Ap [ ],     /* size n+1 */
00152     int Ai [ ],     /* size nz = Ap [n] */
00153 
00154     /* optional input, modified (if present) on output: */
00155     int Q [ ],      /* size n, input column permutation */
00156 
00157     /* output, not defined on input */
00158     int P [ ],      /* size n.  P [k] = j if row and column j are kth row/col
00159          * in permuted matrix. */
00160 
00161     int R [ ],      /* size n+1.  block b is in rows/cols R[b] ... R[b+1]-1 */
00162 
00163     /* workspace, not defined on input or output */
00164     int Work [ ]    /* size 4n */
00165 ) ;
00166 
00167 UF_long amesos_btf_l_strongcomp (UF_long, UF_long *, UF_long *, UF_long *, UF_long *,
00168     UF_long *, UF_long *) ;
00169 
00170 
00171 /* ========================================================================== */
00172 /* === BTF_ORDER ============================================================ */
00173 /* ========================================================================== */
00174 
00175 /* BTF_ORDER permutes a square matrix into upper block triangular form.  It
00176  * does this by first finding a maximum matching (or perhaps a limited matching
00177  * if the work is limited), via the btf_maxtrans function.  If a complete
00178  * matching is not found, BTF_ORDER completes the permutation, but flags the
00179  * columns of P*A*Q to denote which columns are not matched.  If the matrix is
00180  * structurally rank deficient, some of the entries on the diagonal of the
00181  * permuted matrix will be zero.  BTF_ORDER then calls btf_strongcomp to find
00182  * the strongly-connected components.
00183  *
00184  * On output, P and Q are the row and column permutations, where i = P[k] if
00185  * row i of A is the kth row of P*A*Q, and j = BTF_UNFLIP(Q[k]) if column j of
00186  * A is the kth column of P*A*Q.  If Q[k] < 0, then the (k,k)th entry in P*A*Q
00187  * is structurally zero.
00188  *
00189  * The vector R gives the block boundaries, where block b is in rows/columns
00190  * R[b] to R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the
00191  * number of strongly connected components found.
00192  */
00193 
00194 int amesos_btf_order      /* returns number of blocks found */
00195 (
00196     /* --- input, not modified: --- */
00197     int n,      /* A is n-by-n in compressed column form */
00198     int Ap [ ],     /* size n+1 */
00199     int Ai [ ],     /* size nz = Ap [n] */
00200     double maxwork, /* do at most maxwork*nnz(A) work in the maximum
00201          * transversal; no limit if <= 0 */
00202 
00203     /* --- output, not defined on input --- */
00204     double *work,   /* return value from amesos_btf_maxtrans */
00205     int P [ ],      /* size n, row permutation */
00206     int Q [ ],      /* size n, column permutation */
00207     int R [ ],      /* size n+1.  block b is in rows/cols R[b] ... R[b+1]-1 */
00208     int *nmatch,    /* # nonzeros on diagonal of P*A*Q */
00209 
00210     /* --- workspace, not defined on input or output --- */
00211     int Work [ ]    /* size 5n */
00212 ) ;
00213 
00214 UF_long amesos_btf_l_order (UF_long, UF_long *, UF_long *, double , double *,
00215     UF_long *, UF_long *, UF_long *, UF_long *, UF_long *) ;
00216 
00217 
00218 /* ========================================================================== */
00219 /* === BTF marking of singular columns ====================================== */
00220 /* ========================================================================== */
00221 
00222 /* BTF_FLIP is a "negation about -1", and is used to mark an integer j
00223  * that is normally non-negative.  BTF_FLIP (-1) is -1.  BTF_FLIP of
00224  * a number > -1 is negative, and BTF_FLIP of a number < -1 is positive.
00225  * BTF_FLIP (BTF_FLIP (j)) = j for all integers j.  UNFLIP (j) acts
00226  * like an "absolute value" operation, and is always >= -1.  You can test
00227  * whether or not an integer j is "flipped" with the BTF_ISFLIPPED (j)
00228  * macro.
00229  */
00230 
00231 #define BTF_FLIP(j) (-(j)-2)
00232 #define BTF_ISFLIPPED(j) ((j) < -1)
00233 #define BTF_UNFLIP(j) ((BTF_ISFLIPPED (j)) ? BTF_FLIP (j) : (j))
00234 
00235 /* ========================================================================== */
00236 /* === BTF version ========================================================== */
00237 /* ========================================================================== */
00238 
00239 /* All versions of BTF include these definitions.
00240  * As an example, to test if the version you are using is 1.2 or later:
00241  *
00242  *  if (BTF_VERSION >= BTF_VERSION_CODE (1,2)) ...
00243  *
00244  * This also works during compile-time:
00245  *
00246  *  #if (BTF >= BTF_VERSION_CODE (1,2))
00247  *      printf ("This is version 1.2 or later\n") ;
00248  *  #else
00249  *      printf ("This is an early version\n") ;
00250  *  #endif
00251  */
00252 
00253 #define BTF_DATE "May 31, 2007"
00254 #define BTF_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
00255 #define BTF_MAIN_VERSION 1
00256 #define BTF_SUB_VERSION 0
00257 #define BTF_SUBSUB_VERSION 0
00258 #define BTF_VERSION BTF_VERSION_CODE(BTF_MAIN_VERSION,BTF_SUB_VERSION)
00259 
00260 #ifdef __cplusplus
00261 }
00262 #endif
00263 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines