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