Amesos Package Browser (Single Doxygen Collection) Development
amesos_btf_l_maxtrans.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === BTF_MAXTRANS ========================================================= */
00003 /* ========================================================================== */
00004 
00005 /* Finds a column permutation that maximizes the number of entries on the
00006  * diagonal of a sparse matrix.  See btf.h for more information.
00007  *
00008  * This function is identical to cs_maxtrans in CSparse, with the following
00009  * exceptions:
00010  *
00011  *  (1) cs_maxtrans finds both jmatch and imatch, where jmatch [i] = j and
00012  *  imatch [j] = i if row i is matched to column j.  This function returns
00013  *  just jmatch (the Match array).  The MATLAB interface to cs_maxtrans
00014  *  (the single-output cs_dmperm) returns imatch, not jmatch to the MATLAB
00015  *  caller.
00016  *
00017  *  (2) cs_maxtrans includes a pre-pass that counts the number of non-empty
00018  *  rows and columns (m2 and n2, respectively), and computes the matching
00019  *  using the transpose of A if m2 < n2.  cs_maxtrans also returns quickly
00020  *  if the diagonal of the matrix is already zero-free.  This pre-pass
00021  *  allows cs_maxtrans to be much faster than maxtrans, if the use of the
00022  *  transpose is warranted.
00023  *
00024  *  However, for square structurally non-singular matrices with one or more
00025  *  zeros on the diagonal, the pre-pass is a waste of time, and for these
00026  *  matrices, maxtrans can be twice as fast as cs_maxtrans.  Since the
00027  *  maxtrans function is intended primarily for square matrices that are
00028  *  typically structurally nonsingular, the pre-pass is not included here.
00029  *  If this maxtrans function is used on a matrix with many more columns
00030  *  than rows, consider passing the transpose to this function, or use
00031  *  cs_maxtrans instead.
00032  *
00033  *  (3) cs_maxtrans can operate as a randomized algorithm, to help avoid
00034  *  rare cases of excessive run-time.
00035  *
00036  *  (4) this maxtrans function includes an option that limits the total work
00037  *  performed.  If this limit is reached, the maximum transveral might not
00038  *  be found.
00039  *
00040  * Thus, for general usage, cs_maxtrans is preferred.  For square matrices that
00041  * are typically structurally non-singular, maxtrans is preferred.  A partial
00042  * maxtrans can still be very useful when solving a sparse linear system.
00043  *
00044  * Copyright (c) 2004-2007.  Tim Davis, University of Florida,
00045  * with support from Sandia National Laboratories.  All Rights Reserved.
00046  */ 
00047 
00048 /* This file should make the long int version of BTF */
00049 #define DLONG 1
00050 
00051 #include "amesos_btf_decl.h"
00052 #include "amesos_btf_internal.h"
00053 
00054 
00055 /* ========================================================================== */
00056 /* === augment ============================================================== */
00057 /* ========================================================================== */
00058 
00059 /* Perform a depth-first-search starting at column k, to find an augmenting
00060  * path.  An augmenting path is a sequence of row/column pairs (i1,k), (i2,j1),
00061  * (i3,j2), ..., (i(s+1), js), such that all of the following properties hold:
00062  *
00063  *  * column k is not matched to any row
00064  *  * entries in the path are nonzero
00065  *  * the pairs (i1,j1), (i2,j2), (i3,j3) ..., (is,js) have been 
00066  *      previously matched to each other
00067  *  * (i(s+1), js) is nonzero, and row i(s+1) is not matched to any column
00068  *
00069  * Once this path is found, the matching can be changed to the set of pairs
00070  * path.  An augmenting path is a sequence of row/column pairs
00071  *
00072  *  (i1,k), (i2,j1), (i3,j2), ..., (i(s+1), js)
00073  *
00074  * Once a row is matched with a column it remains matched with some column, but
00075  * not necessarily the column it was first matched with.
00076  *
00077  * In the worst case, this function can examine every nonzero in A.  Since it
00078  * is called n times by maxtrans, the total time of maxtrans can be as high as
00079  * O(n*nnz(A)).  To limit this work, pass a value of maxwork > 0.  Then at
00080  * most O((maxwork+1)*nnz(A)) work will be performed; the maximum matching might
00081  * not be found, however.
00082  *
00083  * This routine is very similar to the dfs routine in klu_kernel.c, in the
00084  * KLU sparse LU factorization package.  It is essentially identical to the
00085  * cs_augment routine in CSparse, and its recursive version (augment function
00086  * in cs_maxtransr_mex.c), except that this routine allows for the search to be
00087  * terminated early if too much work is being performed.
00088  *
00089  * The algorithm is based on the paper "On Algorithms for obtaining a maximum
00090  * transversal" by Iain Duff, ACM Trans. Mathematical Software, vol 7, no. 1,
00091  * pp. 315-330, and "Algorithm 575: Permutations for a zero-free diagonal",
00092  * same issue, pp. 387-390.  The code here is a new implementation of that
00093  * algorithm, with different data structures and control flow.  After writing
00094  * this code, I carefully compared my algorithm with MC21A/B (ACM Algorithm 575)
00095  * Some of the comparisons are partial because I didn't dig deeply into all of
00096  * the details of MC21A/B, such as how the stack is maintained.  The following
00097  * arguments are essentially identical between this code and MC21A:
00098  *
00099  * maxtrans MC21A,B
00100  * -------- -------
00101  * n    N     identical
00102  * k    JORD      identical
00103  * Ap   IP      column / row pointers
00104  * Ai   ICN     row / column indices
00105  * Ap[n]  LICN      length of index array (# of nonzeros in A)
00106  * Match  IPERM     output column / row permutation
00107  * nmatch NUMNZ     # of nonzeros on diagonal of permuted matrix
00108  * Flag   CV      mark a node as visited by the depth-first-search
00109  *
00110  * The following are different, but analogous:
00111  *
00112  * Cheap  ARP     indicates what part of the a column / row has
00113  *          already been matched.
00114  *
00115  * The following arguments are very different:
00116  *
00117  * -    LENR      # of entries in each row/column (unused in maxtrans)
00118  * Pstack OUT     Pstack keeps track of where we are in the depth-
00119  *          first-search scan of column j.  I think that OUT
00120  *          plays a similar role in MC21B, but I'm unsure.
00121  * Istack PR      keeps track of the rows in the path.  PR is a link
00122  *          list, though, whereas Istack is a stack.  Maxtrans
00123  *          does not use any link lists.
00124  * Jstack OUT? PR?    the stack for nodes in the path (unsure)
00125  *
00126  * The following control structures are roughly comparable:
00127  *
00128  * maxtrans     MC21B
00129  * --------     -----
00130  * for (k = 0 ; k < n ; k++)  DO 100 JORD=1,N
00131  * while (head >= 0)    DO 70 K=1,JORD
00132  * for (p = Cheap [j] ; ...)  DO 20 II=IN1,IN2
00133  * for (p = head ; ...)   DO 90 K=1,JORD
00134  */
00135 
00136 static Int amesos_augment
00137 (
00138     Int k,    /* which stage of the main loop we're in */
00139     Int Ap [ ],   /* column pointers, size n+1 */
00140     Int Ai [ ],   /* row indices, size nz = Ap [n] */
00141     Int Match [ ],  /* size n,  Match [i] = j if col j matched to i */
00142     Int Cheap [ ],  /* rows Ai [Ap [j] .. Cheap [j]-1] alread matched */
00143     Int Flag [ ], /* Flag [j] = k if j already visited this stage */
00144     Int Istack [ ], /* size n.  Row index stack. */
00145     Int Jstack [ ], /* size n.  Column index stack. */
00146     Int Pstack [ ], /* size n.  Keeps track of position in adjacency list */
00147     double *work, /* work performed by the depth-first-search */
00148     double maxwork  /* maximum work allowed */
00149 )
00150 {
00151     /* local variables, but "global" to all DFS levels: */
00152     Int found ; /* true if match found.  */
00153     Int head ;  /* top of stack */
00154 
00155     /* variables that are purely local to any one DFS level: */
00156     Int j2 ;  /* the next DFS goes to node j2 */
00157     Int pend ;  /* one past the end of the adjacency list for node j */
00158     Int pstart ;
00159     Int quick ;
00160 
00161     /* variables that need to be pushed then popped from the stack: */
00162     Int i ; /* the row tentatively matched to i if DFS successful */
00163     Int j ; /* the DFS is at the current node j */
00164     Int p ; /* current index into the adj. list for node j */
00165     /* the variables i, j, and p are stacked in Istack, Jstack, and Pstack */
00166 
00167     quick = (maxwork > 0) ;
00168 
00169     /* start a DFS to find a match for column k */
00170     found = FALSE ;
00171     i = EMPTY ;
00172     head = 0 ;
00173     Jstack [0] = k ;
00174     ASSERT (Flag [k] != k) ;
00175 
00176     while (head >= 0)
00177     {
00178   j = Jstack [head] ;
00179   pend = Ap [j+1] ;
00180 
00181   if (Flag [j] != k)      /* a node is not yet visited */
00182   {
00183 
00184       /* -------------------------------------------------------------- */
00185       /* prework for node j */
00186       /* -------------------------------------------------------------- */
00187 
00188       /* first time that j has been visited */
00189       Flag [j] = k ;
00190       /* cheap assignment: find the next unmatched row in col j.  This
00191        * loop takes at most O(nnz(A)) time for the sum total of all
00192        * calls to augment. */
00193       for (p = Cheap [j] ; p < pend && !found ; p++)
00194       {
00195     i = Ai [p] ;
00196     found = (Match [i] == EMPTY) ;
00197       }
00198       Cheap [j] = p ;
00199 
00200       /* -------------------------------------------------------------- */
00201 
00202       /* prepare for DFS */
00203       if (found)
00204       {
00205     /* end of augmenting path, column j matched with row i */
00206     Istack [head] = i ;
00207     break ;
00208       }
00209       /* set Pstack [head] to the first entry in column j to scan */
00210       Pstack [head] = Ap [j] ;
00211   }
00212 
00213   /* ------------------------------------------------------------------ */
00214   /* quick return if too much work done */
00215   /* ------------------------------------------------------------------ */
00216 
00217   if (quick && *work > maxwork)
00218   {
00219       /* too much work has been performed; abort the search */
00220       return (EMPTY) ;
00221   }
00222 
00223   /* ------------------------------------------------------------------ */
00224   /* DFS for nodes adjacent to j */
00225   /* ------------------------------------------------------------------ */
00226 
00227   /* If cheap assignment not made, continue the depth-first search.  All
00228    * rows in column j are already matched.  Add the adjacent nodes to the
00229    * stack by iterating through until finding another non-visited node.
00230    *
00231    * It is the following loop that can force maxtrans to take
00232    * O(n*nnz(A)) time. */
00233 
00234   pstart = Pstack [head] ;
00235   for (p = pstart ; p < pend ; p++)
00236   {
00237       i = Ai [p] ;
00238       j2 = Match [i] ;
00239       ASSERT (j2 != EMPTY) ;
00240       if (Flag [j2] != k)
00241       {
00242     /* Node j2 is not yet visited, start a depth-first search on
00243      * node j2.  Keep track of where we left off in the scan of adj
00244      * list of node j so we can restart j where we left off. */
00245     Pstack [head] = p + 1 ;
00246     /* Push j2 onto the stack and immediately break so we can
00247      * recurse on node j2.  Also keep track of row i which (if this
00248      * search for an augmenting path works) will be matched with the
00249      * current node j. */
00250     Istack [head] = i ;
00251     Jstack [++head] = j2 ;
00252     break ;
00253       }
00254   }
00255 
00256   /* ------------------------------------------------------------------ */
00257   /* determine how much work was just performed */
00258   /* ------------------------------------------------------------------ */
00259 
00260   *work += (p - pstart + 1) ;
00261 
00262   /* ------------------------------------------------------------------ */
00263   /* node j is done, but the postwork is postponed - see below */
00264   /* ------------------------------------------------------------------ */
00265 
00266   if (p == pend)
00267   {
00268       /* If all adjacent nodes of j are already visited, pop j from
00269        * stack and continue.  We failed to find a match. */
00270       head-- ;
00271   }
00272     }
00273 
00274     /* postwork for all nodes j in the stack */
00275     /* unwind the path and make the corresponding matches */
00276     if (found)
00277     {
00278   for (p = head ; p >= 0 ; p--)
00279   {
00280       j = Jstack [p] ;
00281       i = Istack [p] ;
00282 
00283       /* -------------------------------------------------------------- */
00284       /* postwork for node j */
00285       /* -------------------------------------------------------------- */
00286       /* if found, match row i with column j */
00287       Match [i] = j ;
00288   }
00289     }
00290     return (found) ;
00291 }
00292 
00293 
00294 /* ========================================================================== */
00295 /* === maxtrans ============================================================= */
00296 /* ========================================================================== */
00297 
00298 Int BTF(maxtrans)   /* returns # of columns in the matching */
00299 (
00300     /* --- input --- */
00301     Int nrow,     /* A is nrow-by-ncol in compressed column form */
00302     Int ncol,
00303     Int Ap [ ],     /* size ncol+1 */
00304     Int Ai [ ],     /* size nz = Ap [ncol] */
00305     double maxwork, /* do at most maxwork*nnz(A) work; no limit if <= 0.  This
00306          * work limit excludes the O(nnz(A)) cheap-match phase. */
00307 
00308     /* --- output --- */
00309     double *work,   /* work = -1 if maxwork > 0 and the total work performed
00310          * reached the maximum of maxwork*nnz(A)).
00311          * Otherwise, work = the total work performed. */
00312 
00313     Int Match [ ],  /* size nrow.  Match [i] = j if column j matched to row i */
00314 
00315     /* --- workspace --- */
00316     Int Work [ ]    /* size 5*ncol */
00317 )
00318 {
00319     Int *Cheap, *Flag, *Istack, *Jstack, *Pstack ;
00320     Int i, j, k, nmatch, work_limit_reached, result ;
00321 
00322     /* ---------------------------------------------------------------------- */
00323     /* get workspace and initialize */
00324     /* ---------------------------------------------------------------------- */
00325 
00326     Cheap  = Work ; Work += ncol ;
00327     Flag   = Work ; Work += ncol ;
00328 
00329     /* stack for non-recursive depth-first search in augment function */
00330     Istack = Work ; Work += ncol ;
00331     Jstack = Work ; Work += ncol ;
00332     Pstack = Work ;
00333 
00334     /* in column j, rows Ai [Ap [j] .. Cheap [j]-1] are known to be matched */
00335     for (j = 0 ; j < ncol ; j++)
00336     {
00337   Cheap [j] = Ap [j] ;
00338   Flag [j] = EMPTY ; 
00339     }
00340 
00341     /* all rows and columns are currently unmatched */
00342     for (i = 0 ; i < nrow ; i++)
00343     {
00344   Match [i] = EMPTY ;
00345     }
00346 
00347     if (maxwork > 0)
00348     {
00349   maxwork *= Ap [ncol] ;
00350     }
00351     *work = 0 ;
00352 
00353     /* ---------------------------------------------------------------------- */
00354     /* find a matching row for each column k */
00355     /* ---------------------------------------------------------------------- */
00356 
00357     nmatch = 0 ;
00358     work_limit_reached = FALSE ;
00359     for (k = 0 ; k < ncol ; k++)
00360     {
00361   /* find an augmenting path to match some row i to column k */
00362   result = amesos_augment (k, Ap, Ai, Match, Cheap, Flag, Istack, Jstack, Pstack,
00363       work, maxwork) ;
00364   if (result == TRUE)
00365   {
00366       /* we found it.  Match [i] = k for some row i has been done. */
00367       nmatch++ ;
00368   }
00369   else if (result == EMPTY)
00370   {
00371       /* augment gave up because of too much work, and no match found */
00372       work_limit_reached = TRUE ;
00373   }
00374     }
00375 
00376     /* ---------------------------------------------------------------------- */
00377     /* return the Match, and the # of matches made */
00378     /* ---------------------------------------------------------------------- */
00379 
00380     /* At this point, row i is matched to j = Match [i] if j >= 0.  i is an
00381      * unmatched row if Match [i] == EMPTY. */
00382 
00383     if (work_limit_reached)
00384     {
00385   /* return -1 if the work limit of maxwork*nnz(A) was reached */
00386   *work = EMPTY ;
00387     }
00388 
00389     return (nmatch) ;
00390 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines