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 /*@ 12 MatReorderForNonzeroDiagonal - Changes matrix ordering to remove 13 zeros from diagonal. This may help in the LU factorization to 14 prevent a zero pivot. 15 16 Collective on Mat 17 18 Input Parameters: 19 + mat - matrix to reorder 20 - rmap,cmap - row and column permutations. Usually obtained from 21 MatGetOrdering(). 22 23 Level: intermediate 24 25 Notes: 26 This is not intended as a replacement for pivoting for matrices that 27 have ``bad'' structure. It is only a stop-gap measure. Should be called 28 after a call to MatGetOrdering(), this routine changes the column 29 ordering defined in cis. 30 31 Only works for SeqAIJ matrices 32 33 Options Database Keys (When using KSP): 34 . -pc_factor_nonzeros_along_diagonal - Reorder to remove zeros from diagonal 35 36 Algorithm Notes: 37 Column pivoting is used. 38 39 1) Choice of column is made by looking at the 40 non-zero elements in the troublesome row for columns that are not yet 41 included (moving from left to right). 42 43 2) If (1) fails we check all the columns to the left of the current row 44 and see if one of them has could be swapped. It can be swapped if 45 its corresponding row has a non-zero in the column it is being 46 swapped with; to make sure the previous nonzero diagonal remains 47 nonzero 48 49 @*/ 50 PetscErrorCode MatReorderForNonzeroDiagonal(Mat mat,PetscReal abstol,IS ris,IS cis) 51 { 52 PetscFunctionBegin; 53 PetscTryMethod(mat,"MatReorderForNonzeroDiagonal_C",(Mat,PetscReal,IS,IS),(mat,abstol,ris,cis)); 54 PetscFunctionReturn(0); 55 } 56 57 PETSC_INTERN PetscErrorCode MatGetRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 58 PETSC_INTERN PetscErrorCode MatRestoreRow_SeqAIJ(Mat,PetscInt,PetscInt*,PetscInt**,PetscScalar**); 59 60 #include <../src/vec/is/is/impls/general/general.h> 61 62 PETSC_INTERN PetscErrorCode MatReorderForNonzeroDiagonal_SeqAIJ(Mat mat,PetscReal abstol,IS ris,IS cis) 63 { 64 PetscInt prow,k,nz,n,repl,*j,*col,*row,m,*icol,nnz,*jj,kk; 65 PetscScalar *v,*vv; 66 PetscReal repla; 67 IS icis; 68 69 PetscFunctionBegin; 70 /* access the indices of the IS directly, because it changes them */ 71 row = ((IS_General*)ris->data)->idx; 72 col = ((IS_General*)cis->data)->idx; 73 PetscCall(ISInvertPermutation(cis,PETSC_DECIDE,&icis)); 74 icol = ((IS_General*)icis->data)->idx; 75 PetscCall(MatGetSize(mat,&m,&n)); 76 77 for (prow=0; prow<n; prow++) { 78 PetscCall(MatGetRow_SeqAIJ(mat,row[prow],&nz,&j,&v)); 79 for (k=0; k<nz; k++) { 80 if (icol[j[k]] == prow) break; 81 } 82 if (k >= nz || PetscAbsScalar(v[k]) <= abstol) { 83 /* Element too small or zero; find the best candidate */ 84 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 85 /* 86 Look for a later column we can swap with this one 87 */ 88 for (k=0; k<nz; k++) { 89 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 90 /* found a suitable later column */ 91 repl = icol[j[k]]; 92 SWAP(icol[col[prow]],icol[col[repl]]); 93 SWAP(col[prow],col[repl]); 94 goto found; 95 } 96 } 97 /* 98 Did not find a suitable later column so look for an earlier column 99 We need to be sure that we don't introduce a zero in a previous 100 diagonal 101 */ 102 for (k=0; k<nz; k++) { 103 if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 104 /* See if this one will work */ 105 repl = icol[j[k]]; 106 PetscCall(MatGetRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv)); 107 for (kk=0; kk<nnz; kk++) { 108 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) { 109 PetscCall(MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv)); 110 SWAP(icol[col[prow]],icol[col[repl]]); 111 SWAP(col[prow],col[repl]); 112 goto found; 113 } 114 } 115 PetscCall(MatRestoreRow_SeqAIJ(mat,row[repl],&nnz,&jj,&vv)); 116 } 117 } 118 /* 119 No column suitable; instead check all future rows 120 Note: this will be very slow 121 */ 122 for (k=prow+1; k<n; k++) { 123 PetscCall(MatGetRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv)); 124 for (kk=0; kk<nnz; kk++) { 125 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > abstol) { 126 /* found a row */ 127 SWAP(row[prow],row[k]); 128 goto found; 129 } 130 } 131 PetscCall(MatRestoreRow_SeqAIJ(mat,row[k],&nnz,&jj,&vv)); 132 } 133 134 found:; 135 } 136 PetscCall(MatRestoreRow_SeqAIJ(mat,row[prow],&nz,&j,&v)); 137 } 138 PetscCall(ISDestroy(&icis)); 139 PetscFunctionReturn(0); 140 } 141