Amesos Package Browser (Single Doxygen Collection) Development
amesos_btf_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 #include "amesos_btf_decl.h"
00030 #include "amesos_btf_internal.h"
00031
00032 /* This function only operates on square matrices (either structurally full-
00033  * rank, or structurally rank deficient). */
00034
00035 Int BTF(order)      /* returns number of blocks found */
00036 (
00037     /* input, not modified: */
00038     Int n,      /* A is n-by-n in compressed column form */
00039     Int Ap [ ],     /* size n+1 */
00040     Int Ai [ ],     /* size nz = Ap [n] */
00041     double maxwork, /* do at most maxwork*nnz(A) work in the maximum
00042          * transversal; no limit if <= 0 */
00043
00044     /* output, not defined on input */
00045     double *work,   /* work performed in maxtrans, or -1 if limit reached */
00046     Int P [ ],      /* size n, row permutation */
00047     Int Q [ ],      /* size n, column permutation */
00048     Int R [ ],      /* size n+1.  block b is in rows/cols R[b] ... R[b+1]-1 */
00049     Int *nmatch,    /* # nonzeros on diagonal of P*A*Q */
00050
00051     /* workspace, not defined on input or output */
00052     Int Work [ ]    /* size 5n */
00053 )
00054 {
00055     Int *Flag ;
00056     Int nblocks, i, j, nbadcol ;
00057
00058     /* ---------------------------------------------------------------------- */
00059     /* compute the maximum matching */
00060     /* ---------------------------------------------------------------------- */
00061
00062     /* if maxwork > 0, then a maximum matching might not be found */
00063
00064     *nmatch = BTF(maxtrans) (n, n, Ap, Ai, maxwork, work, Q, Work) ;
00065
00066     /* ---------------------------------------------------------------------- */
00067     /* complete permutation if the matrix is structurally singular */
00068     /* ---------------------------------------------------------------------- */
00069
00070     /* Since the matrix is square, ensure BTF_UNFLIP(Q[0..n-1]) is a
00071      * permutation of the columns of A so that A has as many nonzeros on the
00072      * diagonal as possible.
00073      */
00074
00075     if (*nmatch < n)
00076     {
00077   /* get a size-n work array */
00078   Flag = Work + n ;
00079   for (j = 0 ; j < n ; j++)
00080   {
00081       Flag [j] = 0 ;
00082   }
00083
00084   /* flag all matched columns */
00085   for (i = 0 ; i < n ; i++)
00086   {
00087       j = Q [i] ;
00088       if (j != EMPTY)
00089       {
00090     /* row i and column j are matched to each other */
00091     Flag [j] = 1 ;
00092       }
00093   }
00094
00095   /* make a list of all unmatched columns, in Work [0..nbadcol-1]  */
00096   nbadcol = 0 ;
00097   for (j = n-1 ; j >= 0 ; j--)
00098   {
00099       if (!Flag [j])
00100       {
00101     /* j is matched to nobody */
00102     Work [nbadcol++] = j ;
00103       }
00104   }
00105   ASSERT (*nmatch + nbadcol == n) ;
00106
00107   /* make an assignment for each unmatched row */
00108   for (i = 0 ; i < n ; i++)
00109   {
00110       if (Q [i] == EMPTY && nbadcol > 0)
00111       {
00112     /* get an unmatched column j */
00113     j = Work [--nbadcol] ;
00114     /* assign j to row i and flag the entry by "flipping" it */
00115     Q [i] = BTF_FLIP (j) ;
00116       }
00117   }
00118     }
00119
00120     /* The permutation of a square matrix can be recovered as follows: Row i is
00121      * matched with column j, where j = BTF_UNFLIP (Q [i]) and where j
00122      * will always be in the valid range 0 to n-1.  The entry A(i,j) is zero
00123      * if BTF_ISFLIPPED (Q [i]) is true, and nonzero otherwise.  nmatch
00124      * is the number of entries in the Q array that are non-negative. */
00125
00126     /* ---------------------------------------------------------------------- */
00127     /* find the strongly connected components */
00128     /* ---------------------------------------------------------------------- */
00129
00130     nblocks = BTF(strongcomp) (n, Ap, Ai, Q, P, R, Work) ;
00131     return (nblocks) ;
00132 }
```