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