Amesos Package Browser (Single Doxygen Collection)
Development
src
SuiteSparse
BTF
Source
amesos_btf_order.c
Go to the documentation of this file.
1
/* ========================================================================== */
2
/* === BTF_ORDER ============================================================ */
3
/* ========================================================================== */
4
5
/* Find a permutation P and Q to permute a square sparse matrix into upper block
6
* triangular form. A(P,Q) will contain a zero-free diagonal if A has
7
* structural full-rank. Otherwise, the number of nonzeros on the diagonal of
8
* A(P,Q) will be maximized, and will equal the structural rank of A.
9
*
10
* Q[k] will be "flipped" if a zero-free diagonal was not found. Q[k] will be
11
* negative, and j = BTF_UNFLIP (Q [k]) gives the corresponding permutation.
12
*
13
* R defines the block boundaries of A(P,Q). The kth block consists of rows
14
* and columns R[k] to R[k+1]-1.
15
*
16
* If maxwork > 0 on input, then the work performed in btf_maxtrans is limited
17
* to maxwork*nnz(A) (excluding the "cheap match" phase, which can take another
18
* nnz(A) work). On output, the work parameter gives the actual work performed,
19
* or -1 if the limit was reached. In the latter case, the diagonal of A(P,Q)
20
* might not be zero-free, and the number of nonzeros on the diagonal of A(P,Q)
21
* might not be equal to the structural rank.
22
*
23
* See btf.h for more details.
24
*
25
* Copyright (c) 2004-2007. Tim Davis, University of Florida,
26
* with support from Sandia National Laboratories. All Rights Reserved.
27
*/
28
29
#include "
amesos_btf_decl.h
"
30
#include "
amesos_btf_internal.h
"
31
32
/* This function only operates on square matrices (either structurally full-
33
* rank, or structurally rank deficient). */
34
35
Int
BTF
(order)
/* returns number of blocks found */
36
(
37
/* input, not modified: */
38
Int
n,
/* A is n-by-n in compressed column form */
39
Int
Ap [ ],
/* size n+1 */
40
Int
Ai [ ],
/* size nz = Ap [n] */
41
double
maxwork,
/* do at most maxwork*nnz(A) work in the maximum
42
* transversal; no limit if <= 0 */
43
44
/* output, not defined on input */
45
double
*work,
/* work performed in maxtrans, or -1 if limit reached */
46
Int
P
[ ],
/* size n, row permutation */
47
Int
Q [ ],
/* size n, column permutation */
48
Int
R [ ],
/* size n+1. block b is in rows/cols R[b] ... R[b+1]-1 */
49
Int
*nmatch,
/* # nonzeros on diagonal of P*A*Q */
50
51
/* workspace, not defined on input or output */
52
Int
Work [ ]
/* size 5n */
53
)
54
{
55
Int
*Flag ;
56
Int
nblocks, i, j, nbadcol ;
57
58
/* ---------------------------------------------------------------------- */
59
/* compute the maximum matching */
60
/* ---------------------------------------------------------------------- */
61
62
/* if maxwork > 0, then a maximum matching might not be found */
63
64
*nmatch =
BTF
(maxtrans) (n, n, Ap, Ai, maxwork, work, Q, Work) ;
65
66
/* ---------------------------------------------------------------------- */
67
/* complete permutation if the matrix is structurally singular */
68
/* ---------------------------------------------------------------------- */
69
70
/* Since the matrix is square, ensure BTF_UNFLIP(Q[0..n-1]) is a
71
* permutation of the columns of A so that A has as many nonzeros on the
72
* diagonal as possible.
73
*/
74
75
if
(*nmatch < n)
76
{
77
/* get a size-n work array */
78
Flag = Work + n ;
79
for
(j = 0 ; j < n ; j++)
80
{
81
Flag [j] = 0 ;
82
}
83
84
/* flag all matched columns */
85
for
(i = 0 ; i < n ; i++)
86
{
87
j = Q [i] ;
88
if
(j !=
EMPTY
)
89
{
90
/* row i and column j are matched to each other */
91
Flag [j] = 1 ;
92
}
93
}
94
95
/* make a list of all unmatched columns, in Work [0..nbadcol-1] */
96
nbadcol = 0 ;
97
for
(j = n-1 ; j >= 0 ; j--)
98
{
99
if
(!Flag [j])
100
{
101
/* j is matched to nobody */
102
Work [nbadcol++] = j ;
103
}
104
}
105
ASSERT
(*nmatch + nbadcol == n) ;
106
107
/* make an assignment for each unmatched row */
108
for
(i = 0 ; i < n ; i++)
109
{
110
if
(Q [i] ==
EMPTY
&& nbadcol > 0)
111
{
112
/* get an unmatched column j */
113
j = Work [--nbadcol] ;
114
/* assign j to row i and flag the entry by "flipping" it */
115
Q [i] =
BTF_FLIP
(j) ;
116
}
117
}
118
}
119
120
/* The permutation of a square matrix can be recovered as follows: Row i is
121
* matched with column j, where j = BTF_UNFLIP (Q [i]) and where j
122
* will always be in the valid range 0 to n-1. The entry A(i,j) is zero
123
* if BTF_ISFLIPPED (Q [i]) is true, and nonzero otherwise. nmatch
124
* is the number of entries in the Q array that are non-negative. */
125
126
/* ---------------------------------------------------------------------- */
127
/* find the strongly connected components */
128
/* ---------------------------------------------------------------------- */
129
130
nblocks =
BTF
(strongcomp) (n, Ap, Ai, Q,
P
, R, Work) ;
131
return
(nblocks) ;
132
}
EMPTY
#define EMPTY
Definition:
amesos_amd_internal.h:144
Int
#define Int
Definition:
amesos_amd_internal.h:190
amesos_btf_internal.h
P
#define P(k)
Definition:
amesos_cholmod_solve.c:76
ASSERT
#define ASSERT(expression)
Definition:
amesos_amd_internal.h:331
BTF
Int BTF(order)
Definition:
amesos_btf_order.c:35
BTF_FLIP
#define BTF_FLIP(j)
Definition:
amesos_btf_decl.h:231
amesos_btf_decl.h
Generated by
1.8.14