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