Actual source code: zerodiag.c
2: /*
3: This file contains routines to reorder a matrix so that the diagonal
4: elements are nonzero.
5: */
7: #include <petsc/private/matimpl.h>
9: #define SWAP(a, b) \
10: { \
11: PetscInt _t; \
12: _t = a; \
13: a = b; \
14: b = _t; \
15: }
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.
22: Collective
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
30: Level: intermediate
32: Options Database Key:
33: . -pc_factor_nonzeros_along_diagonal - Reorder to remove zeros from diagonal
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.
39: Should be called
40: after a call to `MatGetOrdering()`.
42: Only works for `MATSEQAIJ` matrices
44: Developer Notes:
45: Column pivoting is used.
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).
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
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: }
66: PETSC_INTERN PetscErrorCode MatGetRow_SeqAIJ(Mat, PetscInt, PetscInt *, PetscInt **, PetscScalar **);
67: PETSC_INTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat, PetscInt, PetscInt *, PetscInt **, PetscScalar **);
69: #include <../src/vec/is/is/impls/general/general.h>
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;
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));
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: }
143: found:;
144: }
145: PetscCall(MatRestoreRow_SeqAIJ(mat, row[prow], &nz, &j, &v));
146: }
147: PetscCall(ISDestroy(&icis));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }