1 2 /* 3 This file contains routines to reorder a matrix so that the diagonal 4 elements are nonzero. 5 */ 6 7 #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/ 8 9 #define SWAP(a,b) {PetscInt _t; _t = a; a = b; b = _t; } 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatReorderForNonzeroDiagonal" 13 /*@ 14 MatReorderForNonzeroDiagonal - Changes matrix ordering to remove 15 zeros from diagonal. This may help in the LU factorization to 16 prevent a zero pivot. 17 18 Collective on Mat 19 20 Input Parameters: 21 + mat - matrix to reorder 22 - rmap,cmap - row and column permutations. Usually obtained from 23 MatGetOrdering(). 24 25 Level: intermediate 26 27 Notes: 28 This is not intended as a replacement for pivoting for matrices that 29 have ``bad'' structure. It is only a stop-gap measure. Should be called 30 after a call to MatGetOrdering(), this routine changes the column 31 ordering defined in cis. 32 33 Only works for SeqAIJ matrices 34 35 Options Database Keys (When using KSP): 36 . -pc_factor_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 PetscErrorCode MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis) 54 { 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 ierr = PetscTryMethod(mat,"MatReorderForNonzeroDiagonal_C",(Mat,PetscReal,IS,IS),(mat,abstol,ris,cis));CHKERRQ(ierr); 59 PetscFunctionReturn(0); 60 } 61 62 PETSC_INTERN PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 63 PETSC_INTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 64 65 #include <../src/vec/is/is/impls/general/general.h> 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "MatReorderForNonzeroDiagonal_SeqAIJ" 69 PETSC_INTERN PetscErrorCode MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis) 70 { 71 PetscErrorCode ierr; 72 PetscInt prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk; 73 PetscScalar *v,*vv; 74 PetscReal repla; 75 IS icis; 76 77 PetscFunctionBegin; 78 /* access the indices of the IS directly, because it changes them */ 79 row = ((IS_General*)ris->data)->idx; 80 col = ((IS_General*)cis->data)->idx; 81 ierr = ISInvertPermutation(cis,PETSC_DECIDE,&icis);CHKERRQ(ierr); 82 icol = ((IS_General*)icis->data)->idx; 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++) { 88 if (icol[j[k]] == prow) break; 89 } 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 ierr = MatGetRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 115 for (kk=0; kk<nnz; kk++) { 116 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) { 117 ierr = MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 118 SWAP(icol[col[prow]],icol[col[repl]]); 119 SWAP(col[prow],col[repl]); 120 goto found; 121 } 122 } 123 ierr = MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv);CHKERRQ(ierr); 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 ierr = MatGetRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 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 ierr = MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv);CHKERRQ(ierr); 140 } 141 142 found:; 143 } 144 ierr = MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v);CHKERRQ(ierr); 145 } 146 ierr = ISDestroy(&icis);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 151 152