Amesos Package Browser (Single Doxygen Collection)
Development
src
SuiteSparse
BTF
Source
amesos_btf_l_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
/* This file should make the long int version of BTF */
30
#define DLONG 1
31
32
#include "
amesos_btf_decl.h
"
33
#include "
amesos_btf_internal.h
"
34
35
/* This function only operates on square matrices (either structurally full-
36
* rank, or structurally rank deficient). */
37
38
Int
BTF
(order)
/* returns number of blocks found */
39
(
40
/* input, not modified: */
41
Int
n
,
/* A is n-by-n in compressed column form */
42
Int
Ap [ ],
/* size n+1 */
43
Int
Ai [ ],
/* size nz = Ap [n] */
44
double
maxwork,
/* do at most maxwork*nnz(A) work in the maximum
45
* transversal; no limit if <= 0 */
46
47
/* output, not defined on input */
48
double
*work,
/* work performed in maxtrans, or -1 if limit reached */
49
Int
P
[ ],
/* size n, row permutation */
50
Int
Q [ ],
/* size n, column permutation */
51
Int
R [ ],
/* size n+1. block b is in rows/cols R[b] ... R[b+1]-1 */
52
Int
*nmatch,
/* # nonzeros on diagonal of P*A*Q */
53
54
/* workspace, not defined on input or output */
55
Int
Work [ ]
/* size 5n */
56
)
57
{
58
Int
*Flag ;
59
Int
nblocks, i, j, nbadcol ;
60
61
/* ---------------------------------------------------------------------- */
62
/* compute the maximum matching */
63
/* ---------------------------------------------------------------------- */
64
65
/* if maxwork > 0, then a maximum matching might not be found */
66
67
*nmatch =
BTF
(maxtrans) (
n
,
n
, Ap, Ai, maxwork, work, Q, Work) ;
68
69
/* ---------------------------------------------------------------------- */
70
/* complete permutation if the matrix is structurally singular */
71
/* ---------------------------------------------------------------------- */
72
73
/* Since the matrix is square, ensure BTF_UNFLIP(Q[0..n-1]) is a
74
* permutation of the columns of A so that A has as many nonzeros on the
75
* diagonal as possible.
76
*/
77
78
if
(*nmatch <
n
)
79
{
80
/* get a size-n work array */
81
Flag = Work +
n
;
82
for
(j = 0 ; j <
n
; j++)
83
{
84
Flag [j] = 0 ;
85
}
86
87
/* flag all matched columns */
88
for
(i = 0 ; i <
n
; i++)
89
{
90
j = Q [i] ;
91
if
(j !=
EMPTY
)
92
{
93
/* row i and column j are matched to each other */
94
Flag [j] = 1 ;
95
}
96
}
97
98
/* make a list of all unmatched columns, in Work [0..nbadcol-1] */
99
nbadcol = 0 ;
100
for
(j =
n
-1 ; j >= 0 ; j--)
101
{
102
if
(!Flag [j])
103
{
104
/* j is matched to nobody */
105
Work [nbadcol++] = j ;
106
}
107
}
108
ASSERT
(*nmatch + nbadcol ==
n
) ;
109
110
/* make an assignment for each unmatched row */
111
for
(i = 0 ; i <
n
; i++)
112
{
113
if
(Q [i] ==
EMPTY
&& nbadcol > 0)
114
{
115
/* get an unmatched column j */
116
j = Work [--nbadcol] ;
117
/* assign j to row i and flag the entry by "flipping" it */
118
Q [i] =
BTF_FLIP
(j) ;
119
}
120
}
121
}
122
123
/* The permutation of a square matrix can be recovered as follows: Row i is
124
* matched with column j, where j = BTF_UNFLIP (Q [i]) and where j
125
* will always be in the valid range 0 to n-1. The entry A(i,j) is zero
126
* if BTF_ISFLIPPED (Q [i]) is true, and nonzero otherwise. nmatch
127
* is the number of entries in the Q array that are non-negative. */
128
129
/* ---------------------------------------------------------------------- */
130
/* find the strongly connected components */
131
/* ---------------------------------------------------------------------- */
132
133
nblocks =
BTF
(strongcomp) (
n
, Ap, Ai, Q,
P
, R, Work) ;
134
return
(nblocks) ;
135
}
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_FLIP
#define BTF_FLIP(j)
Definition:
amesos_btf_decl.h:231
amesos_btf_decl.h
BTF
Int BTF(order)
Definition:
amesos_btf_l_order.c:38
n
int n
Generated by
1.8.14