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