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 include/private/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_factor_nonzeros_along_diagonal
39: Algorithm Notes:
40: Column pivoting is used.
42: 1) Choice of column is made by looking at the
43: non-zero elements in the troublesome row for columns that are not yet
44: included (moving from left to right).
45:
46: 2) If (1) fails we check all the columns to the left of the current row
47: and see if one of them has could be swapped. It can be swapped if
48: its corresponding row has a non-zero in the column it is being
49: swapped with; to make sure the previous nonzero diagonal remains
50: nonzero
53: @*/
54: PetscErrorCode MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis)
55: {
56: PetscErrorCode ierr,(*f)(Mat,PetscReal,IS,IS);
59: PetscObjectQueryFunction((PetscObject)mat,"MatReorderForNonzeroDiagonal_C",(void (**)(void))&f);
60: if (f) {
61: (*f)(mat,abstol,ris,cis);
62: }
63: return(0);
64: }
66: EXTERN PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);
67: EXTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**);
72: PetscErrorCode MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis)
73: {
75: PetscInt prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk;
76: PetscScalar *v,*vv;
77: PetscReal repla;
78: IS icis;
81: ISGetIndices(ris,&row);
82: ISGetIndices(cis,&col);
83: ISInvertPermutation(cis,PETSC_DECIDE,&icis);
84: ISGetIndices(icis,&icol);
85: MatGetSize(mat,&m,&n);
87: for (prow=0; prow<n; prow++) {
88: MatGetRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
89: for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;}
90: if (k >= nz || PetscAbsScalar(v[k]) <= abstol) {
91: /* Element too small or zero; find the best candidate */
92: repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]);
93: /*
94: Look for a later column we can swap with this one
95: */
96: for (k=0; k<nz; k++) {
97: if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) {
98: /* found a suitable later column */
99: repl = icol[j[k]];
100: SWAP(icol[col[prow]],icol[col[repl]]);
101: SWAP(col[prow],col[repl]);
102: goto found;
103: }
104: }
105: /*
106: Did not find a suitable later column so look for an earlier column
107: We need to be sure that we don't introduce a zero in a previous
108: diagonal
109: */
110: for (k=0; k<nz; k++) {
111: if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) {
112: /* See if this one will work */
113: repl = icol[j[k]];
114: MatGetRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
115: for (kk=0; kk<nnz; kk++) {
116: if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
117: MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
118: SWAP(icol[col[prow]],icol[col[repl]]);
119: SWAP(col[prow],col[repl]);
120: goto found;
121: }
122: }
123: MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);
124: }
125: }
126: /*
127: No column suitable; instead check all future rows
128: Note: this will be very slow
129: */
130: for (k=prow+1; k<n; k++) {
131: MatGetRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);
132: for (kk=0; kk<nnz; kk++) {
133: if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) {
134: /* found a row */
135: SWAP(row[prow],row[k]);
136: goto found;
137: }
138: }
139: MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);
140: }
142: found:;
143: }
144: MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);
145: }
146: ISRestoreIndices(ris,&row);
147: ISRestoreIndices(cis,&col);
148: ISRestoreIndices(icis,&icol);
149: ISDestroy(icis);
150: return(0);
151: }