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) {PetscInt _t; _t = a; a = b; b = _t; }

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

 16:     Collective on Mat

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

 23:     Level: intermediate

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

 31:     Only works for SeqAIJ matrices

 33:     Options Database Keys (When using KSP):
 34: .      -pc_factor_nonzeros_along_diagonal

 36:     Algorithm Notes:
 37:     Column pivoting is used.

 39:     1) Choice of column is made by looking at the
 40:        non-zero elements in the troublesome row for columns that are not yet
 41:        included (moving from left to right).

 43:     2) If (1) fails we check all the columns to the left of the current row
 44:        and see if one of them has could be swapped. It can be swapped if
 45:        its corresponding row has a non-zero in the column it is being
 46:        swapped with; to make sure the previous nonzero diagonal remains
 47:        nonzero

 49: @*/
 50: PetscErrorCode  MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis)
 51: {

 55:   PetscTryMethod(mat,"MatReorderForNonzeroDiagonal_C",(Mat,PetscReal,IS,IS),(mat,abstol,ris,cis));
 56:   return(0);
 57: }

 59: PETSC_INTERN PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);
 60: PETSC_INTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);

 62: #include <../src/vec/is/is/impls/general/general.h>

 64: PETSC_INTERN PetscErrorCode  MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis)
 65: {
 67:   PetscInt       prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk;
 68:   PetscScalar    *v,*vv;
 69:   PetscReal      repla;
 70:   IS             icis;

 73:   /* access the indices of the IS directly, because it changes them */
 74:   row  = ((IS_General*)ris->data)->idx;
 75:   col  = ((IS_General*)cis->data)->idx;
 76:   ISInvertPermutation(cis,PETSC_DECIDE,&icis);
 77:   icol = ((IS_General*)icis->data)->idx;
 78:   MatGetSize(mat,&m,&n);

 80:   for (prow=0; prow<n; prow++) {
 81:     MatGetRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
 82:     for (k=0; k<nz; k++) {
 83:       if (icol[j[k]] == prow) break;
 84:     }
 85:     if (k >= nz || PetscAbsScalar(v[k]) <= abstol) {
 86:       /* Element too small or zero; find the best candidate */
 87:       repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]);
 88:       /*
 89:           Look for a later column we can swap with this one
 90:       */
 91:       for (k=0; k<nz; k++) {
 92:         if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) {
 93:           /* found a suitable later column */
 94:           repl = icol[j[k]];
 95:           SWAP(icol[col[prow]],icol[col[repl]]);
 96:           SWAP(col[prow],col[repl]);
 97:           goto found;
 98:         }
 99:       }
100:       /*
101:            Did not find a suitable later column so look for an earlier column
102:            We need to be sure that we don't introduce a zero in a previous
103:            diagonal
104:       */
105:       for (k=0; k<nz; k++) {
106:         if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) {
107:           /* See if this one will work */
108:           repl = icol[j[k]];
109:           MatGetRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
110:           for (kk=0; kk<nnz; kk++) {
111:             if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
112:               MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
113:               SWAP(icol[col[prow]],icol[col[repl]]);
114:               SWAP(col[prow],col[repl]);
115:               goto found;
116:             }
117:           }
118:           MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
119:         }
120:       }
121:       /*
122:           No column  suitable; instead check all future rows
123:           Note: this will be very slow
124:       */
125:       for (k=prow+1; k<n; k++) {
126:         MatGetRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);
127:         for (kk=0; kk<nnz; kk++) {
128:           if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
129:             /* found a row */
130:             SWAP(row[prow],row[k]);
131:             goto found;
132:           }
133:         }
134:         MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);
135:       }

137: found:;
138:     }
139:     MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
140:   }
141:   ISDestroy(&icis);
142:   return(0);
143: }