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: - rmap,cmap - row and column permutations. Usually obtained from
27: `MatGetOrdering()`.
29: Level: intermediate
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.
37: Only works for `MATSEQAIJ` matrices
39: Options Database Keys (When using `KSP`):
40: . -pc_factor_nonzeros_along_diagonal - Reorder to remove zeros from diagonal
42: Algorithm Notes:
43: Column pivoting is used.
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).
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
55: @*/
56: PetscErrorCode MatReorderForNonzeroDiagonal(Mat mat, PetscReal abstol, IS ris, IS cis)
57: {
58: PetscTryMethod(mat, "MatReorderForNonzeroDiagonal_C", (Mat, PetscReal, IS, IS), (mat, abstol, ris, cis));
59: return 0;
60: }
62: PETSC_INTERN PetscErrorCode MatGetRow_SeqAIJ(Mat, PetscInt, PetscInt *, PetscInt **, PetscScalar **);
63: PETSC_INTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat, PetscInt, PetscInt *, PetscInt **, PetscScalar **);
65: #include <../src/vec/is/is/impls/general/general.h>
67: PETSC_INTERN PetscErrorCode MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat, PetscReal abstol, IS ris, IS cis)
68: {
69: PetscInt prow, k, nz, n, repl, *j, *col, *row, m, *icol, nnz, *jj, kk;
70: PetscScalar *v, *vv;
71: PetscReal repla;
72: IS icis;
74: /* access the indices of the IS directly, because it changes them */
75: row = ((IS_General *)ris->data)->idx;
76: col = ((IS_General *)cis->data)->idx;
77: ISInvertPermutation(cis, PETSC_DECIDE, &icis);
78: icol = ((IS_General *)icis->data)->idx;
79: MatGetSize(mat, &m, &n);
81: for (prow = 0; prow < n; prow++) {
82: MatGetRow_SeqAIJ(mat, row[prow], &nz, &j, &v);
83: for (k = 0; k < nz; k++) {
84: if (icol[j[k]] == prow) break;
85: }
86: if (k >= nz || PetscAbsScalar(v[k]) <= abstol) {
87: /* Element too small or zero; find the best candidate */
88: repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]);
89: /*
90: Look for a later column we can swap with this one
91: */
92: for (k = 0; k < nz; k++) {
93: if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) {
94: /* found a suitable later column */
95: repl = icol[j[k]];
96: SWAP(icol[col[prow]], icol[col[repl]]);
97: SWAP(col[prow], col[repl]);
98: goto found;
99: }
100: }
101: /*
102: Did not find a suitable later column so look for an earlier column
103: We need to be sure that we don't introduce a zero in a previous
104: diagonal
105: */
106: for (k = 0; k < nz; k++) {
107: if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) {
108: /* See if this one will work */
109: repl = icol[j[k]];
110: MatGetRow_SeqAIJ(mat, row[repl], &nnz, &jj, &vv);
111: for (kk = 0; kk < nnz; kk++) {
112: if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
113: MatRestoreRow_SeqAIJ(mat, row[repl], &nnz, &jj, &vv);
114: SWAP(icol[col[prow]], icol[col[repl]]);
115: SWAP(col[prow], col[repl]);
116: goto found;
117: }
118: }
119: MatRestoreRow_SeqAIJ(mat, row[repl], &nnz, &jj, &vv);
120: }
121: }
122: /*
123: No column suitable; instead check all future rows
124: Note: this will be very slow
125: */
126: for (k = prow + 1; k < n; k++) {
127: MatGetRow_SeqAIJ(mat, row[k], &nnz, &jj, &vv);
128: for (kk = 0; kk < nnz; kk++) {
129: if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
130: /* found a row */
131: SWAP(row[prow], row[k]);
132: goto found;
133: }
134: }
135: MatRestoreRow_SeqAIJ(mat, row[k], &nnz, &jj, &vv);
136: }
138: found:;
139: }
140: MatRestoreRow_SeqAIJ(mat, row[prow], &nz, &j, &v);
141: }
142: ISDestroy(&icis);
143: return 0;
144: }