2
* -- SuperLU routine (version 2.0) --
3
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
4
* and Lawrence Berkeley National Lab.
11
extern int genmmd_(int *, int *, int *, int *, int *, int *, int *,
12
int *, int *, int *, int *, int *);
16
const int m, /* number of rows in matrix A. */
17
const int n, /* number of columns in matrix A. */
18
const int nnz,/* number of nonzeros in matrix A. */
19
int *colptr, /* column pointer of size n+1 for matrix A. */
20
int *rowind, /* row indices of size nz for matrix A. */
21
int *perm_c /* out - the column permutation vector. */
24
int Alen, *A, i, info, *p;
27
Alen = colamd_recommended(nnz, m, n);
29
if ( !(knobs = (double *) SUPERLU_MALLOC(COLAMD_KNOBS * sizeof(double))) )
30
ABORT("Malloc fails for knobs");
31
colamd_set_defaults(knobs);
33
if (!(A = (int *) SUPERLU_MALLOC(Alen * sizeof(int))) )
34
ABORT("Malloc fails for A[]");
35
if (!(p = (int *) SUPERLU_MALLOC((n+1) * sizeof(int))) )
36
ABORT("Malloc fails for p[]");
37
for (i = 0; i <= n; ++i) p[i] = colptr[i];
38
for (i = 0; i < nnz; ++i) A[i] = rowind[i];
39
info = colamd(m, n, Alen, A, p, knobs);
40
if ( info == FALSE ) ABORT("COLAMD failed");
42
for (i = 0; i < n; ++i) perm_c[p[i]] = i;
51
const int m, /* number of rows in matrix A. */
52
const int n, /* number of columns in matrix A. */
53
const int nz, /* number of nonzeros in matrix A */
54
int *colptr, /* column pointer of size n+1 for matrix A. */
55
int *rowind, /* row indices of size nz for matrix A. */
56
int *atanz, /* out - on exit, returns the actual number of
57
nonzeros in matrix A'*A. */
58
int **ata_colptr, /* out - size n+1 */
59
int **ata_rowind /* out - size *atanz */
65
* Form the structure of A'*A. A is an m-by-n matrix in column oriented
66
* format represented by (colptr, rowind). The output A'*A is in column
67
* oriented format (symmetrically, also row oriented), represented by
68
* (ata_colptr, ata_rowind).
70
* This routine is modified from GETATA routine by Tim Davis.
71
* The complexity of this algorithm is: SUM_{i=1,m} r(i)^2,
72
* i.e., the sum of the square of the row counts.
76
* o Do I need to withhold the *dense* rows?
77
* o How do I know the number of nonzeros in A'*A?
81
register int i, j, k, col, num_nz, ti, trow;
82
int *marker, *b_colptr, *b_rowind;
83
int *t_colptr, *t_rowind; /* a column oriented form of T = A' */
85
if ( !(marker = (int*) SUPERLU_MALLOC((SUPERLU_MAX(m,n)+1)*sizeof(int))) )
86
ABORT("SUPERLU_MALLOC fails for marker[]");
87
if ( !(t_colptr = (int*) SUPERLU_MALLOC((m+1) * sizeof(int))) )
88
ABORT("SUPERLU_MALLOC t_colptr[]");
89
if ( !(t_rowind = (int*) SUPERLU_MALLOC(nz * sizeof(int))) )
90
ABORT("SUPERLU_MALLOC fails for t_rowind[]");
93
/* Get counts of each column of T, and set up column pointers */
94
for (i = 0; i < m; ++i) marker[i] = 0;
95
for (j = 0; j < n; ++j) {
96
for (i = colptr[j]; i < colptr[j+1]; ++i)
100
for (i = 0; i < m; ++i) {
101
t_colptr[i+1] = t_colptr[i] + marker[i];
102
marker[i] = t_colptr[i];
105
/* Transpose the matrix from A to T */
106
for (j = 0; j < n; ++j)
107
for (i = colptr[j]; i < colptr[j+1]; ++i) {
109
t_rowind[marker[col]] = j;
114
/* ----------------------------------------------------------------
115
compute B = T * A, where column j of B is:
117
Struct (B_*j) = UNION ( Struct (T_*k) )
120
do not include the diagonal entry
122
( Partition A as: A = (A_*1, ..., A_*n)
123
Then B = T * A = (T * A_*1, ..., T * A_*n), where
124
T * A_*j = (T_*1, ..., T_*m) * A_*j. )
125
---------------------------------------------------------------- */
127
/* Zero the diagonal flag */
128
for (i = 0; i < n; ++i) marker[i] = -1;
130
/* First pass determines number of nonzeros in B */
132
for (j = 0; j < n; ++j) {
133
/* Flag the diagonal so it's not included in the B matrix */
136
for (i = colptr[j]; i < colptr[j+1]; ++i) {
137
/* A_kj is nonzero, add pattern of column T_*k to B_*j */
139
for (ti = t_colptr[k]; ti < t_colptr[k+1]; ++ti) {
141
if ( marker[trow] != j ) {
150
/* Allocate storage for A'*A */
151
if ( !(*ata_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
152
ABORT("SUPERLU_MALLOC fails for ata_colptr[]");
154
if ( !(*ata_rowind = (int*) SUPERLU_MALLOC( *atanz * sizeof(int)) ) )
155
ABORT("SUPERLU_MALLOC fails for ata_rowind[]");
157
b_colptr = *ata_colptr; /* aliasing */
158
b_rowind = *ata_rowind;
160
/* Zero the diagonal flag */
161
for (i = 0; i < n; ++i) marker[i] = -1;
163
/* Compute each column of B, one at a time */
165
for (j = 0; j < n; ++j) {
166
b_colptr[j] = num_nz;
168
/* Flag the diagonal so it's not included in the B matrix */
171
for (i = colptr[j]; i < colptr[j+1]; ++i) {
172
/* A_kj is nonzero, add pattern of column T_*k to B_*j */
174
for (ti = t_colptr[k]; ti < t_colptr[k+1]; ++ti) {
176
if ( marker[trow] != j ) {
178
b_rowind[num_nz++] = trow;
183
b_colptr[n] = num_nz;
185
SUPERLU_FREE(marker);
186
SUPERLU_FREE(t_colptr);
187
SUPERLU_FREE(t_rowind);
193
const int n, /* number of columns in matrix A. */
194
const int nz, /* number of nonzeros in matrix A */
195
int *colptr, /* column pointer of size n+1 for matrix A. */
196
int *rowind, /* row indices of size nz for matrix A. */
197
int *bnz, /* out - on exit, returns the actual number of
198
nonzeros in matrix A'*A. */
199
int **b_colptr, /* out - size n+1 */
200
int **b_rowind /* out - size *bnz */
207
* Form the structure of A'+A. A is an n-by-n matrix in column oriented
208
* format represented by (colptr, rowind). The output A'+A is in column
209
* oriented format (symmetrically, also row oriented), represented by
210
* (b_colptr, b_rowind).
213
register int i, j, k, col, num_nz;
214
int *t_colptr, *t_rowind; /* a column oriented form of T = A' */
217
if ( !(marker = (int*) SUPERLU_MALLOC( n * sizeof(int)) ) )
218
ABORT("SUPERLU_MALLOC fails for marker[]");
219
if ( !(t_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
220
ABORT("SUPERLU_MALLOC fails for t_colptr[]");
221
if ( !(t_rowind = (int*) SUPERLU_MALLOC( nz * sizeof(int)) ) )
222
ABORT("SUPERLU_MALLOC fails t_rowind[]");
225
/* Get counts of each column of T, and set up column pointers */
226
for (i = 0; i < n; ++i) marker[i] = 0;
227
for (j = 0; j < n; ++j) {
228
for (i = colptr[j]; i < colptr[j+1]; ++i)
232
for (i = 0; i < n; ++i) {
233
t_colptr[i+1] = t_colptr[i] + marker[i];
234
marker[i] = t_colptr[i];
237
/* Transpose the matrix from A to T */
238
for (j = 0; j < n; ++j)
239
for (i = colptr[j]; i < colptr[j+1]; ++i) {
241
t_rowind[marker[col]] = j;
246
/* ----------------------------------------------------------------
247
compute B = A + T, where column j of B is:
249
Struct (B_*j) = Struct (A_*k) UNION Struct (T_*k)
251
do not include the diagonal entry
252
---------------------------------------------------------------- */
254
/* Zero the diagonal flag */
255
for (i = 0; i < n; ++i) marker[i] = -1;
257
/* First pass determines number of nonzeros in B */
259
for (j = 0; j < n; ++j) {
260
/* Flag the diagonal so it's not included in the B matrix */
263
/* Add pattern of column A_*k to B_*j */
264
for (i = colptr[j]; i < colptr[j+1]; ++i) {
266
if ( marker[k] != j ) {
272
/* Add pattern of column T_*k to B_*j */
273
for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
275
if ( marker[k] != j ) {
283
/* Allocate storage for A+A' */
284
if ( !(*b_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
285
ABORT("SUPERLU_MALLOC fails for b_colptr[]");
287
if ( !(*b_rowind = (int*) SUPERLU_MALLOC( *bnz * sizeof(int)) ) )
288
ABORT("SUPERLU_MALLOC fails for b_rowind[]");
291
/* Zero the diagonal flag */
292
for (i = 0; i < n; ++i) marker[i] = -1;
294
/* Compute each column of B, one at a time */
296
for (j = 0; j < n; ++j) {
297
(*b_colptr)[j] = num_nz;
299
/* Flag the diagonal so it's not included in the B matrix */
302
/* Add pattern of column A_*k to B_*j */
303
for (i = colptr[j]; i < colptr[j+1]; ++i) {
305
if ( marker[k] != j ) {
307
(*b_rowind)[num_nz++] = k;
311
/* Add pattern of column T_*k to B_*j */
312
for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
314
if ( marker[k] != j ) {
316
(*b_rowind)[num_nz++] = k;
320
(*b_colptr)[n] = num_nz;
322
SUPERLU_FREE(marker);
323
SUPERLU_FREE(t_colptr);
324
SUPERLU_FREE(t_rowind);
328
get_perm_c(int ispec, SuperMatrix *A, int *perm_c)
333
* GET_PERM_C obtains a permutation matrix Pc, by applying the multiple
334
* minimum degree ordering code by Joseph Liu to matrix A'*A or A+A'.
335
* or using approximate minimum degree column ordering by Davis et. al.
336
* The LU factorization of A*Pc tends to have less fill than the LU
337
* factorization of A.
343
* Specifies the type of column ordering to reduce fill:
344
* = 1: minimum degree on the structure of A^T * A
345
* = 2: minimum degree on the structure of A^T + A
346
* = 3: approximate minimum degree for unsymmetric matrices
347
* If ispec == 0, the natural ordering (i.e., Pc = I) is returned.
349
* A (input) SuperMatrix*
350
* Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
351
* of the linear equations is A->nrow. Currently, the type of A
352
* can be: Stype = NC; Dtype = _D; Mtype = GE. In the future,
353
* more general A can be handled.
355
* perm_c (output) int*
356
* Column permutation vector of size A->ncol, which defines the
357
* permutation matrix Pc; perm_c[i] = j means column i of A is
358
* in position j in A*Pc.
362
NCformat *Astore = A->Store;
363
int m, n, bnz, *b_colptr, i;
364
int delta, maxint, nofsub, *invp;
365
int *b_rowind, *dhead, *qsize, *llist, *marker;
366
double t, SuperLU_timer_();
371
t = SuperLU_timer_();
373
case 0: /* Natural ordering */
374
for (i = 0; i < n; ++i) perm_c[i] = i;
375
printf("Use natural column ordering.\n");
377
case 1: /* Minimum degree ordering on A'*A */
378
getata(m, n, Astore->nnz, Astore->colptr, Astore->rowind,
379
&bnz, &b_colptr, &b_rowind);
380
printf("Use minimum degree ordering on A'*A.\n");
381
t = SuperLU_timer_() - t;
382
/*printf("Form A'*A time = %8.3f\n", t);*/
384
case 2: /* Minimum degree ordering on A'+A */
385
if ( m != n ) ABORT("Matrix is not square");
386
at_plus_a(n, Astore->nnz, Astore->colptr, Astore->rowind,
387
&bnz, &b_colptr, &b_rowind);
388
printf("Use minimum degree ordering on A'+A.\n");
389
t = SuperLU_timer_() - t;
390
/*printf("Form A'+A time = %8.3f\n", t);*/
392
case 3: /* Approximate minimum degree column ordering. */
393
get_colamd(m, n, Astore->nnz, Astore->colptr, Astore->rowind,
395
printf(".. Use approximate minimum degree column ordering.\n");
398
ABORT("Invalid ISPEC");
402
t = SuperLU_timer_();
404
/* Initialize and allocate storage for GENMMD. */
405
delta = 1; /* DELTA is a parameter to allow the choice of nodes
406
whose degree <= min-degree + DELTA. */
407
maxint = 2147483647; /* 2**31 - 1 */
408
invp = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
409
if ( !invp ) ABORT("SUPERLU_MALLOC fails for invp.");
410
dhead = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
411
if ( !dhead ) ABORT("SUPERLU_MALLOC fails for dhead.");
412
qsize = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
413
if ( !qsize ) ABORT("SUPERLU_MALLOC fails for qsize.");
414
llist = (int *) SUPERLU_MALLOC(n*sizeof(int));
415
if ( !llist ) ABORT("SUPERLU_MALLOC fails for llist.");
416
marker = (int *) SUPERLU_MALLOC(n*sizeof(int));
417
if ( !marker ) ABORT("SUPERLU_MALLOC fails for marker.");
419
/* Transform adjacency list into 1-based indexing required by GENMMD.*/
420
for (i = 0; i <= n; ++i) ++b_colptr[i];
421
for (i = 0; i < bnz; ++i) ++b_rowind[i];
423
genmmd_(&n, b_colptr, b_rowind, perm_c, invp, &delta, dhead,
424
qsize, llist, marker, &maxint, &nofsub);
426
/* Transform perm_c into 0-based indexing. */
427
for (i = 0; i < n; ++i) --perm_c[i];
429
SUPERLU_FREE(b_colptr);
430
SUPERLU_FREE(b_rowind);
435
SUPERLU_FREE(marker);
437
t = SuperLU_timer_() - t;
438
/* printf("call GENMMD time = %8.3f\n", t);*/
440
} else { /* Empty adjacency structure */
441
for (i = 0; i < n; ++i) perm_c[i] = i;
added
removed
Lib/sparse/SuperLU/SRC/get_perm_c.c
