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