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