1 2 3 4 #ifdef PETSC_RCS_HEADER 5 static char vcid[] = "$Id: zerodiag.c,v 1.26 1998/11/04 19:40:26 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 we can 89 90 91 @*/ 92 int MatReorderForNonzeroDiagonal(Mat mat,double atol,IS ris,IS cis ) 93 { 94 int ierr, prow, k, nz, n, repl, *j, *col, *row, m, *irow, *icol; 95 Scalar *v; 96 double repla; 97 IS icis,iris; 98 99 PetscFunctionBegin; 100 PetscValidHeaderSpecific(mat,MAT_COOKIE); 101 PetscValidHeaderSpecific(ris,IS_COOKIE); 102 PetscValidHeaderSpecific(cis,IS_COOKIE); 103 104 ierr = ISGetIndices(ris,&row); CHKERRQ(ierr); 105 ierr = ISGetIndices(cis,&col); CHKERRQ(ierr); 106 ierr = ISInvertPermutation(cis,&icis);CHKERRQ(ierr); 107 ierr = ISInvertPermutation(ris,&iris);CHKERRQ(ierr); 108 ierr = ISGetIndices(icis,&icol); CHKERRQ(ierr); 109 ierr = MatGetSize(mat,&m,&n); CHKERRQ(ierr); 110 111 for (prow=0; prow<n; prow++) { 112 ierr = MatGetRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 113 for (k=0; k<nz; k++) {if (icol[j[k]] == prow) break;} 114 if (k >= nz || PetscAbsScalar(v[k]) <= atol) { 115 /* Element too small or zero; find the best candidate */ 116 repl = prow; 117 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 118 /* 119 Here one could sort the icol[j[k]] list to try to select the 120 column closest to the diagonal in the new ordering. (Note have 121 to permute the v[k] values as well, and use a fixed bound on the 122 quality of repla rather then looking for the absolute largest. 123 */ 124 for (k=0; k<nz; k++) { 125 if (icol[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 126 repl = icol[j[k]]; 127 repla = PetscAbsScalar(v[k]); 128 } 129 } 130 if (prow == repl) { 131 /* 132 Look for an element that allows us 133 to pivot with a previous column. To do this, we need 134 to be sure that we don't introduce a zero in a previous 135 diagonal 136 */ 137 if (!MatZeroFindPre_Private(mat,prow,row,col,icol,repla,atol,&repl,&repla,nz,j,v)){ 138 (*PetscErrorPrintf)("Permuted row number %d\n",prow); 139 SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Cannot reorder matrix to eliminate zero diagonal entry"); 140 } 141 } 142 SWAP(icol[col[prow]],icol[col[repl]]); 143 SWAP(col[prow],col[repl]); 144 } 145 ierr = MatRestoreRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 146 } 147 ierr = ISRestoreIndices(ris,&row); CHKERRQ(ierr); 148 ierr = ISRestoreIndices(cis,&col); CHKERRQ(ierr); 149 ierr = ISRestoreIndices(icis,&icol); CHKERRQ(ierr); 150 ierr = ISDestroy(icis); CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 154 155 156