1 2 3 4 #ifdef PETSC_RCS_HEADER 5 static char vcid[] = "$Id: zerodiag.c,v 1.27 1998/11/04 19:45:13 bsmith Exp bsmith $"; 6 #endif 7 8 /* 9 This file contains routines to reorder a matrix so that the diagonal 10 elements are nonzero. 11 */ 12 13 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 14 #include <math.h> 15 16 #define SWAP(a,b) {int _t; _t = a; a = b; b = _t; } 17 18 #undef __FUNC__ 19 #define __FUNC__ "MatZeroFindPre_Private" 20 /* 21 Given a current row and current permutation, find a column permutation 22 that removes a zero diagonal. 23 */ 24 int MatZeroFindPre_Private(Mat mat,int prow,int* row,int* col,int *icol,double repla, 25 double atol,int* rc,double* rcv,int nz, int *j, Scalar *v) 26 { 27 int k, repl, kk, nnz, *jj,ierr; 28 Scalar *vv; 29 30 PetscFunctionBegin; 31 /* 32 Here one could sort the col[j[k]] to try to select the column closest 33 to the diagonal (in the new ordering) that satisfies the criteria 34 */ 35 for (k=0; k<nz; k++) { 36 if (icol[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 37 /* See if this one will work */ 38 repl = icol[j[k]]; 39 ierr = MatGetRow( mat, row[repl], &nnz, &jj, &vv ); CHKERRQ(ierr); 40 for (kk=0; kk<nnz; kk++) { 41 if (icol[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 42 *rcv = PetscAbsScalar(v[k]); 43 *rc = repl; 44 ierr = MatRestoreRow( mat, row[repl], &nnz, &jj, &vv ); CHKERRQ(ierr); 45 PetscFunctionReturn(1); 46 } 47 } 48 ierr = MatRestoreRow( mat, row[repl], &nnz, &jj, &vv ); CHKERRQ(ierr); 49 } 50 } 51 52 SETERRQ(1,1," "); 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNC__ 57 #define __FUNC__ "MatReorderForNonzeroDiagonal" 58 /*@ 59 MatReorderForNonzeroDiagonal - Changes matrix ordering to remove 60 zeros from diagonal. This may help in the LU factorization to 61 prevent a zero pivot. 62 63 Collective on Mat 64 65 Input Parameters: 66 + mat - matrix to reorder 67 - rmap,cmap - row and column permutations. Usually obtained from 68 MatGetReordering(). 69 70 Notes: 71 This is not intended as a replacement for pivoting for matrices that 72 have ``bad'' structure. It is only a stop-gap measure. Should be called 73 after a call to MatGetReordering(), this routine changes the column 74 ordering defined in cis. 75 76 Options Database Keys (When using SLES): 77 + -pc_ilu_nonzeros_along_diagonal 78 - -pc_lu_nonzeros_along_diagonal 79 80 Algorithm Notes: 81 Column pivoting is used. 82 83 1) Choice of column is made by looking at the 84 non-zero elements in the troublesome row for columns that are not yet 85 included (moving from left to right). 86 87 2) If (1) fails we check all the columns to the left of the current row 88 and see if one of them has could be swapped. It can be swapped if 89 its corresponding row has a non-zero in the column it is being 90 swapped with; to make sure the previous nonzero diagonal remains 91 nonzero 92 93 94 @*/ 95 int MatReorderForNonzeroDiagonal(Mat mat,double atol,IS ris,IS cis ) 96 { 97 int ierr, prow, k, nz, n, repl, *j, *col, *row, m, *irow, *icol; 98 Scalar *v; 99 double repla; 100 IS icis,iris; 101 102 PetscFunctionBegin; 103 PetscValidHeaderSpecific(mat,MAT_COOKIE); 104 PetscValidHeaderSpecific(ris,IS_COOKIE); 105 PetscValidHeaderSpecific(cis,IS_COOKIE); 106 107 ierr = ISGetIndices(ris,&row); CHKERRQ(ierr); 108 ierr = ISGetIndices(cis,&col); CHKERRQ(ierr); 109 ierr = ISInvertPermutation(cis,&icis);CHKERRQ(ierr); 110 ierr = ISInvertPermutation(ris,&iris);CHKERRQ(ierr); 111 ierr = ISGetIndices(icis,&icol); CHKERRQ(ierr); 112 ierr = MatGetSize(mat,&m,&n); CHKERRQ(ierr); 113 114 for (prow=0; prow<n; prow++) { 115 ierr = MatGetRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 116 for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;} 117 if (k >= nz || PetscAbsScalar(v[k]) <= atol) { 118 /* Element too small or zero; find the best candidate */ 119 repl = prow; 120 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 121 /* 122 Here one could sort the icol[j[k]] list to try to select the 123 column closest to the diagonal in the new ordering. (Note have 124 to permute the v[k] values as well, and use a fixed bound on the 125 quality of repla rather then looking for the absolute largest. 126 */ 127 for (k=0; k<nz; k++) { 128 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 129 repl = icol[j[k]]; 130 repla = PetscAbsScalar(v[k]); 131 } 132 } 133 if (prow == repl) { 134 /* 135 Look for an element that allows us 136 to pivot with a previous column. To do this, we need 137 to be sure that we don't introduce a zero in a previous 138 diagonal 139 */ 140 if (!MatZeroFindPre_Private(mat,prow,row,col,icol,repla,atol,&repl,&repla,nz,j,v)){ 141 (*PetscErrorPrintf)("Permuted row number %d\n",prow); 142 SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Cannot reorder matrix to eliminate zero diagonal entry"); 143 } 144 } 145 SWAP(icol[col[prow]],icol[col[repl]]); 146 SWAP(col[prow],col[repl]); 147 } 148 ierr = MatRestoreRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 149 } 150 ierr = ISRestoreIndices(ris,&row); CHKERRQ(ierr); 151 ierr = ISRestoreIndices(cis,&col); CHKERRQ(ierr); 152 ierr = ISRestoreIndices(icis,&icol); CHKERRQ(ierr); 153 ierr = ISDestroy(icis); CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 158 159