Amesos Package Browser (Single Doxygen Collection) Development
amesos_btf_l_order.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === BTF_ORDER ============================================================ */
00003 /* ========================================================================== */
00004 
00005 /* Find a permutation P and Q to permute a square sparse matrix into upper block
00006  * triangular form.  A(P,Q) will contain a zero-free diagonal if A has
00007  * structural full-rank.  Otherwise, the number of nonzeros on the diagonal of
00008  * A(P,Q) will be maximized, and will equal the structural rank of A.
00009  *
00010  * Q[k] will be "flipped" if a zero-free diagonal was not found.  Q[k] will be
00011  * negative, and j = BTF_UNFLIP (Q [k]) gives the corresponding permutation.
00012  *
00013  * R defines the block boundaries of A(P,Q).  The kth block consists of rows
00014  * and columns R[k] to R[k+1]-1.
00015  *
00016  * If maxwork > 0 on input, then the work performed in btf_maxtrans is limited
00017  * to maxwork*nnz(A) (excluding the "cheap match" phase, which can take another
00018  * nnz(A) work).  On output, the work parameter gives the actual work performed,
00019  * or -1 if the limit was reached.  In the latter case, the diagonal of A(P,Q)
00020  * might not be zero-free, and the number of nonzeros on the diagonal of A(P,Q)
00021  * might not be equal to the structural rank.
00022  *
00023  * See btf.h for more details.
00024  *
00025  * Copyright (c) 2004-2007.  Tim Davis, University of Florida,
00026  * with support from Sandia National Laboratories.  All Rights Reserved.
00027  */
00028 
00029 /* This file should make the long int version of BTF */
00030 #define DLONG 1
00031 
00032 #include "amesos_btf_decl.h"
00033 #include "amesos_btf_internal.h"
00034 
00035 /* This function only operates on square matrices (either structurally full-
00036  * rank, or structurally rank deficient). */
00037 
00038 Int BTF(order)      /* returns number of blocks found */
00039 (
00040     /* input, not modified: */
00041     Int n,      /* A is n-by-n in compressed column form */
00042     Int Ap [ ],     /* size n+1 */
00043     Int Ai [ ],     /* size nz = Ap [n] */
00044     double maxwork, /* do at most maxwork*nnz(A) work in the maximum
00045          * transversal; no limit if <= 0 */
00046 
00047     /* output, not defined on input */
00048     double *work,   /* work performed in maxtrans, or -1 if limit reached */
00049     Int P [ ],      /* size n, row permutation */
00050     Int Q [ ],      /* size n, column permutation */
00051     Int R [ ],      /* size n+1.  block b is in rows/cols R[b] ... R[b+1]-1 */
00052     Int *nmatch,    /* # nonzeros on diagonal of P*A*Q */
00053 
00054     /* workspace, not defined on input or output */
00055     Int Work [ ]    /* size 5n */
00056 )
00057 {
00058     Int *Flag ;
00059     Int nblocks, i, j, nbadcol ;
00060 
00061     /* ---------------------------------------------------------------------- */
00062     /* compute the maximum matching */
00063     /* ---------------------------------------------------------------------- */
00064 
00065     /* if maxwork > 0, then a maximum matching might not be found */
00066 
00067     *nmatch = BTF(maxtrans) (n, n, Ap, Ai, maxwork, work, Q, Work) ;
00068 
00069     /* ---------------------------------------------------------------------- */
00070     /* complete permutation if the matrix is structurally singular */
00071     /* ---------------------------------------------------------------------- */
00072 
00073     /* Since the matrix is square, ensure BTF_UNFLIP(Q[0..n-1]) is a
00074      * permutation of the columns of A so that A has as many nonzeros on the
00075      * diagonal as possible.
00076      */
00077 
00078     if (*nmatch < n)
00079     {
00080   /* get a size-n work array */
00081   Flag = Work + n ;
00082   for (j = 0 ; j < n ; j++)
00083   {
00084       Flag [j] = 0 ;
00085   }
00086 
00087   /* flag all matched columns */
00088   for (i = 0 ; i < n ; i++)
00089   {
00090       j = Q [i] ;
00091       if (j != EMPTY)
00092       {
00093     /* row i and column j are matched to each other */
00094     Flag [j] = 1 ;
00095       }
00096   }
00097 
00098   /* make a list of all unmatched columns, in Work [0..nbadcol-1]  */
00099   nbadcol = 0 ;
00100   for (j = n-1 ; j >= 0 ; j--)
00101   {
00102       if (!Flag [j])
00103       {
00104     /* j is matched to nobody */
00105     Work [nbadcol++] = j ;
00106       }
00107   }
00108   ASSERT (*nmatch + nbadcol == n) ;
00109 
00110   /* make an assignment for each unmatched row */
00111   for (i = 0 ; i < n ; i++)
00112   {
00113       if (Q [i] == EMPTY && nbadcol > 0)
00114       {
00115     /* get an unmatched column j */
00116     j = Work [--nbadcol] ;
00117     /* assign j to row i and flag the entry by "flipping" it */
00118     Q [i] = BTF_FLIP (j) ;
00119       }
00120   }
00121     }
00122 
00123     /* The permutation of a square matrix can be recovered as follows: Row i is
00124      * matched with column j, where j = BTF_UNFLIP (Q [i]) and where j
00125      * will always be in the valid range 0 to n-1.  The entry A(i,j) is zero
00126      * if BTF_ISFLIPPED (Q [i]) is true, and nonzero otherwise.  nmatch
00127      * is the number of entries in the Q array that are non-negative. */
00128 
00129     /* ---------------------------------------------------------------------- */
00130     /* find the strongly connected components */
00131     /* ---------------------------------------------------------------------- */
00132 
00133     nblocks = BTF(strongcomp) (n, Ap, Ai, Q, P, R, Work) ;
00134     return (nblocks) ;
00135 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines