diff options
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers')
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 74 |
1 files changed, 38 insertions, 36 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index c413e9e1a..6d63d45e4 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -25,14 +25,14 @@ namespace internal { * \param ind The array of index for the elements in @p row * \param ncut The number of largest elements to keep **/ -template <typename VectorV, typename VectorI, typename Index> +template <typename VectorV, typename VectorI> Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) { typedef typename VectorV::RealScalar RealScalar; using std::swap; using std::abs; Index mid; - Index n = convert_index<Index>(row.size()); /* length of the vector */ + Index n = row.size(); /* length of the vector */ Index first, last ; ncut--; /* to fit the zero-based indices */ @@ -124,9 +124,9 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar> > compute(mat); } - StorageIndex rows() const { return m_lu.rows(); } + Index rows() const { return m_lu.rows(); } - StorageIndex cols() const { return m_lu.cols(); } + Index cols() const { return m_lu.cols(); } /** \brief Reports whether previous computation was successful. * @@ -239,9 +239,10 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) using std::sqrt; using std::swap; using std::abs; + using internal::convert_index; eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix"); - StorageIndex n = amat.cols(); // Size of the matrix + Index n = amat.cols(); // Size of the matrix m_lu.resize(n,n); // Declare Working vectors and variables Vector u(n) ; // real values of the row -- maximum size is n -- @@ -259,36 +260,36 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) u.fill(0); // number of largest elements to keep in each row: - StorageIndex fill_in = static_cast<StorageIndex> (amat.nonZeros()*m_fillfactor)/n+1; + Index fill_in = (amat.nonZeros()*m_fillfactor)/n + 1; if (fill_in > n) fill_in = n; // number of largest nonzero elements to keep in the L and the U part of the current row: - StorageIndex nnzL = fill_in/2; - StorageIndex nnzU = nnzL; + Index nnzL = fill_in/2; + Index nnzU = nnzL; m_lu.reserve(n * (nnzL + nnzU + 1)); // global loop over the rows of the sparse matrix - for (StorageIndex ii = 0; ii < n; ii++) + for (Index ii = 0; ii < n; ii++) { // 1 - copy the lower and the upper part of the row i of mat in the working vector u - StorageIndex sizeu = 1; // number of nonzero elements in the upper part of the current row - StorageIndex sizel = 0; // number of nonzero elements in the lower part of the current row - ju(ii) = ii; + Index sizeu = 1; // number of nonzero elements in the upper part of the current row + Index sizel = 0; // number of nonzero elements in the lower part of the current row + ju(ii) = convert_index<StorageIndex>(ii); u(ii) = 0; - jr(ii) = ii; + jr(ii) = convert_index<StorageIndex>(ii); RealScalar rownorm = 0; typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii for (; j_it; ++j_it) { - StorageIndex k = j_it.index(); + Index k = j_it.index(); if (k < ii) { // copy the lower part - ju(sizel) = k; + ju(sizel) = convert_index<StorageIndex>(k); u(sizel) = j_it.value(); - jr(k) = sizel; + jr(k) = convert_index<StorageIndex>(sizel); ++sizel; } else if (k == ii) @@ -298,10 +299,10 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) else { // copy the upper part - StorageIndex jpos = ii + sizeu; - ju(jpos) = k; + Index jpos = ii + sizeu; + ju(jpos) = convert_index<StorageIndex>(k); u(jpos) = j_it.value(); - jr(k) = jpos; + jr(k) = convert_index<StorageIndex>(jpos); ++sizeu; } rownorm += numext::abs2(j_it.value()); @@ -317,21 +318,22 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) rownorm = sqrt(rownorm); // 3 - eliminate the previous nonzero rows - StorageIndex jj = 0; - StorageIndex len = 0; + Index jj = 0; + Index len = 0; while (jj < sizel) { // In order to eliminate in the correct order, // we must select first the smallest column index among ju(jj:sizel) - StorageIndex k; - StorageIndex minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment + Index k; + Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment k += jj; if (minrow != ju(jj)) { // swap the two locations - StorageIndex j = ju(jj); + Index j = ju(jj); swap(ju(jj), ju(k)); - jr(minrow) = jj; jr(j) = k; + jr(minrow) = convert_index<StorageIndex>(jj); + jr(j) = convert_index<StorageIndex>(k); swap(u(jj), u(k)); } // Reset this location @@ -355,11 +357,11 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) for (; ki_it; ++ki_it) { Scalar prod = fact * ki_it.value(); - StorageIndex j = ki_it.index(); - StorageIndex jpos = jr(j); + Index j = ki_it.index(); + Index jpos = jr(j); if (jpos == -1) // fill-in element { - StorageIndex newpos; + Index newpos; if (j >= ii) // dealing with the upper part { newpos = ii + sizeu; @@ -372,23 +374,23 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) sizel++; eigen_internal_assert(sizel<=ii); } - ju(newpos) = j; + ju(newpos) = convert_index<StorageIndex>(j); u(newpos) = -prod; - jr(j) = newpos; + jr(j) = convert_index<StorageIndex>(newpos); } else u(jpos) -= prod; } // store the pivot element - u(len) = fact; - ju(len) = minrow; + u(len) = fact; + ju(len) = convert_index<StorageIndex>(minrow); ++len; jj++; } // end of the elimination on the row ii // reset the upper part of the pointer jr to zero - for(StorageIndex k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1; + for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1; // 4 - partially sort and insert the elements in the m_lu matrix @@ -401,7 +403,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) // store the largest m_fill elements of the L part m_lu.startVec(ii); - for(StorageIndex k = 0; k < len; k++) + for(Index k = 0; k < len; k++) m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); // store the diagonal element @@ -413,7 +415,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) // sort the U-part of the row // apply the dropping rule first len = 0; - for(StorageIndex k = 1; k < sizeu; k++) + for(Index k = 1; k < sizeu; k++) { if(abs(u(ii+k)) > m_droptol * rownorm ) { @@ -429,7 +431,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) internal::QuickSplit(uu, juu, len); // store the largest elements of the U part - for(StorageIndex k = ii + 1; k < ii + len; k++) + for(Index k = ii + 1; k < ii + len; k++) m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); } |