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