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