1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
|
/*! \file
Copyright (c) 2003, The Regents of the University of California, through
Lawrence Berkeley National Laboratory (subject to receipt of any required
approvals from U.S. Dept. of Energy)
All rights reserved.
The source code is distributed under BSD license, see the file License.txt
at the top-level directory.
*/
/*! @file
* \brief Finds a row permutation so that the matrix has large entries on the diagonal
*
* <pre>
* -- SuperLU routine (version 4.0) --
* Lawrence Berkeley National Laboratory.
* June 30, 2009
* </pre>
*/
#include "slu_cdefs.h"
extern int_t mc64id_(int_t*);
extern int_t mc64ad_(int_t*, int_t*, int_t*, int_t [], int_t [], double [],
int_t*, int_t [], int_t*, int_t[], int_t*, double [],
int_t [], int_t []);
/*! \brief
*
* <pre>
* Purpose
* =======
*
* CLDPERM finds a row permutation so that the matrix has large
* entries on the diagonal.
*
* Arguments
* =========
*
* job (input) int
* Control the action. Possible values for JOB are:
* = 1 : Compute a row permutation of the matrix so that the
* permuted matrix has as many entries on its diagonal as
* possible. The values on the diagonal are of arbitrary size.
* HSL subroutine MC21A/AD is used for this.
* = 2 : Compute a row permutation of the matrix so that the smallest
* value on the diagonal of the permuted matrix is maximized.
* = 3 : Compute a row permutation of the matrix so that the smallest
* value on the diagonal of the permuted matrix is maximized.
* The algorithm differs from the one used for JOB = 2 and may
* have quite a different performance.
* = 4 : Compute a row permutation of the matrix so that the sum
* of the diagonal entries of the permuted matrix is maximized.
* = 5 : Compute a row permutation of the matrix so that the product
* of the diagonal entries of the permuted matrix is maximized
* and vectors to scale the matrix so that the nonzero diagonal
* entries of the permuted matrix are one in absolute value and
* all the off-diagonal entries are less than or equal to one in
* absolute value.
* Restriction: 1 <= JOB <= 5.
*
* n (input) int
* The order of the matrix.
*
* nnz (input) int
* The number of nonzeros in the matrix.
*
* adjncy (input) int*, of size nnz
* The adjacency structure of the matrix, which contains the row
* indices of the nonzeros.
*
* colptr (input) int*, of size n+1
* The pointers to the beginning of each column in ADJNCY.
*
* nzval (input) complex*, of size nnz
* The nonzero values of the matrix. nzval[k] is the value of
* the entry corresponding to adjncy[k].
* It is not used if job = 1.
*
* perm (output) int*, of size n
* The permutation vector. perm[i] = j means row i in the
* original matrix is in row j of the permuted matrix.
*
* u (output) double*, of size n
* If job = 5, the natural logarithms of the row scaling factors.
*
* v (output) double*, of size n
* If job = 5, the natural logarithms of the column scaling factors.
* The scaled matrix B has entries b_ij = a_ij * exp(u_i + v_j).
* </pre>
*/
int
cldperm(int_t job, int_t n, int_t nnz, int_t colptr[], int_t adjncy[],
complex nzval[], int_t *perm, float u[], float v[])
{
int_t i, liw, ldw, num;
int_t *iw, icntl[10], info[10];
double *dw;
double *nzval_d = (double *) SUPERLU_MALLOC(nnz * sizeof(double));
#if ( DEBUGlevel>=1 )
CHECK_MALLOC("Enter cldperm()");
#endif
liw = 5*n;
if ( job == 3 ) liw = 10*n + nnz;
if ( !(iw = intMalloc(liw)) ) ABORT("Malloc fails for iw[]");
ldw = 3*n + nnz;
if ( !(dw = (double*) SUPERLU_MALLOC(ldw * sizeof(double))) )
ABORT("Malloc fails for dw[]");
/* Increment one to get 1-based indexing. */
for (i = 0; i <= n; ++i) ++colptr[i];
for (i = 0; i < nnz; ++i) ++adjncy[i];
#if ( DEBUGlevel>=2 )
printf("LDPERM(): n %d, nnz %d\n", n, nnz);
slu_PrintInt10("colptr", n+1, colptr);
slu_PrintInt10("adjncy", nnz, adjncy);
#endif
/*
* NOTE:
* =====
*
* MC64AD assumes that column permutation vector is defined as:
* perm(i) = j means column i of permuted A is in column j of original A.
*
* Since a symmetric permutation preserves the diagonal entries. Then
* by the following relation:
* P'(A*P')P = P'A
* we can apply inverse(perm) to rows of A to get large diagonal entries.
* But, since 'perm' defined in MC64AD happens to be the reverse of
* SuperLU's definition of permutation vector, therefore, it is already
* an inverse for our purpose. We will thus use it directly.
*
*/
mc64id_(icntl);
#if 0
/* Suppress error and warning messages. */
icntl[0] = -1;
icntl[1] = -1;
#endif
for (i = 0; i < nnz; ++i) nzval_d[i] = c_abs1(&nzval[i]);
mc64ad_(&job, &n, &nnz, colptr, adjncy, nzval_d, &num, perm,
&liw, iw, &ldw, dw, icntl, info);
#if ( DEBUGlevel>=2 )
slu_PrintInt10("perm", n, perm);
printf(".. After MC64AD info %d\tsize of matching %d\n", info[0], num);
#endif
if ( info[0] == 1 ) { /* Structurally singular */
printf(".. The last %d permutations:\n", n-num);
slu_PrintInt10("perm", n-num, &perm[num]);
}
/* Restore to 0-based indexing. */
for (i = 0; i <= n; ++i) --colptr[i];
for (i = 0; i < nnz; ++i) --adjncy[i];
for (i = 0; i < n; ++i) --perm[i];
if ( job == 5 )
for (i = 0; i < n; ++i) {
u[i] = dw[i];
v[i] = dw[n+i];
}
SUPERLU_FREE(iw);
SUPERLU_FREE(dw);
SUPERLU_FREE(nzval_d);
#if ( DEBUGlevel>=1 )
CHECK_MALLOC("Exit cldperm()");
#endif
return info[0];
}
|
