1 #define PETSCMAT_DLL 2 3 /* 4 This file contains routines to reorder a matrix so that the diagonal 5 elements are nonzero. 6 */ 7 8 #include "include/private/matimpl.h" /*I "petscmat.h" I*/ 9 10 #define SWAP(a,b) {PetscInt _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 KSP): 37 . -pc_factor_nonzeros_along_diagonal 38 39 Algorithm Notes: 40 Column pivoting is used. 41 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 51 52 53 @*/ 54 PetscErrorCode PETSCMAT_DLLEXPORT MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis) 55 { 56 PetscErrorCode ierr,(*f)(Mat,PetscReal,IS,IS); 57 58 PetscFunctionBegin; 59 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatReorderForNonzeroDiagonal_C",(void (**)(void))&f);CHKERRQ(ierr); 60 if (f) { 61 ierr = (*f)(mat,abstol,ris,cis);CHKERRQ(ierr); 62 } 63 PetscFunctionReturn(0); 64 } 65 66 EXTERN PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 67 EXTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 68 69 #include "src/vec/is/impls/general/general.h" 70 71 EXTERN_C_BEGIN 72 #undef __FUNCT__ 73 #define __FUNCT__ "MatReorderForNonzeroDiagonal_SeqAIJ" 74 PetscErrorCode PETSCMAT_DLLEXPORT MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis) 75 { 76 PetscErrorCode ierr; 77 PetscInt prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk; 78 PetscScalar *v,*vv; 79 PetscReal repla; 80 IS icis; 81 82 PetscFunctionBegin; 83 /* access the indices of the IS directly, because it changes them */ 84 row = ((IS_General*)ris->data)->idx; 85 col = ((IS_General*)cis->data)->idx; 86 ierr = ISInvertPermutation(cis,PETSC_DECIDE,&icis);CHKERRQ(ierr); 87 icol = ((IS_General*)icis->data)->idx; 88 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 89 90 for (prow=0; prow<n; prow++) { 91 ierr = MatGetRow_SeqAIJ(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 92 for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;} 93 if (k >= nz || PetscAbsScalar(v[k]) <= abstol) { 94 /* Element too small or zero; find the best candidate */ 95 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 96 /* 97 Look for a later column we can swap with this one 98 */ 99 for (k=0; k<nz; k++) { 100 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 101 /* found a suitable later column */ 102 repl = icol[j[k]]; 103 SWAP(icol[col[prow]],icol[col[repl]]); 104 SWAP(col[prow],col[repl]); 105 goto found; 106 } 107 } 108 /* 109 Did not find a suitable later column so look for an earlier column 110 We need to be sure that we don't introduce a zero in a previous 111 diagonal 112 */ 113 for (k=0; k<nz; k++) { 114 if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 115 /* See if this one will work */ 116 repl = icol[j[k]]; 117 ierr = MatGetRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 118 for (kk=0; kk<nnz; kk++) { 119 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) { 120 ierr = MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 121 SWAP(icol[col[prow]],icol[col[repl]]); 122 SWAP(col[prow],col[repl]); 123 goto found; 124 } 125 } 126 ierr = MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 127 } 128 } 129 /* 130 No column suitable; instead check all future rows 131 Note: this will be very slow 132 */ 133 for (k=prow+1; k<n; k++) { 134 ierr = MatGetRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 135 for (kk=0; kk<nnz; kk++) { 136 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) { 137 /* found a row */ 138 SWAP(row[prow],row[k]); 139 goto found; 140 } 141 } 142 ierr = MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 143 } 144 145 found:; 146 } 147 ierr = MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 148 } 149 ierr = ISDestroy(icis);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 EXTERN_C_END 153 154 155