1 /*$Id: zerodiag.c,v 1.38 2000/04/12 04:24:13 bsmith Exp balay $*/ 2 3 /* 4 This file contains routines to reorder a matrix so that the diagonal 5 elements are nonzero. 6 */ 7 8 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 9 10 #define SWAP(a,b) {int _t; _t = a; a = b; b = _t; } 11 12 #undef __FUNC__ 13 #define __FUNC__ /*<a name=""></a>*/"MatReorderForNonzeroDiagonal" 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. 18 19 Collective on Mat 20 21 Input Parameters: 22 + mat - matrix to reorder 23 - rmap,cmap - row and column permutations. Usually obtained from 24 MatGetOrdering(). 25 26 Level: intermediate 27 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. 33 34 Options Database Keys (When using SLES): 35 + -pc_ilu_nonzeros_along_diagonal 36 - -pc_lu_nonzeros_along_diagonal 37 38 Algorithm Notes: 39 Column pivoting is used. 40 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 50 51 52 @*/ 53 int MatReorderForNonzeroDiagonal(Mat mat,PetscReal atol,IS ris,IS cis) 54 { 55 int ierr,prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk; 56 Scalar *v,*vv; 57 PetscReal repla; 58 IS icis; 59 60 PetscFunctionBegin; 61 PetscValidHeaderSpecific(mat,MAT_COOKIE); 62 PetscValidHeaderSpecific(ris,IS_COOKIE); 63 PetscValidHeaderSpecific(cis,IS_COOKIE); 64 65 ierr = ISGetIndices(ris,&row);CHKERRQ(ierr); 66 ierr = ISGetIndices(cis,&col);CHKERRQ(ierr); 67 ierr = ISInvertPermutation(cis,PETSC_DECIDE,&icis);CHKERRQ(ierr); 68 ierr = ISGetIndices(icis,&icol);CHKERRQ(ierr); 69 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 70 71 for (prow=0; prow<n; prow++) { 72 ierr = MatGetRow(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 73 for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;} 74 if (k >= nz || PetscAbsScalar(v[k]) <= atol) { 75 /* Element too small or zero; find the best candidate */ 76 repl = prow; 77 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 78 /* 79 Look for a later column we can swap with this one 80 */ 81 for (k=0; k<nz; k++) { 82 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 83 /* found a suitable later column */ 84 repl = icol[j[k]]; 85 repla = PetscAbsScalar(v[k]); 86 SWAP(icol[col[prow]],icol[col[repl]]); 87 SWAP(col[prow],col[repl]); 88 goto found; 89 } 90 } 91 /* 92 Did not find a suitable later column so look for an earlier column 93 We need to be sure that we don't introduce a zero in a previous 94 diagonal 95 */ 96 for (k=0; k<nz; k++) { 97 if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 98 /* See if this one will work */ 99 repl = icol[j[k]]; 100 ierr = MatGetRow(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 101 for (kk=0; kk<nnz; kk++) { 102 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 103 repla = PetscAbsScalar(v[k]); 104 ierr = MatRestoreRow(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 105 SWAP(icol[col[prow]],icol[col[repl]]); 106 SWAP(col[prow],col[repl]); 107 goto found; 108 } 109 } 110 ierr = MatRestoreRow(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 111 } 112 } 113 /* 114 No column suitable; instead check all future rows 115 Note: this will be very slow 116 */ 117 for (k=prow+1; k<n; k++) { 118 ierr = MatGetRow(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 119 for (kk=0; kk<nnz; kk++) { 120 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 121 /* found a row */ 122 SWAP(row[prow],row[k]); 123 goto found; 124 } 125 } 126 ierr = MatRestoreRow(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 127 } 128 129 found:; 130 } 131 ierr = MatRestoreRow(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 132 } 133 ierr = ISRestoreIndices(ris,&row);CHKERRQ(ierr); 134 ierr = ISRestoreIndices(cis,&col);CHKERRQ(ierr); 135 ierr = ISRestoreIndices(icis,&icol);CHKERRQ(ierr); 136 ierr = ISDestroy(icis);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 141 142