1 #ifndef lint 2 static char vcid[] = "$Id: zerodiag.c,v 1.1 1994/03/18 00:27:07 gropp Exp $"; 3 #endif 4 5 /* 6 This file contains routines to reorder a matrix so that the diagonal 7 elements meet a simple criteria. The most common use is expected to 8 be reordering matrices with zero diagonal elements for incomplete 9 factorizations. A pivoting (partial or pairwise) routine is planned. 10 */ 11 12 #include "tools.h" 13 #include <math.h> 14 #include "sparse/spmat.h" 15 #include "sparse/sppriv.h" 16 #define SWAP(a,b) {int _t; _t = a; a = b; b = _t; } 17 #define SWAPD(a,b) {double _t; _t = a; a = b; b = _t; } 18 19 int SpiZeroFindPre( ); 20 21 /*@ 22 SpUnSymmetricReorderForZeroDiagonal - Reorder the matrix so that no 23 zeros (or small elements) are on the diagonal 24 25 Input Parameters: 26 . mat - matrix to reorder 27 . atol - elements smaller than this in magnitude are candidates for 28 reordering 29 . rmap,cmap - row and column permutations. On entrance, they should be 30 either the identity ordering or copies of the matrix's 31 row and column permutations. On exit, they will be the 32 modified permutations. 33 34 Error Conditions: 35 $ No memory 36 $ Unable to find suitable reording 37 38 Notes: 39 This is not intended as a replacement for pivoting for matrices that 40 have ``bad'' structure; use the pivoting factorization routine in that 41 case. 42 43 After this routine has been used, you MUST recompute the icolmap 44 mapping (this is the inverse column mapping; since this routine 45 modifies the column mapping, the inverse column mapping becomes 46 invalid. The reason that it is not "corrected" here is because 47 additional orderings can be applied before the inverse column 48 mapping is needed. THIS DECISION MAY CHANGE if it proves a poor 49 choice. Use the routine SpInverse to recompute the inverse column mapping. 50 51 Algorithm: 52 Column pivoting is used. Choice of column is made by looking at the 53 non-zero elements in the row. This algorithm is simple and fast but 54 does NOT guarentee that a non-singular or well conditioned 55 principle submatrix will be produced. 56 @*/ 57 void SpUnSymmetricReorderForZeroDiagonal( mat, atol, rmap, cmap ) 58 SpMat *mat; 59 double atol; 60 int *rmap, *cmap; 61 { 62 int prow, k, nz, n, repl, *j, *col, *row; 63 double *v, repla; 64 65 SPLITTOMAT(mat); 66 67 col = cmap; 68 row = rmap; 69 n = mat->rows; 70 71 for (prow=0; prow<n; prow++) { 72 SpScatterFromRow( mat, row[prow], &nz, &j, &v ); 73 for (k=0; k<nz; k++) 74 if (col[j[k]] == prow) break; 75 if (k >= nz || fabs(v[k]) <= atol) { 76 /* Element too small or zero; find the best candidate */ 77 repl = prow; 78 repla = (k >= nz) ? 0.0 : fabs(v[k]); 79 for (k=0; k<nz; k++) 80 if (col[j[k]] > prow && fabs(v[k]) > repla) { 81 repl = col[j[k]]; 82 repla = fabs(v[k]); 83 } 84 if (prow == repl) { 85 /* Now we need to look for an element that allows us 86 to pivot with a previous column. To do this, we need 87 to be sure that we don't introduce a zero in a previous 88 diagonal */ 89 if (!SpiZeroFindPre( mat, prow, row, col, repla, atol, 90 &repl, &repla )) { 91 SETERRC(1,"Can not reorder matrix"); 92 break; 93 } 94 } 95 SWAP(col[prow],col[repl]); 96 } 97 } 98 } 99 100 /*@ 101 SpSymmetricReorderForZeroDiagonal - Reorder the matrix so that no 102 zeros (or small elements) are on the diagonal 103 104 Input Parameters: 105 . mat - matrix to reorder 106 . atol - elements smaller than this in magnitude are candidates for 107 reordering 108 . rmap,cmap - row and column permutations. On entrance, they should be 109 either the identity ordering or copies of the matrix's 110 row and column permutations. On exit, they will be the 111 moidified permutations. 112 113 Error Conditions: 114 $ No memory 115 $ Unable to find suitable reording 116 117 Notes: 118 This is not intended as a replacement for pivoting for matrices that 119 have ``bad'' structure; use the pivoting factorization routine in that 120 case (those routines are not yet available). 121 122 After this routine has been used, you MUST recompute the icolmap 123 mapping (this is the inverse column mapping; since this routine 124 modifies the column mapping, the inverse column mapping becomes 125 invalid. The reason that it is not "corrected" here is because 126 additional orderings can be applied before the inverse column 127 mapping is needed. THIS DECISION MAY CHANGE if it proves a poor 128 choice. Use the routine SpInverse to recompute the inverse column 129 mapping. 130 131 Algorithm: 132 This is more complicated than the unsymmetric version, since a zero 133 diagonal element can not be moved off the diagonal by symmetric 134 permutations. 135 Rather, we try to insure that a zero diagonal will be subject to some 136 fill before it is encountered. We do this as follows 137 138 For each zero diagonal element, look down that row for a non-zero entry 139 such that the diagonal under (in the same column) is non-zero. The 140 criteria is to pick the largest corresponding diagonal (other criteria 141 are possible, such as a balance or a "no-fill" estimate of the 142 resulting pivot sizes). To speed this computation up, we gather some 143 auxillery data: location of the diagonal elements and their values. 144 We have to update these structures as we choose a reordering. 145 @*/ 146 void SpSymmetricReorderForZeroDiagonal( mat, atol, rmap, cmap ) 147 SpMat *mat; 148 double atol; 149 int *rmap, *cmap; 150 { 151 int prow, k, nz, n, repl, *j, *col, *row, cloc; 152 double *v, repla; 153 double *diagv, offd, rval; 154 155 SPLITTOMAT(mat); 156 157 n = mat->rows; 158 col = cmap; 159 row = rmap; 160 diagv = (double *)MALLOC( n * sizeof(double) ); CHKPTR(diagv); 161 /* load the diagonal values */ 162 for (prow=0; prow<n; prow++) { 163 SpScatterFromRow( mat, row[prow], &nz, &j, &v ); 164 rval = 0.0; 165 for (k=0; k<nz; k++) 166 if (col[j[k]] == prow) { 167 rval = v[k]; 168 break; 169 } 170 if (k >= nz) { 171 /* No diagonal element. Add it */ 172 SpAddValue( mat, 0.0, prow, prow ); 173 } 174 diagv[row[prow]] = fabs(rval); 175 } 176 177 for (prow=0; prow<n; prow++) { 178 if (diagv[prow] > atol) continue; 179 /* If we knew that the column indices in a row were sorted, we'd only 180 have to look at the values after the diagonal's location. Since 181 we are not making that assumption, we have to look at every element. 182 This lets us use this routine after applying a fill-reducing ordering 183 (though in that case, the diagonal elements should already have been 184 introduced into the structure). */ 185 /* For each column, look for the largest diagonal element */ 186 SpScatterFromRow( mat, row[prow], &nz, &j, &v ); 187 repl = prow; 188 repla = diagv[row[prow]]; 189 for (k=0; k<nz; k++) 190 if ((cloc = col[j[k]]) > prow && diagv[cloc] > repla) { 191 repl = j[k]; 192 repla = diagv[cloc]; 193 offd = v[k]; 194 } 195 if (prow == repl) { 196 SETERRC(1,"Can not reorder matrix"); 197 return; 198 } 199 /* Compute an estimate of the replaced diagonal */ 200 diagv[row[prow]] = fabs( offd * offd / repla ); 201 SWAPD(diagv[row[prow]],diagv[row[repl]]); 202 203 SWAP(col[prow],col[repl]); 204 if (col != row) 205 SWAP(row[prow],row[repl]); 206 } 207 FREE(diagv); 208 } 209 210 /* Given a current row and current permutation, find a column permutation 211 that removes a zero diagonal */ 212 int SpiZeroFindPre( mat, prow, row, col, repla, atol, rc, rcv ) 213 SpMat *mat; 214 int prow, *row, *col, *rc; 215 double repla, atol, *rcv; 216 { 217 int k, nz, repl, *j, kk, nnz, *jj; 218 double *v, *vv; 219 220 SpScatterFromRow( mat, row[prow], &nz, &j, &v ); 221 for (k=0; k<nz; k++) { 222 if (col[j[k]] < prow && fabs(v[k]) > repla) { 223 /* See if this one will work */ 224 repl = col[j[k]]; 225 SpScatterFromRow( mat, row[repl], &nnz, &jj, &vv ); 226 for (kk=0; kk<nnz; kk++) { 227 if (col[jj[kk]] == prow && fabs(vv[kk]) > atol) { 228 *rcv = fabs(v[k]); 229 *rc = repl; 230 return 1; 231 } 232 } 233 } 234 } 235 return 0; 236 } 237