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