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