30 #ifndef SPARSELU_PIVOTL_H 31 #define SPARSELU_PIVOTL_H 59 template <
typename Scalar,
typename Index>
60 Index SparseLUImpl<Scalar,Index>::pivotL(
const Index jcol,
const RealScalar& diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, Index& pivrow, GlobalLU_t& glu)
63 Index fsupc = (glu.xsup)((glu.supno)(jcol));
64 Index nsupc = jcol - fsupc;
65 Index lptr = glu.xlsub(fsupc);
66 Index nsupr = glu.xlsub(fsupc+1) - lptr;
67 Index lda = glu.xlusup(fsupc+1) - glu.xlusup(fsupc);
68 Scalar* lu_sup_ptr = &(glu.lusup.data()[glu.xlusup(fsupc)]);
69 Scalar* lu_col_ptr = &(glu.lusup.data()[glu.xlusup(jcol)]);
70 Index* lsub_ptr = &(glu.lsub.data()[lptr]);
73 Index diagind = iperm_c(jcol);
74 RealScalar pivmax = 0.0;
76 Index diag = emptyIdxLU;
78 Index isub, icol, itemp, k;
79 for (isub = nsupc; isub < nsupr; ++isub) {
80 rtemp = std::abs(lu_col_ptr[isub]);
85 if (lsub_ptr[isub] == diagind) diag = isub;
89 if ( pivmax == 0.0 ) {
90 pivrow = lsub_ptr[pivptr];
91 perm_r(pivrow) = jcol;
95 RealScalar thresh = diagpivotthresh * pivmax;
104 rtemp = std::abs(lu_col_ptr[diag]);
105 if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag;
107 pivrow = lsub_ptr[pivptr];
111 perm_r(pivrow) = jcol;
113 if (pivptr != nsupc )
115 std::swap( lsub_ptr[pivptr], lsub_ptr[nsupc] );
118 for (icol = 0; icol <= nsupc; icol++)
120 itemp = pivptr + icol * lda;
121 std::swap(lu_sup_ptr[itemp], lu_sup_ptr[nsupc + icol * lda]);
125 Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc];
126 for (k = nsupc+1; k < nsupr; k++)
127 lu_col_ptr[k] *= temp;
134 #endif // SPARSELU_PIVOTL_H
Definition: Eigen_Colamd.h:54