| 58 | */ |
| 59 | template <typename Scalar, typename StorageIndex> |
| 60 | Index SparseLUImpl<Scalar,StorageIndex>::pivotL(const Index jcol, const RealScalar& diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, Index& pivrow, GlobalLU_t& glu) |
| 61 | { |
| 62 | |
| 63 | Index fsupc = (glu.xsup)((glu.supno)(jcol)); // First column in the supernode containing the column jcol |
| 64 | Index nsupc = jcol - fsupc; // Number of columns in the supernode portion, excluding jcol; nsupc >=0 |
| 65 | Index lptr = glu.xlsub(fsupc); // pointer to the starting location of the row subscripts for this supernode portion |
| 66 | Index nsupr = glu.xlsub(fsupc+1) - lptr; // Number of rows in the supernode |
| 67 | Index lda = glu.xlusup(fsupc+1) - glu.xlusup(fsupc); // leading dimension |
| 68 | Scalar* lu_sup_ptr = &(glu.lusup.data()[glu.xlusup(fsupc)]); // Start of the current supernode |
| 69 | Scalar* lu_col_ptr = &(glu.lusup.data()[glu.xlusup(jcol)]); // Start of jcol in the supernode |
| 70 | StorageIndex* lsub_ptr = &(glu.lsub.data()[lptr]); // Start of row indices of the supernode |
| 71 | |
| 72 | // Determine the largest abs numerical value for partial pivoting |
| 73 | Index diagind = iperm_c(jcol); // diagonal index |
| 74 | RealScalar pivmax(-1.0); |
| 75 | Index pivptr = nsupc; |
| 76 | Index diag = emptyIdxLU; |
| 77 | RealScalar rtemp; |
| 78 | Index isub, icol, itemp, k; |
| 79 | for (isub = nsupc; isub < nsupr; ++isub) { |
| 80 | using std::abs; |
| 81 | rtemp = abs(lu_col_ptr[isub]); |
| 82 | if (rtemp > pivmax) { |
| 83 | pivmax = rtemp; |
| 84 | pivptr = isub; |
| 85 | } |
| 86 | if (lsub_ptr[isub] == diagind) diag = isub; |
| 87 | } |
| 88 | |
| 89 | // Test for singularity |
| 90 | if ( pivmax <= RealScalar(0.0) ) { |
| 91 | // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero |
| 92 | pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr]; |
| 93 | perm_r(pivrow) = StorageIndex(jcol); |
| 94 | return (jcol+1); |
| 95 | } |
| 96 | |
| 97 | RealScalar thresh = diagpivotthresh * pivmax; |
| 98 | |
| 99 | // Choose appropriate pivotal element |
| 100 | |
| 101 | { |
| 102 | // Test if the diagonal element can be used as a pivot (given the threshold value) |
| 103 | if (diag >= 0 ) |
| 104 | { |
| 105 | // Diagonal element exists |
| 106 | using std::abs; |
| 107 | rtemp = abs(lu_col_ptr[diag]); |
| 108 | if (rtemp != RealScalar(0.0) && rtemp >= thresh) pivptr = diag; |
| 109 | } |
| 110 | pivrow = lsub_ptr[pivptr]; |
| 111 | } |
| 112 | |
| 113 | // Record pivot row |
| 114 | perm_r(pivrow) = StorageIndex(jcol); |
| 115 | // Interchange row subscripts |
| 116 | if (pivptr != nsupc ) |
| 117 | { |