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