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 @*/
MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis)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
MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis)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