1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: zerodiag.c,v 1.21 1998/04/15 22:52:39 curfman 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,double repla, 22 double atol,int* rc,double* rcv ) 23 { 24 int k, nz, repl, *j, kk, nnz, *jj,ierr; 25 Scalar *v, *vv; 26 27 PetscFunctionBegin; 28 ierr = MatGetRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 29 /* 30 Here one could sort the col[j[k]] to try to select the column closest 31 to the diagonal (in the new ordering) that satisfies the criteria 32 */ 33 for (k=0; k<nz; k++) { 34 if (col[j[k]] < prow && PetscAbsScalar(v[k]) > repla) { 35 /* See if this one will work */ 36 repl = col[j[k]]; 37 ierr = MatGetRow( mat, row[repl], &nnz, &jj, &vv ); CHKERRQ(ierr); 38 for (kk=0; kk<nnz; kk++) { 39 if (col[jj[kk]] == prow && PetscAbsScalar(vv[kk]) > atol) { 40 *rcv = PetscAbsScalar(v[k]); 41 *rc = repl; 42 ierr = MatRestoreRow( mat, row[repl], &nnz, &jj, &vv ); CHKERRQ(ierr); 43 ierr = MatRestoreRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 44 PetscFunctionReturn(1); 45 } 46 } 47 ierr = MatRestoreRow( mat, row[repl], &nnz, &jj, &vv ); CHKERRQ(ierr); 48 } 49 } 50 ierr = MatRestoreRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNC__ 55 #define __FUNC__ "MatReorderForNonzeroDiagonal" 56 /*@ 57 MatReorderForNonzeroDiagonal - Changes matrix ordering to remove 58 zeros from diagonal. This may help in the LU factorization to 59 prevent a zero pivot. 60 61 Collective on Mat 62 63 Input Parameters: 64 + mat - matrix to reorder 65 - rmap,cmap - row and column permutations. Usually obtained from 66 MatGetReordering(). 67 68 Notes: 69 This is not intended as a replacement for pivoting for matrices that 70 have ``bad'' structure. It is only a stop-gap measure. Should be called 71 after a call to MatGetReordering(), this routine changes the column 72 ordering defined in cis. 73 74 Options Database Keys (When using SLES): 75 + -pc_ilu_nonzeros_along_diagonal 76 - -pc_lu_nonzeros_along_diagonal 77 78 Algorithm: 79 Column pivoting is used. Choice of column is made by looking at the 80 non-zero elements in the row. This algorithm is simple and fast but 81 does NOT guarantee that a non-singular or well conditioned 82 principle submatrix will be produced. 83 84 @*/ 85 int MatReorderForNonzeroDiagonal(Mat mat,double atol,IS ris,IS cis ) 86 { 87 int ierr, prow, k, nz, n, repl, *j, *col, *row, m; 88 Scalar *v; 89 double repla; 90 91 PetscFunctionBegin; 92 PetscValidHeaderSpecific(mat,MAT_COOKIE); 93 PetscValidHeaderSpecific(ris,IS_COOKIE); 94 PetscValidHeaderSpecific(cis,IS_COOKIE); 95 96 ierr = ISGetIndices(ris,&row); CHKERRQ(ierr); 97 ierr = ISGetIndices(cis,&col); CHKERRQ(ierr); 98 ierr = MatGetSize(mat,&m,&n); CHKERRQ(ierr); 99 100 for (prow=0; prow<n; prow++) { 101 ierr = MatGetRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 102 for (k=0; k<nz; k++) {if (col[j[k]] == prow) break;} 103 if (k >= nz || PetscAbsScalar(v[k]) <= atol) { 104 /* Element too small or zero; find the best candidate */ 105 repl = prow; 106 repla = (k >= nz) ? 0.0 : PetscAbsScalar(v[k]); 107 /* 108 Here one could sort the col[j[k]] list to try to select the 109 column closest to the diagonal in the new ordering. (Note have 110 to permute the v[k] values as well, and use a fixed bound on the 111 quality of repla rather then looking for the absolute largest. 112 */ 113 for (k=0; k<nz; k++) { 114 if (col[j[k]] > prow && PetscAbsScalar(v[k]) > repla) { 115 repl = col[j[k]]; 116 repla = PetscAbsScalar(v[k]); 117 } 118 } 119 if (prow == repl) { 120 /* Look for an element that allows us 121 to pivot with a previous column. To do this, we need 122 to be sure that we don't introduce a zero in a previous 123 diagonal */ 124 if (!MatZeroFindPre_Private(mat,prow,row,col,repla,atol,&repl,&repla)){ 125 SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Cannot reorder matrix to eliminate zero diagonal entry"); 126 } 127 } 128 SWAP(col[prow],col[repl]); 129 } 130 ierr = MatRestoreRow( mat, row[prow], &nz, &j, &v ); CHKERRQ(ierr); 131 } 132 ierr = ISRestoreIndices(ris,&row); CHKERRQ(ierr); 133 ierr = ISRestoreIndices(cis,&col); CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 138 139