1 2 /* 3 This file contains routines to reorder a matrix so that the diagonal 4 elements are nonzero. 5 */ 6 7 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 8 9 #define SWAP(a, b) \ 10 { \ 11 PetscInt _t; \ 12 _t = a; \ 13 a = b; \ 14 b = _t; \ 15 } 16 17 /*@ 18 MatReorderForNonzeroDiagonal - Changes matrix ordering to remove 19 zeros from diagonal. This may help in the `PCLU` factorization to 20 prevent a zero pivot. 21 22 Collective 23 24 Input Parameters: 25 + mat - matrix to reorder 26 . abstol - absolute tolerance, it attempts to move all values smaller off the diagonal 27 . ris - the row reordering 28 - cis - the column reordering; this may be changed 29 30 Level: intermediate 31 32 Options Database Key: 33 . -pc_factor_nonzeros_along_diagonal - Reorder to remove zeros from diagonal 34 35 Notes: 36 This is not intended as a replacement for pivoting for matrices that 37 have ``bad'' structure. It is only a stop-gap measure. 38 39 Should be called 40 after a call to `MatGetOrdering()`. 41 42 Only works for `MATSEQAIJ` matrices 43 44 Developer Notes: 45 Column pivoting is used. 46 47 1) Choice of column is made by looking at the 48 non-zero elements in the troublesome row for columns that are not yet 49 included (moving from left to right). 50 51 2) If (1) fails we check all the columns to the left of the current row 52 and see if one of them has could be swapped. It can be swapped if 53 its corresponding row has a non-zero in the column it is being 54 swapped with; to make sure the previous nonzero diagonal remains 55 nonzero 56 57 .seealso: `Mat`, `MatGetFactor()`, `MatGetOrdering()` 58 @*/ 59 PetscErrorCode MatReorderForNonzeroDiagonal(Mat mat, PetscReal abstol, IS ris, IS cis) 60 { 61 PetscFunctionBegin; 62 PetscTryMethod(mat, "MatReorderForNonzeroDiagonal_C", (Mat, PetscReal, IS, IS), (mat, abstol, ris, cis)); 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 66 PETSC_INTERN PetscErrorCode MatGetRow_SeqAIJ(Mat, PetscInt, PetscInt *, PetscInt **, PetscScalar **); 67 PETSC_INTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat, PetscInt, PetscInt *, PetscInt **, PetscScalar **); 68 69 #include <../src/vec/is/is/impls/general/general.h> 70 71 PETSC_INTERN PetscErrorCode MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat, PetscReal abstol, IS ris, IS cis) 72 { 73 PetscInt prow, k, nz, n, repl, *j, *col, *row, m, *icol, nnz, *jj, kk; 74 PetscScalar *v, *vv; 75 PetscReal repla; 76 IS icis; 77 78 PetscFunctionBegin; 79 /* access the indices of the IS directly, because it changes them */ 80 row = ((IS_General *)ris->data)->idx; 81 col = ((IS_General *)cis->data)->idx; 82 PetscCall(ISInvertPermutation(cis, PETSC_DECIDE, &icis)); 83 icol = ((IS_General *)icis->data)->idx; 84 PetscCall(MatGetSize(mat, &m, &n)); 85 86 for (prow = 0; prow < n; prow++) { 87 PetscCall(MatGetRow_SeqAIJ(mat, row[prow], &nz, &j, &v)); 88 for (k = 0; k < nz; k++) { 89 if (icol[j[k]] == prow) break; 90 } 91 if (k >= nz || PetscAbsScalar(v[k]) <= abstol) { 92 /* Element too small or zero; find the best candidate */ 93 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 94 /* 95 Look for a later column we can swap with this one 96 */ 97 for (k = 0; k < nz; k++) { 98 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 99 /* found a suitable later column */ 100 repl = icol[j[k]]; 101 SWAP(icol[col[prow]], icol[col[repl]]); 102 SWAP(col[prow], col[repl]); 103 goto found; 104 } 105 } 106 /* 107 Did not find a suitable later column so look for an earlier column 108 We need to be sure that we don't introduce a zero in a previous 109 diagonal 110 */ 111 for (k = 0; k < nz; k++) { 112 if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 113 /* See if this one will work */ 114 repl = icol[j[k]]; 115 PetscCall(MatGetRow_SeqAIJ(mat, row[repl], &nnz, &jj, &vv)); 116 for (kk = 0; kk < nnz; kk++) { 117 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) { 118 PetscCall(MatRestoreRow_SeqAIJ(mat, row[repl], &nnz, &jj, &vv)); 119 SWAP(icol[col[prow]], icol[col[repl]]); 120 SWAP(col[prow], col[repl]); 121 goto found; 122 } 123 } 124 PetscCall(MatRestoreRow_SeqAIJ(mat, row[repl], &nnz, &jj, &vv)); 125 } 126 } 127 /* 128 No column suitable; instead check all future rows 129 Note: this will be very slow 130 */ 131 for (k = prow + 1; k < n; k++) { 132 PetscCall(MatGetRow_SeqAIJ(mat, row[k], &nnz, &jj, &vv)); 133 for (kk = 0; kk < nnz; kk++) { 134 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) { 135 /* found a row */ 136 SWAP(row[prow], row[k]); 137 goto found; 138 } 139 } 140 PetscCall(MatRestoreRow_SeqAIJ(mat, row[k], &nnz, &jj, &vv)); 141 } 142 143 found:; 144 } 145 PetscCall(MatRestoreRow_SeqAIJ(mat, row[prow], &nz, &j, &v)); 146 } 147 PetscCall(ISDestroy(&icis)); 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150