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