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
50: @*/
51: PetscErrorCodeMatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis) 52: {
56: PetscTryMethod(mat,"MatReorderForNonzeroDiagonal_C",(Mat,PetscReal,IS,IS),(mat,abstol,ris,cis));
57: return(0);
58: }
60: PETSC_INTERN PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);
61: PETSC_INTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);
63: #include <../src/vec/is/is/impls/general/general.h> 65: PETSC_INTERN PetscErrorCode MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis) 66: {
68: PetscInt prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk;
69: PetscScalar *v,*vv;
70: PetscReal repla;
71: IS icis;
74: /* access the indices of the IS directly, because it changes them */
75: row = ((IS_General*)ris->data)->idx;
76: col = ((IS_General*)cis->data)->idx;
77: ISInvertPermutation(cis,PETSC_DECIDE,&icis);
78: icol = ((IS_General*)icis->data)->idx;
79: MatGetSize(mat,&m,&n);
81: for (prow=0; prow<n; prow++) {
82: MatGetRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
83: for (k=0; k<nz; k++) {
84: if (icol[j[k]] == prow) break;
85: }
86: if (k >= nz || PetscAbsScalar(v[k]) <= abstol) {
87: /* Element too small or zero; find the best candidate */
88: repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]);
89: /*
90: Look for a later column we can swap with this one
91: */
92: for (k=0; k<nz; k++) {
93: if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) {
94: /* found a suitable later column */
95: repl = icol[j[k]];
96: SWAP(icol[col[prow]],icol[col[repl]]);
97: SWAP(col[prow],col[repl]);
98: goto found;
99: }
100: }
101: /*
102: Did not find a suitable later column so look for an earlier column
103: We need to be sure that we don't introduce a zero in a previous
104: diagonal
105: */
106: for (k=0; k<nz; k++) {
107: if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) {
108: /* See if this one will work */
109: repl = icol[j[k]];
110: MatGetRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
111: for (kk=0; kk<nnz; kk++) {
112: if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
113: MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
114: SWAP(icol[col[prow]],icol[col[repl]]);
115: SWAP(col[prow],col[repl]);
116: goto found;
117: }
118: }
119: MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
120: }
121: }
122: /*
123: No column suitable; instead check all future rows
124: Note: this will be very slow
125: */
126: for (k=prow+1; k<n; k++) {
127: MatGetRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);
128: for (kk=0; kk<nnz; kk++) {
129: if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
130: /* found a row */
131: SWAP(row[prow],row[k]);
132: goto found;
133: }
134: }
135: MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);
136: }
138: found:;
139: }
140: MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
141: }
142: ISDestroy(&icis);
143: return(0);
144: }