1 /* 2 This file contains routines to reorder a matrix so that the diagonal 3 elements are nonzero. 4 */ 5 6 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 7 8 #define SWAP(a,b) {int _t; _t = a; a = b; b = _t; } 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "MatReorderForNonzeroDiagonal" 12 /*@ 13 MatReorderForNonzeroDiagonal - Changes matrix ordering to remove 14 zeros from diagonal. This may help in the LU factorization to 15 prevent a zero pivot. 16 17 Collective on Mat 18 19 Input Parameters: 20 + mat - matrix to reorder 21 - rmap,cmap - row and column permutations. Usually obtained from 22 MatGetOrdering(). 23 24 Level: intermediate 25 26 Notes: 27 This is not intended as a replacement for pivoting for matrices that 28 have ``bad'' structure. It is only a stop-gap measure. Should be called 29 after a call to MatGetOrdering(), this routine changes the column 30 ordering defined in cis. 31 32 Only works for SeqAIJ matrices 33 34 Options Database Keys (When using KSP): 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,(*f)(Mat,PetscReal,IS,IS); 56 57 PetscFunctionBegin; 58 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatReorderForNonzeroDiagonal_C",(void (**)(void))&f);CHKERRQ(ierr); 59 if (f) { 60 ierr = (*f)(mat,atol,ris,cis);CHKERRQ(ierr); 61 } 62 PetscFunctionReturn(0); 63 } 64 65 EXTERN int MatGetRow_SeqAIJ(Mat,int,int*,int**,PetscScalar**); 66 EXTERN int MatRestoreRow_SeqAIJ(Mat,int,int*,int**,PetscScalar**); 67 68 EXTERN_C_BEGIN 69 #undef __FUNCT__ 70 #define __FUNCT__ "MatReorderForNonzeroDiagonal_SeqAIJ" 71 int MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal atol,IS ris,IS cis) 72 { 73 int ierr,prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk; 74 PetscScalar *v,*vv; 75 PetscReal repla; 76 IS icis; 77 78 PetscFunctionBegin; 79 ierr = ISGetIndices(ris,&row);CHKERRQ(ierr); 80 ierr = ISGetIndices(cis,&col);CHKERRQ(ierr); 81 ierr = ISInvertPermutation(cis,PETSC_DECIDE,&icis);CHKERRQ(ierr); 82 ierr = ISGetIndices(icis,&icol);CHKERRQ(ierr); 83 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 84 85 for (prow=0; prow<n; prow++) { 86 ierr = MatGetRow_SeqAIJ(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 87 for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;} 88 if (k >= nz || PetscAbsScalar(v[k]) <= atol) { 89 /* Element too small or zero; find the best candidate */ 90 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 91 /* 92 Look for a later column we can swap with this one 93 */ 94 for (k=0; k<nz; k++) { 95 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 96 /* found a suitable later column */ 97 repl = icol[j[k]]; 98 SWAP(icol[col[prow]],icol[col[repl]]); 99 SWAP(col[prow],col[repl]); 100 goto found; 101 } 102 } 103 /* 104 Did not find a suitable later column so look for an earlier column 105 We need to be sure that we don't introduce a zero in a previous 106 diagonal 107 */ 108 for (k=0; k<nz; k++) { 109 if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 110 /* See if this one will work */ 111 repl = icol[j[k]]; 112 ierr = MatGetRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 113 for (kk=0; kk<nnz; kk++) { 114 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 115 ierr = MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 116 SWAP(icol[col[prow]],icol[col[repl]]); 117 SWAP(col[prow],col[repl]); 118 goto found; 119 } 120 } 121 ierr = MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 122 } 123 } 124 /* 125 No column suitable; instead check all future rows 126 Note: this will be very slow 127 */ 128 for (k=prow+1; k<n; k++) { 129 ierr = MatGetRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 130 for (kk=0; kk<nnz; kk++) { 131 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 132 /* found a row */ 133 SWAP(row[prow],row[k]); 134 goto found; 135 } 136 } 137 ierr = MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 138 } 139 140 found:; 141 } 142 ierr = MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 143 } 144 ierr = ISRestoreIndices(ris,&row);CHKERRQ(ierr); 145 ierr = ISRestoreIndices(cis,&col);CHKERRQ(ierr); 146 ierr = ISRestoreIndices(icis,&icol);CHKERRQ(ierr); 147 ierr = ISDestroy(icis);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 EXTERN_C_END 151 152 153