Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/Eigen/src/SparseCholesky/CMakeLists.txt b/Eigen/src/SparseCholesky/CMakeLists.txt
deleted file mode 100644
index 375a59d..0000000
--- a/Eigen/src/SparseCholesky/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_SparseCholesky_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_SparseCholesky_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseCholesky COMPONENT Devel
- )
diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h
index e1f96ba..2907f65 100644
--- a/Eigen/src/SparseCholesky/SimplicialCholesky.h
+++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h
@@ -17,43 +17,74 @@
SimplicialCholeskyLDLT
};
+namespace internal {
+ template<typename CholMatrixType, typename InputMatrixType>
+ struct simplicial_cholesky_grab_input {
+ typedef CholMatrixType const * ConstCholMatrixPtr;
+ static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
+ {
+ tmp = input;
+ pmat = &tmp;
+ }
+ };
+
+ template<typename MatrixType>
+ struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
+ typedef MatrixType const * ConstMatrixPtr;
+ static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/)
+ {
+ pmat = &input;
+ }
+ };
+} // end namespace internal
+
/** \ingroup SparseCholesky_Module
- * \brief A direct sparse Cholesky factorizations
+ * \brief A base class for direct sparse Cholesky factorizations
*
- * These classes provide LL^T and LDL^T Cholesky factorizations of sparse matrices that are
- * selfadjoint and positive definite. The factorization allows for solving A.X = B where
+ * This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are
+ * selfadjoint and positive definite. These factorizations allow for solving A.X = B where
* X and B can be either dense or sparse.
*
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
* such that the factorized matrix is P A P^-1.
*
- * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
- * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
- * or Upper. Default is Lower.
+ * \tparam Derived the type of the derived class, that is the actual factorization type.
*
*/
template<typename Derived>
-class SimplicialCholeskyBase : internal::noncopyable
+class SimplicialCholeskyBase : public SparseSolverBase<Derived>
{
+ typedef SparseSolverBase<Derived> Base;
+ using Base::m_isInitialized;
+
public:
typedef typename internal::traits<Derived>::MatrixType MatrixType;
typedef typename internal::traits<Derived>::OrderingType OrderingType;
enum { UpLo = internal::traits<Derived>::UpLo };
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
+ typedef CholMatrixType const * ConstCholMatrixPtr;
typedef Matrix<Scalar,Dynamic,1> VectorType;
+ typedef Matrix<StorageIndex,Dynamic,1> VectorI;
+
+ enum {
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+ };
public:
+
+ using Base::derived;
/** Default constructor */
SimplicialCholeskyBase()
- : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
+ : m_info(Success), m_shiftOffset(0), m_shiftScale(1)
{}
- SimplicialCholeskyBase(const MatrixType& matrix)
- : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
+ explicit SimplicialCholeskyBase(const MatrixType& matrix)
+ : m_info(Success), m_shiftOffset(0), m_shiftScale(1)
{
derived().compute(matrix);
}
@@ -79,42 +110,14 @@
return m_info;
}
- /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
- *
- * \sa compute()
- */
- template<typename Rhs>
- inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
- solve(const MatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
- eigen_assert(rows()==b.rows()
- && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
- }
-
- /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
- *
- * \sa compute()
- */
- template<typename Rhs>
- inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
- solve(const SparseMatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
- eigen_assert(rows()==b.rows()
- && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
- return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
- }
-
/** \returns the permutation P
* \sa permutationPinv() */
- const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
+ const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationP() const
{ return m_P; }
/** \returns the inverse P^-1 of the permutation P
* \sa permutationP() */
- const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const
+ const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationPinv() const
{ return m_Pinv; }
/** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization.
@@ -150,7 +153,7 @@
/** \internal */
template<typename Rhs,typename Dest>
- void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
+ void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
{
eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
eigen_assert(m_matrix.rows()==b.rows());
@@ -175,6 +178,12 @@
if(m_P.size()>0)
dest = m_Pinv * dest;
}
+
+ template<typename Rhs,typename Dest>
+ void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
+ {
+ internal::solve_sparse_through_dense_panels(derived(), b, dest);
+ }
#endif // EIGEN_PARSED_BY_DOXYGEN
@@ -186,20 +195,33 @@
{
eigen_assert(matrix.rows()==matrix.cols());
Index size = matrix.cols();
- CholMatrixType ap(size,size);
- ordering(matrix, ap);
- analyzePattern_preordered(ap, DoLDLT);
- factorize_preordered<DoLDLT>(ap);
+ CholMatrixType tmp(size,size);
+ ConstCholMatrixPtr pmat;
+ ordering(matrix, pmat, tmp);
+ analyzePattern_preordered(*pmat, DoLDLT);
+ factorize_preordered<DoLDLT>(*pmat);
}
template<bool DoLDLT>
void factorize(const MatrixType& a)
{
eigen_assert(a.rows()==a.cols());
- int size = a.cols();
- CholMatrixType ap(size,size);
- ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
- factorize_preordered<DoLDLT>(ap);
+ Index size = a.cols();
+ CholMatrixType tmp(size,size);
+ ConstCholMatrixPtr pmat;
+
+ if(m_P.size()==0 && (UpLo&Upper)==Upper)
+ {
+ // If there is no ordering, try to directly use the input matrix without any copy
+ internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
+ }
+ else
+ {
+ tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
+ pmat = &tmp;
+ }
+
+ factorize_preordered<DoLDLT>(*pmat);
}
template<bool DoLDLT>
@@ -208,14 +230,15 @@
void analyzePattern(const MatrixType& a, bool doLDLT)
{
eigen_assert(a.rows()==a.cols());
- int size = a.cols();
- CholMatrixType ap(size,size);
- ordering(a, ap);
- analyzePattern_preordered(ap,doLDLT);
+ Index size = a.cols();
+ CholMatrixType tmp(size,size);
+ ConstCholMatrixPtr pmat;
+ ordering(a, pmat, tmp);
+ analyzePattern_preordered(*pmat,doLDLT);
}
void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
- void ordering(const MatrixType& a, CholMatrixType& ap);
+ void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
/** keeps off-diagonal entries; drops diagonal entries */
struct keep_diag {
@@ -226,24 +249,23 @@
};
mutable ComputationInfo m_info;
- bool m_isInitialized;
bool m_factorizationIsOk;
bool m_analysisIsOk;
CholMatrixType m_matrix;
VectorType m_diag; // the diagonal coefficients (LDLT mode)
- VectorXi m_parent; // elimination tree
- VectorXi m_nonZerosPerCol;
- PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation
- PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation
+ VectorI m_parent; // elimination tree
+ VectorI m_nonZerosPerCol;
+ PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // the permutation
+ PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation
RealScalar m_shiftOffset;
RealScalar m_shiftScale;
};
-template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
-template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
-template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky;
+template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT;
+template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT;
+template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky;
namespace internal {
@@ -253,12 +275,12 @@
typedef _Ordering OrderingType;
enum { UpLo = _UpLo };
typedef typename MatrixType::Scalar Scalar;
- typedef typename MatrixType::Index Index;
- typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
- typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
- typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
- static inline MatrixL getL(const MatrixType& m) { return m; }
- static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
+ typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
+ typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
+ static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
+ static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
};
template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
@@ -267,12 +289,12 @@
typedef _Ordering OrderingType;
enum { UpLo = _UpLo };
typedef typename MatrixType::Scalar Scalar;
- typedef typename MatrixType::Index Index;
- typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
- typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
- typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
- static inline MatrixL getL(const MatrixType& m) { return m; }
- static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
+ typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
+ typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
+ static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
+ static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
};
template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
@@ -300,6 +322,8 @@
* or Upper. Default is Lower.
* \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
*
+ * \implsparsesolverconcept
+ *
* \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering
*/
template<typename _MatrixType, int _UpLo, typename _Ordering>
@@ -311,7 +335,7 @@
typedef SimplicialCholeskyBase<SimplicialLLT> Base;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::StorageIndex StorageIndex;
typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef internal::traits<SimplicialLLT> Traits;
@@ -321,7 +345,7 @@
/** Default constructor */
SimplicialLLT() : Base() {}
/** Constructs and performs the LLT factorization of \a matrix */
- SimplicialLLT(const MatrixType& matrix)
+ explicit SimplicialLLT(const MatrixType& matrix)
: Base(matrix) {}
/** \returns an expression of the factor L */
@@ -389,6 +413,8 @@
* or Upper. Default is Lower.
* \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
*
+ * \implsparsesolverconcept
+ *
* \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering
*/
template<typename _MatrixType, int _UpLo, typename _Ordering>
@@ -400,8 +426,8 @@
typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef internal::traits<SimplicialLDLT> Traits;
typedef typename Traits::MatrixL MatrixL;
@@ -411,7 +437,7 @@
SimplicialLDLT() : Base() {}
/** Constructs and performs the LLT factorization of \a matrix */
- SimplicialLDLT(const MatrixType& matrix)
+ explicit SimplicialLDLT(const MatrixType& matrix)
: Base(matrix) {}
/** \returns a vector expression of the diagonal D */
@@ -482,8 +508,8 @@
typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef internal::traits<SimplicialCholesky> Traits;
typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
@@ -491,7 +517,7 @@
public:
SimplicialCholesky() : Base(), m_LDLT(true) {}
- SimplicialCholesky(const MatrixType& matrix)
+ explicit SimplicialCholesky(const MatrixType& matrix)
: Base(), m_LDLT(true)
{
compute(matrix);
@@ -560,7 +586,7 @@
/** \internal */
template<typename Rhs,typename Dest>
- void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
+ void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
{
eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
eigen_assert(Base::m_matrix.rows()==b.rows());
@@ -596,6 +622,13 @@
dest = Base::m_Pinv * dest;
}
+ /** \internal */
+ template<typename Rhs,typename Dest>
+ void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
+ {
+ internal::solve_sparse_through_dense_panels(*this, b, dest);
+ }
+
Scalar determinant() const
{
if(m_LDLT)
@@ -614,58 +647,43 @@
};
template<typename Derived>
-void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
+void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
{
eigen_assert(a.rows()==a.cols());
const Index size = a.rows();
- // Note that amd compute the inverse permutation
+ pmat = ≈
+ // Note that ordering methods compute the inverse permutation
+ if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
{
- CholMatrixType C;
- C = a.template selfadjointView<UpLo>();
+ {
+ CholMatrixType C;
+ C = a.template selfadjointView<UpLo>();
+
+ OrderingType ordering;
+ ordering(C,m_Pinv);
+ }
+
+ if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
+ else m_P.resize(0);
- OrderingType ordering;
- ordering(C,m_Pinv);
+ ap.resize(size,size);
+ ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
}
-
- if(m_Pinv.size()>0)
- m_P = m_Pinv.inverse();
else
+ {
+ m_Pinv.resize(0);
m_P.resize(0);
-
- ap.resize(size,size);
- ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
+ if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
+ {
+ // we have to transpose the lower part to to the upper one
+ ap.resize(size,size);
+ ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
+ }
+ else
+ internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
+ }
}
-namespace internal {
-
-template<typename Derived, typename Rhs>
-struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
- : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
-{
- typedef SimplicialCholeskyBase<Derived> Dec;
- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec().derived()._solve(rhs(),dst);
- }
-};
-
-template<typename Derived, typename Rhs>
-struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
- : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
-{
- typedef SimplicialCholeskyBase<Derived> Dec;
- EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- this->defaultEvalTo(dst);
- }
-};
-
-} // end namespace internal
-
} // end namespace Eigen
#endif // EIGEN_SIMPLICIAL_CHOLESKY_H
diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
index 7aaf702..31e0699 100644
--- a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
+++ b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
@@ -50,14 +50,14 @@
template<typename Derived>
void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT)
{
- const Index size = ap.rows();
+ const StorageIndex size = StorageIndex(ap.rows());
m_matrix.resize(size, size);
m_parent.resize(size);
m_nonZerosPerCol.resize(size);
- ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
+ ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);
- for(Index k = 0; k < size; ++k)
+ for(StorageIndex k = 0; k < size; ++k)
{
/* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
m_parent[k] = -1; /* parent of k is not yet known */
@@ -65,7 +65,7 @@
m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
{
- Index i = it.index();
+ StorageIndex i = it.index();
if(i < k)
{
/* follow path from i to root of etree, stop at flagged node */
@@ -82,9 +82,9 @@
}
/* construct Lp index array from m_nonZerosPerCol column counts */
- Index* Lp = m_matrix.outerIndexPtr();
+ StorageIndex* Lp = m_matrix.outerIndexPtr();
Lp[0] = 0;
- for(Index k = 0; k < size; ++k)
+ for(StorageIndex k = 0; k < size; ++k)
Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
m_matrix.resizeNonZeros(Lp[size]);
@@ -104,31 +104,31 @@
eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
eigen_assert(ap.rows()==ap.cols());
- const Index size = ap.rows();
- eigen_assert(m_parent.size()==size);
- eigen_assert(m_nonZerosPerCol.size()==size);
+ eigen_assert(m_parent.size()==ap.rows());
+ eigen_assert(m_nonZerosPerCol.size()==ap.rows());
- const Index* Lp = m_matrix.outerIndexPtr();
- Index* Li = m_matrix.innerIndexPtr();
+ const StorageIndex size = StorageIndex(ap.rows());
+ const StorageIndex* Lp = m_matrix.outerIndexPtr();
+ StorageIndex* Li = m_matrix.innerIndexPtr();
Scalar* Lx = m_matrix.valuePtr();
ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
- ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0);
- ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
+ ei_declare_aligned_stack_constructed_variable(StorageIndex, pattern, size, 0);
+ ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);
bool ok = true;
m_diag.resize(DoLDLT ? size : 0);
- for(Index k = 0; k < size; ++k)
+ for(StorageIndex k = 0; k < size; ++k)
{
// compute nonzero pattern of kth row of L, in topological order
y[k] = 0.0; // Y(0:k) is now all zero
- Index top = size; // stack for pattern is empty
+ StorageIndex top = size; // stack for pattern is empty
tags[k] = k; // mark node k as visited
m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L
- for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
+ for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
{
- Index i = it.index();
+ StorageIndex i = it.index();
if(i <= k)
{
y[i] += numext::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */