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