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: }