1 /*$Id: zerodiag.c,v 1.44 2001/08/06 21:16:10 bsmith Exp $*/ 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 __FUNCT__ 13 #define __FUNCT__ "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 Only works for SeqAIJ matrices 35 36 Options Database Keys (When using SLES): 37 + -pc_ilu_nonzeros_along_diagonal 38 - -pc_lu_nonzeros_along_diagonal 39 40 Algorithm Notes: 41 Column pivoting is used. 42 43 1) Choice of column is made by looking at the 44 non-zero elements in the troublesome row for columns that are not yet 45 included (moving from left to right). 46 47 2) If (1) fails we check all the columns to the left of the current row 48 and see if one of them has could be swapped. It can be swapped if 49 its corresponding row has a non-zero in the column it is being 50 swapped with; to make sure the previous nonzero diagonal remains 51 nonzero 52 53 54 @*/ 55 int MatReorderForNonzeroDiagonal(Mat mat,PetscReal atol,IS ris,IS cis) 56 { 57 int ierr,(*f)(Mat,PetscReal,IS,IS); 58 59 PetscFunctionBegin; 60 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatReorderForNonzeroDiagonal_C",(void (**)(void))&f);CHKERRQ(ierr); 61 if (f) { 62 ierr = (*f)(mat,atol,ris,cis);CHKERRQ(ierr); 63 } 64 PetscFunctionReturn(0); 65 } 66 67 EXTERN_C_BEGIN 68 #undef __FUNCT__ 69 #define __FUNCT__ "MatReorderForNonzeroDiagonal_SeqAIJ" 70 int MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal atol,IS ris,IS cis) 71 { 72 int ierr,prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk; 73 PetscScalar *v,*vv; 74 PetscReal repla; 75 IS icis; 76 PetscTruth flg; 77 78 PetscFunctionBegin; 79 PetscValidHeaderSpecific(mat,MAT_COOKIE); 80 PetscValidHeaderSpecific(ris,IS_COOKIE); 81 PetscValidHeaderSpecific(cis,IS_COOKIE); 82 83 ierr = ISGetIndices(ris,&row);CHKERRQ(ierr); 84 ierr = ISGetIndices(cis,&col);CHKERRQ(ierr); 85 ierr = ISInvertPermutation(cis,PETSC_DECIDE,&icis);CHKERRQ(ierr); 86 ierr = ISGetIndices(icis,&icol);CHKERRQ(ierr); 87 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 88 89 for (prow=0; prow<n; prow++) { 90 ierr = MatGetRow(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 91 for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;} 92 if (k >= nz || PetscAbsScalar(v[k]) <= atol) { 93 /* Element too small or zero; find the best candidate */ 94 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 95 /* 96 Look for a later column we can swap with this one 97 */ 98 for (k=0; k<nz; k++) { 99 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 100 /* found a suitable later column */ 101 repl = icol[j[k]]; 102 SWAP(icol[col[prow]],icol[col[repl]]); 103 SWAP(col[prow],col[repl]); 104 goto found; 105 } 106 } 107 /* 108 Did not find a suitable later column so look for an earlier column 109 We need to be sure that we don't introduce a zero in a previous 110 diagonal 111 */ 112 for (k=0; k<nz; k++) { 113 if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 114 /* See if this one will work */ 115 repl = icol[j[k]]; 116 ierr = MatGetRow(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 117 for (kk=0; kk<nnz; kk++) { 118 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 119 ierr = MatRestoreRow(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 120 SWAP(icol[col[prow]],icol[col[repl]]); 121 SWAP(col[prow],col[repl]); 122 goto found; 123 } 124 } 125 ierr = MatRestoreRow(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 126 } 127 } 128 /* 129 No column suitable; instead check all future rows 130 Note: this will be very slow 131 */ 132 for (k=prow+1; k<n; k++) { 133 ierr = MatGetRow(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 134 for (kk=0; kk<nnz; kk++) { 135 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 136 /* found a row */ 137 SWAP(row[prow],row[k]); 138 goto found; 139 } 140 } 141 ierr = MatRestoreRow(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 142 } 143 144 found:; 145 } 146 ierr = MatRestoreRow(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 147 } 148 ierr = ISRestoreIndices(ris,&row);CHKERRQ(ierr); 149 ierr = ISRestoreIndices(cis,&col);CHKERRQ(ierr); 150 ierr = ISRestoreIndices(icis,&icol);CHKERRQ(ierr); 151 ierr = ISDestroy(icis);CHKERRQ(ierr); 152 PetscFunctionReturn(0); 153 } 154 EXTERN_C_END 155 156 157