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>       /*I  "petscmat.h"  I*/

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

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

 18:     Collective on Mat

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

 25:     Level: intermediate

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

 33:     Only works for SeqAIJ matrices

 35:     Options Database Keys (When using KSP):
 36: .      -pc_factor_nonzeros_along_diagonal

 38:     Algorithm Notes:
 39:     Column pivoting is used. 

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


 52: @*/
 53: PetscErrorCode  MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis)
 54: {

 58:   PetscTryMethod(mat,"MatReorderForNonzeroDiagonal_C",(Mat,PetscReal,IS,IS),(mat,abstol,ris,cis));
 59:   return(0);
 60: }

 62: extern PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);
 63: extern PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);

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

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

 79:   /* access the indices of the IS directly, because it changes them */
 80:   row  = ((IS_General*)ris->data)->idx;
 81:   col  = ((IS_General*)cis->data)->idx;
 82:   ISInvertPermutation(cis,PETSC_DECIDE,&icis);
 83:   icol  = ((IS_General*)icis->data)->idx;
 84:   MatGetSize(mat,&m,&n);

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

141:       found:;
142:     }
143:     MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
144:   }
145:   ISDestroy(&icis);
146:   return(0);
147: }
148: EXTERN_C_END