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