Actual source code: zerodiag.c

      1: #define PETSCMAT_DLL

      3: /*
      4:     This file contains routines to reorder a matrix so that the diagonal
      5:     elements are nonzero.
      6:  */

     8:  #include src/mat/matimpl.h

     10: #define SWAP(a,b) {PetscInt _t; _t = a; a = b; b = _t; }

     14: /*@
     15:     MatReorderForNonzeroDiagonal - Changes matrix ordering to remove
     16:     zeros from diagonal. This may help in the LU factorization to
     17:     prevent a zero pivot.

     19:     Collective on Mat

     21:     Input Parameters:
     22: +   mat  - matrix to reorder
     23: -   rmap,cmap - row and column permutations.  Usually obtained from
     24:                MatGetOrdering().

     26:     Level: intermediate

     28:     Notes:
     29:     This is not intended as a replacement for pivoting for matrices that
     30:     have ``bad'' structure. It is only a stop-gap measure. Should be called
     31:     after a call to MatGetOrdering(), this routine changes the column
     32:     ordering defined in cis.

     34:     Only works for SeqAIJ matrices

     36:     Options Database Keys (When using KSP):
     37: +      -pc_ilu_nonzeros_along_diagonal
     38: -      -pc_lu_nonzeros_along_diagonal

     40:     Algorithm Notes:
     41:     Column pivoting is used.

     43:     1) Choice of column is made by looking at the
     44:        non-zero elements in the troublesome row for columns that are not yet
     45:        included (moving from left to right).
     46:
     47:     2) If (1) fails we check all the columns to the left of the current row
     48:        and see if one of them has could be swapped. It can be swapped if
     49:        its corresponding row has a non-zero in the column it is being
     50:        swapped with; to make sure the previous nonzero diagonal remains
     51:        nonzero

     54: @*/
     55: PetscErrorCode PETSCMAT_DLLEXPORT MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis)
     56: {
     57:   PetscErrorCode ierr,(*f)(Mat,PetscReal,IS,IS);

     60:   PetscObjectQueryFunction((PetscObject)mat,"MatReorderForNonzeroDiagonal_C",(void (**)(void))&f);
     61:   if (f) {
     62:     (*f)(mat,abstol,ris,cis);
     63:   }
     64:   return(0);
     65: }

     67: EXTERN PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);
     68: EXTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);

     73: PetscErrorCode PETSCMAT_DLLEXPORT MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis)
     74: {
     76:   PetscInt       prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk;
     77:   PetscScalar    *v,*vv;
     78:   PetscReal      repla;
     79:   IS             icis;

     82:   ISGetIndices(ris,&row);
     83:   ISGetIndices(cis,&col);
     84:   ISInvertPermutation(cis,PETSC_DECIDE,&icis);
     85:   ISGetIndices(icis,&icol);
     86:   MatGetSize(mat,&m,&n);

     88:   for (prow=0; prow<n; prow++) {
     89:     MatGetRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
     90:     for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;}
     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:           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:               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:           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:         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:         MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);
    141:       }

    143:       found:;
    144:     }
    145:     MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
    146:   }
    147:   ISRestoreIndices(ris,&row);
    148:   ISRestoreIndices(cis,&col);
    149:   ISRestoreIndices(icis,&icol);
    150:   ISDestroy(icis);
    151:   return(0);
    152: }