aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/Cholesky6
-rw-r--r--Eigen/src/Cholesky/Cholesky.h56
-rw-r--r--Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h65
-rw-r--r--Eigen/src/Cholesky/LDLT.h202
-rw-r--r--Eigen/src/Cholesky/LLT.h186
-rw-r--r--Eigen/src/Core/MatrixBase.h3
-rw-r--r--Eigen/src/Core/util/ForwardDeclarations.h3
-rw-r--r--Eigen/src/Core/util/Macros.h6
-rw-r--r--Eigen/src/LU/LU.h3
-rw-r--r--Eigen/src/QR/QrInstantiations.cpp2
-rw-r--r--Eigen/src/QR/SelfAdjointEigenSolver.h2
-rw-r--r--Eigen/src/SVD/SVD.h2
-rw-r--r--Eigen/src/Sparse/AmbiVector.h12
-rw-r--r--Eigen/src/Sparse/CholmodSupport.h22
-rw-r--r--Eigen/src/Sparse/SparseCholesky.h72
-rw-r--r--Eigen/src/Sparse/TaucsSupport.h27
-rw-r--r--Eigen/src/Sparse/TriangularSolver.h2
-rw-r--r--bench/BenchSparseUtil.h2
-rw-r--r--bench/benchCholesky.cpp32
-rw-r--r--bench/sparse_product.cpp50
-rw-r--r--doc/snippets/LLT_solve.cpp (renamed from doc/snippets/Cholesky_solve.cpp)6
-rw-r--r--test/cholesky.cpp30
-rw-r--r--test/sparse.cpp2
23 files changed, 609 insertions, 184 deletions
diff --git a/Eigen/Cholesky b/Eigen/Cholesky
index f5bcec52d..ddde61a31 100644
--- a/Eigen/Cholesky
+++ b/Eigen/Cholesky
@@ -17,8 +17,8 @@ namespace Eigen {
/** \defgroup Cholesky_Module Cholesky module
* This module provides two variants of the Cholesky decomposition for selfadjoint (hermitian) matrices.
* Those decompositions are accessible via the following MatrixBase methods:
- * - MatrixBase::cholesky(),
- * - MatrixBase::choleskyNoSqrt()
+ * - MatrixBase::llt(),
+ * - MatrixBase::ldlt()
*
* \code
* #include <Eigen/Cholesky>
@@ -27,6 +27,8 @@ namespace Eigen {
#include "src/Array/CwiseOperators.h"
#include "src/Array/Functors.h"
+#include "src/Cholesky/LLT.h"
+#include "src/Cholesky/LDLT.h"
#include "src/Cholesky/Cholesky.h"
#include "src/Cholesky/CholeskyWithoutSquareRoot.h"
diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h
index b94fea8dc..ada413b33 100644
--- a/Eigen/src/Cholesky/Cholesky.h
+++ b/Eigen/src/Cholesky/Cholesky.h
@@ -29,22 +29,7 @@
*
* \class Cholesky
*
- * \brief Standard Cholesky decomposition of a matrix and associated features
- *
- * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
- *
- * This class performs a standard Cholesky decomposition of a symmetric, positive definite
- * matrix A such that A = LL^* = U^*U, where L is lower triangular.
- *
- * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like D^*D x = b,
- * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
- * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
- * situations like generalised eigen problems with hermitian matrices.
- *
- * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
- * the strict lower part does not have to store correct values.
- *
- * \sa MatrixBase::cholesky(), class CholeskyWithoutSquareRoot
+ * \deprecated this class has been renamed LLT
*/
template<typename MatrixType> class Cholesky
{
@@ -72,7 +57,10 @@ template<typename MatrixType> class Cholesky
inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
template<typename Derived>
- typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
+ typename Derived::Eval solve(const MatrixBase<Derived> &b) const EIGEN_DEPRECATED;
+
+ template<typename RhsDerived, typename ResDerived>
+ bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &bAndX) const;
@@ -128,16 +116,7 @@ void Cholesky<MatrixType>::compute(const MatrixType& a)
}
}
-/** \returns the solution of \f$ A x = b \f$ using the current decomposition of A.
- * In other words, it returns \f$ A^{-1} b \f$ computing
- * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
- * \param b the column vector \f$ b \f$, which can also be a matrix.
- *
- * Example: \include Cholesky_solve.cpp
- * Output: \verbinclude Cholesky_solve.out
- *
- * \sa MatrixBase::cholesky(), CholeskyWithoutSquareRoot::solve()
- */
+/** \deprecated */
template<typename MatrixType>
template<typename Derived>
typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b) const
@@ -147,7 +126,6 @@ typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b)
typename ei_eval_to_column_major<Derived>::type x(b);
solveInPlace(x);
return x;
- //return m_matrix.adjoint().template part<Upper>().solveTriangular(matrixL().solveTriangular(b));
}
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
@@ -162,7 +140,25 @@ typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b)
* Example: \include Cholesky_solve.cpp
* Output: \verbinclude Cholesky_solve.out
*
- * \sa MatrixBase::cholesky(), Cholesky::solve()
+ * \sa MatrixBase::cholesky(), Cholesky::solveInPlace()
+ */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool Cholesky<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+ const int size = m_matrix.rows();
+ ei_assert(size==b.rows() && "Cholesky::solve(): invalid number of rows of the right hand side matrix b");
+ return solveInPlace((*result) = b);
+}
+
+/** This is the \em in-place version of solve().
+ *
+ * \param bAndX represents both the right-hand side matrix b and result x.
+ *
+ * This version avoids a copy when the right hand side matrix b is not
+ * needed anymore.
+ *
+ * \sa Cholesky::solve(), MatrixBase::cholesky()
*/
template<typename MatrixType>
template<typename Derived>
@@ -178,7 +174,7 @@ bool Cholesky<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
}
/** \cholesky_module
- * \returns the Cholesky decomposition of \c *this
+ * \deprecated has been renamed llt()
*/
template<typename Derived>
inline const Cholesky<typename MatrixBase<Derived>::EvalType>
diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
index fd111fb1f..af44634a0 100644
--- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
+++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
@@ -25,25 +25,11 @@
#ifndef EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
#define EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
-/** \ingroup Cholesky_Module
+/** \deprecated \ingroup Cholesky_Module
*
* \class CholeskyWithoutSquareRoot
*
- * \brief Robust Cholesky decomposition of a matrix and associated features
- *
- * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
- *
- * This class performs a Cholesky decomposition without square root of a symmetric, positive definite
- * matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal
- * and D is a diagonal matrix.
- *
- * Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more
- * stable computation.
- *
- * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
- * the strict lower part does not have to store correct values.
- *
- * \sa MatrixBase::choleskyNoSqrt(), class Cholesky
+ * \deprecated this class has been renamed LDLT
*/
template<typename MatrixType> class CholeskyWithoutSquareRoot
{
@@ -69,7 +55,10 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
template<typename Derived>
- typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
+ typename Derived::Eval solve(const MatrixBase<Derived> &b) const EIGEN_DEPRECATED;
+
+ template<typename RhsDerived, typename ResDerived>
+ bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &bAndX) const;
@@ -141,15 +130,7 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
}
}
-/** \returns the solution of \f$ A x = b \f$ using the current decomposition of A.
- * In other words, it returns \f$ A^{-1} b \f$ computing
- * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
- * \param b the column vector \f$ b \f$, which can also be a matrix.
- *
- * See Cholesky::solve() for a example.
- *
- * \sa CholeskyWithoutSquareRoot::solveInPlace(), MatrixBase::choleskyNoSqrt()
- */
+/** \deprecated */
template<typename MatrixType>
template<typename Derived>
typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(const MatrixBase<Derived> &b) const
@@ -173,10 +154,30 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(const Matrix
* \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
* \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
*
- * Example: \include Cholesky_solve.cpp
- * Output: \verbinclude Cholesky_solve.out
+ * Example: \include CholeskyCholeskyWithoutSquareRoot_solve.cpp
+ * Output: \verbinclude CholeskyCholeskyWithoutSquareRoot_solve.out
+ *
+ * \sa CholeskyWithoutSquareRoot::solveInPlace(), MatrixBase::choleskyNoSqrt()
+ */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool CholeskyWithoutSquareRoot<MatrixType>
+::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+ const int size = m_matrix.rows();
+ ei_assert(size==b.rows() && "Cholesky::solve(): invalid number of rows of the right hand side matrix b");
+ *result = b;
+ return solveInPlace(*result);
+}
+
+/** This is the \em in-place version of solve().
+ *
+ * \param bAndX represents both the right-hand side matrix b and result x.
+ *
+ * This version avoids a copy when the right hand side matrix b is not
+ * needed anymore.
*
- * \sa MatrixBase::cholesky(), CholeskyWithoutSquareRoot::solve()
+ * \sa CholeskyWithoutSquareRoot::solve(), MatrixBase::choleskyNoSqrt()
*/
template<typename MatrixType>
template<typename Derived>
@@ -187,13 +188,13 @@ bool CholeskyWithoutSquareRoot<MatrixType>::solveInPlace(MatrixBase<Derived> &bA
if (!m_isPositiveDefinite)
return false;
matrixL().solveTriangularInPlace(bAndX);
- bAndX *= m_matrix.cwise().inverse().template part<Diagonal>();
+ bAndX = (m_matrix.cwise().inverse().template part<Diagonal>() * bAndX).lazy();
m_matrix.adjoint().template part<UnitUpper>().solveTriangularInPlace(bAndX);
return true;
}
-/** \cholesky_module
- * \returns the Cholesky decomposition without square root of \c *this
+/** \deprecated \cholesky_module
+ * \deprecated has been renamed ldlt()
*/
template<typename Derived>
inline const CholeskyWithoutSquareRoot<typename MatrixBase<Derived>::EvalType>
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
new file mode 100644
index 000000000..e70a324f6
--- /dev/null
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -0,0 +1,202 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_LDLT_H
+#define EIGEN_LDLT_H
+
+/** \ingroup cholesky_Module
+ *
+ * \class LDLT
+ *
+ * \brief Robust Cholesky decomposition of a matrix and associated features
+ *
+ * \param MatrixType the type of the matrix of which we are computing the LDL^T Cholesky decomposition
+ *
+ * This class performs a Cholesky decomposition without square root of a symmetric, positive definite
+ * matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal
+ * and D is a diagonal matrix.
+ *
+ * Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more
+ * stable computation.
+ *
+ * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+ * the strict lower part does not have to store correct values.
+ *
+ * \sa MatrixBase::ldlt(), class LLT
+ */
+template<typename MatrixType> class LDLT
+{
+ public:
+
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+ typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
+
+ LDLT(const MatrixType& matrix)
+ : m_matrix(matrix.rows(), matrix.cols())
+ {
+ compute(matrix);
+ }
+
+ /** \returns the lower triangular matrix L */
+ inline Part<MatrixType, UnitLower> matrixL(void) const { return m_matrix; }
+
+ /** \returns the coefficients of the diagonal matrix D */
+ inline DiagonalCoeffs<MatrixType> vectorD(void) const { return m_matrix.diagonal(); }
+
+ /** \returns true if the matrix is positive definite */
+ inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
+
+ template<typename RhsDerived, typename ResDerived>
+ bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
+
+ template<typename Derived>
+ bool solveInPlace(MatrixBase<Derived> &bAndX) const;
+
+ void compute(const MatrixType& matrix);
+
+ protected:
+ /** \internal
+ * Used to compute and store the cholesky decomposition A = L D L^* = U^* D U.
+ * The strict upper part is used during the decomposition, the strict lower
+ * part correspond to the coefficients of L (its diagonal is equal to 1 and
+ * is not stored), and the diagonal entries correspond to D.
+ */
+ MatrixType m_matrix;
+
+ bool m_isPositiveDefinite;
+};
+
+/** Compute / recompute the LLT decomposition A = L D L^* = U^* D U of \a matrix
+ */
+template<typename MatrixType>
+void LDLT<MatrixType>::compute(const MatrixType& a)
+{
+ assert(a.rows()==a.cols());
+ const int size = a.rows();
+ m_matrix.resize(size, size);
+ m_isPositiveDefinite = true;
+ const RealScalar eps = ei_sqrt(precision<Scalar>());
+
+ if (size<=1)
+ {
+ m_matrix = a;
+ return;
+ }
+
+ // Let's preallocate a temporay vector to evaluate the matrix-vector product into it.
+ // Unlike the standard LLT decomposition, here we cannot evaluate it to the destination
+ // matrix because it a sub-row which is not compatible suitable for efficient packet evaluation.
+ // (at least if we assume the matrix is col-major)
+ Matrix<Scalar,MatrixType::RowsAtCompileTime,1> _temporary(size);
+
+ // Note that, in this algorithm the rows of the strict upper part of m_matrix is used to store
+ // column vector, thus the strange .conjugate() and .transpose()...
+
+ m_matrix.row(0) = a.row(0).conjugate();
+ m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
+ for (int j = 1; j < size; ++j)
+ {
+ RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
+ m_matrix.coeffRef(j,j) = tmp;
+
+ if (tmp < eps)
+ {
+ m_isPositiveDefinite = false;
+ return;
+ }
+
+ int endSize = size-j-1;
+ if (endSize>0)
+ {
+ _temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
+ * m_matrix.col(j).start(j).conjugate() ).lazy();
+
+ m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
+ - _temporary.end(endSize).transpose();
+
+ m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
+ }
+ }
+}
+
+/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
+ * The result is stored in \a bAndx
+ *
+ * \returns true in case of success, false otherwise.
+ *
+ * In other words, it computes \f$ b = A^{-1} b \f$ with
+ * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
+ * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
+ *
+ * Example: \include LLTLDLT_solve.cpp
+ * Output: \verbinclude LLTLDLT_solve.out
+ *
+ * \sa LDLT::solveInPlace(), MatrixBase::ldlt()
+ */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool LDLT<MatrixType>
+::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+ const int size = m_matrix.rows();
+ ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
+ *result = b;
+ return solveInPlace(*result);
+}
+
+/** This is the \em in-place version of solve().
+ *
+ * \param bAndX represents both the right-hand side matrix b and result x.
+ *
+ * This version avoids a copy when the right hand side matrix b is not
+ * needed anymore.
+ *
+ * \sa LDLT::solve(), MatrixBase::ldlt()
+ */
+template<typename MatrixType>
+template<typename Derived>
+bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
+{
+ const int size = m_matrix.rows();
+ ei_assert(size==bAndX.rows());
+ if (!m_isPositiveDefinite)
+ return false;
+ matrixL().solveTriangularInPlace(bAndX);
+ bAndX = (m_matrix.cwise().inverse().template part<Diagonal>() * bAndX).lazy();
+ m_matrix.adjoint().template part<UnitUpper>().solveTriangularInPlace(bAndX);
+ return true;
+}
+
+/** \cholesky_module
+ * \returns the Cholesky decomposition without square root of \c *this
+ */
+template<typename Derived>
+inline const LDLT<typename MatrixBase<Derived>::EvalType>
+MatrixBase<Derived>::ldlt() const
+{
+ return derived();
+}
+
+#endif // EIGEN_LDLT_H
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
new file mode 100644
index 000000000..8d4a1a752
--- /dev/null
+++ b/Eigen/src/Cholesky/LLT.h
@@ -0,0 +1,186 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_LLT_H
+#define EIGEN_LLT_H
+
+/** \ingroup cholesky_Module
+ *
+ * \class LLT
+ *
+ * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
+ *
+ * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
+ *
+ * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
+ * matrix A such that A = LL^* = U^*U, where L is lower triangular.
+ *
+ * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like D^*D x = b,
+ * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
+ * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
+ * situations like generalised eigen problems with hermitian matrices.
+ *
+ * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+ * the strict lower part does not have to store correct values.
+ *
+ * \sa MatrixBase::llt(), class LDLT
+ */
+template<typename MatrixType> class LLT
+{
+ private:
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+ typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
+
+ enum {
+ PacketSize = ei_packet_traits<Scalar>::size,
+ AlignmentMask = int(PacketSize)-1
+ };
+
+ public:
+
+ LLT(const MatrixType& matrix)
+ : m_matrix(matrix.rows(), matrix.cols())
+ {
+ compute(matrix);
+ }
+
+ inline Part<MatrixType, Lower> matrixL(void) const { return m_matrix; }
+
+ /** \returns true if the matrix is positive definite */
+ inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
+
+ template<typename RhsDerived, typename ResDerived>
+ bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
+
+ template<typename Derived>
+ bool solveInPlace(MatrixBase<Derived> &bAndX) const;
+
+ void compute(const MatrixType& matrix);
+
+ protected:
+ /** \internal
+ * Used to compute and store L
+ * The strict upper part is not used and even not initialized.
+ */
+ MatrixType m_matrix;
+ bool m_isPositiveDefinite;
+};
+
+/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
+ */
+template<typename MatrixType>
+void LLT<MatrixType>::compute(const MatrixType& a)
+{
+ assert(a.rows()==a.cols());
+ const int size = a.rows();
+ m_matrix.resize(size, size);
+ const RealScalar eps = ei_sqrt(precision<Scalar>());
+
+ RealScalar x;
+ x = ei_real(a.coeff(0,0));
+ m_isPositiveDefinite = x > eps && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1));
+ m_matrix.coeffRef(0,0) = ei_sqrt(x);
+ m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
+ for (int j = 1; j < size; ++j)
+ {
+ Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).norm2();
+ x = ei_real(tmp);
+ if (x < eps || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
+ {
+ m_isPositiveDefinite = false;
+ return;
+ }
+ m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
+
+ int endSize = size-j-1;
+ if (endSize>0) {
+ // Note that when all matrix columns have good alignment, then the following
+ // product is guaranteed to be optimal with respect to alignment.
+ m_matrix.col(j).end(endSize) =
+ (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy();
+
+ // FIXME could use a.col instead of a.row
+ m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint()
+ - m_matrix.col(j).end(endSize) ) / x;
+ }
+ }
+}
+
+/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
+ * The result is stored in \a bAndx
+ *
+ * \returns true in case of success, false otherwise.
+ *
+ * In other words, it computes \f$ b = A^{-1} b \f$ with
+ * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
+ * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
+ *
+ * Example: \include LLT_solve.cpp
+ * Output: \verbinclude LLT_solve.out
+ *
+ * \sa LLT::solveInPlace(), MatrixBase::llt()
+ */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+ const int size = m_matrix.rows();
+ ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
+ return solveInPlace((*result) = b);
+}
+
+/** This is the \em in-place version of solve().
+ *
+ * \param bAndX represents both the right-hand side matrix b and result x.
+ *
+ * This version avoids a copy when the right hand side matrix b is not
+ * needed anymore.
+ *
+ * \sa LLT::solve(), MatrixBase::llt()
+ */
+template<typename MatrixType>
+template<typename Derived>
+bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
+{
+ const int size = m_matrix.rows();
+ ei_assert(size==bAndX.rows());
+ if (!m_isPositiveDefinite)
+ return false;
+ matrixL().solveTriangularInPlace(bAndX);
+ m_matrix.adjoint().template part<Upper>().solveTriangularInPlace(bAndX);
+ return true;
+}
+
+/** \cholesky_module
+ * \returns the LLT decomposition of \c *this
+ */
+template<typename Derived>
+inline const LLT<typename MatrixBase<Derived>::EvalType>
+MatrixBase<Derived>::llt() const
+{
+ return LLT<typename ei_eval<Derived>::type>(derived());
+}
+
+#endif // EIGEN_LLT_H
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 944d353d8..8b0129060 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -563,6 +563,9 @@ template<typename Derived> class MatrixBase
/////////// Cholesky module ///////////
+ const LLT<EvalType> llt() const;
+ const LDLT<EvalType> ldlt() const;
+ // deprecated:
const Cholesky<EvalType> cholesky() const;
const CholeskyWithoutSquareRoot<EvalType> choleskyNoSqrt() const;
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 55016e05d..98d15b415 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -102,6 +102,9 @@ template<typename ExpressionType, int Direction> class PartialRedux;
template<typename MatrixType> class LU;
template<typename MatrixType> class QR;
template<typename MatrixType> class SVD;
+template<typename MatrixType> class LLT;
+template<typename MatrixType> class LDLT;
+// deprecated:
template<typename MatrixType> class Cholesky;
template<typename MatrixType> class CholeskyWithoutSquareRoot;
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index e3429e146..348f313d2 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -98,6 +98,12 @@ using Eigen::ei_cos;
#endif
#if (defined __GNUC__)
+#define EIGEN_DEPRECATED __attribute__((deprecated))
+#else
+#define EIGEN_DEPRECATED
+#endif
+
+#if (defined __GNUC__)
#define EIGEN_ALIGN_128 __attribute__ ((aligned(16)))
#else
#define EIGEN_ALIGN_128
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index f7ed1b412..554d8bd3c 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -181,7 +181,8 @@ template<typename MatrixType> class LU
* \returns true if any solution exists, false if no solution exists.
*
* \note If there exist more than one solution, this method will arbitrarily choose one.
- * If you need a complete analysis of the space of solutions, take the one solution obtained * by this method and add to it elements of the kernel, as determined by kernel().
+ * If you need a complete analysis of the space of solutions, take the one solution obtained
+ * by this method and add to it elements of the kernel, as determined by kernel().
*
* Example: \include LU_solve.cpp
* Output: \verbinclude LU_solve.out
diff --git a/Eigen/src/QR/QrInstantiations.cpp b/Eigen/src/QR/QrInstantiations.cpp
index 3b2594e34..dacb05d3d 100644
--- a/Eigen/src/QR/QrInstantiations.cpp
+++ b/Eigen/src/QR/QrInstantiations.cpp
@@ -26,8 +26,6 @@
#define EIGEN_EXTERN_INSTANTIATIONS
#endif
#include "../../Core"
-// commented because of -pedantic
-// #include "../../Cholesky"
#undef EIGEN_EXTERN_INSTANTIATIONS
#include "../../QR"
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index 9e929620c..fdb2764c4 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -227,7 +227,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
// Compute the cholesky decomposition of matB = L L'
- Cholesky<MatrixType> cholB(matB);
+ LLT<MatrixType> cholB(matB);
// compute C = inv(L) A inv(L')
MatrixType matC = matA;
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index c3f3bb235..100ca9310 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -505,7 +505,7 @@ SVD<MatrixType>& SVD<MatrixType>::sort()
/** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A.
* The parts of the solution corresponding to zero singular values are ignored.
*
- * \sa MatrixBase::svd(), LU::solve(), Cholesky::solve()
+ * \sa MatrixBase::svd(), LU::solve(), LLT::solve()
*/
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h
index 0c4bd1af1..44697bceb 100644
--- a/Eigen/src/Sparse/AmbiVector.h
+++ b/Eigen/src/Sparse/AmbiVector.h
@@ -28,7 +28,7 @@
/** \internal
* Hybrid sparse/dense vector class designed for intensive read-write operations.
*
- * See BasicSparseCholesky and SparseProduct for usage examples.
+ * See BasicSparseLLT and SparseProduct for usage examples.
*/
template<typename _Scalar> class AmbiVector
{
@@ -269,12 +269,14 @@ class AmbiVector<_Scalar>::Iterator
/** Default constructor
* \param vec the vector on which we iterate
- * \param nonZeroReferenceValue reference value used to prune zero coefficients.
- * In practice, the coefficient are compared to \a nonZeroReferenceValue * precision<Scalar>().
+ * \param epsilon the minimal value used to prune zero coefficients.
+ * In practice, all coefficients having a magnitude smaller than \a epsilon
+ * are skipped.
*/
- Iterator(const AmbiVector& vec, RealScalar nonZeroReferenceValue = RealScalar(0.1)) : m_vector(vec)
+ Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*precision<RealScalar>())
+ : m_vector(vec)
{
- m_epsilon = nonZeroReferenceValue * precision<Scalar>();
+ m_epsilon = epsilon;
m_isDense = m_vector.m_mode==IsDense;
if (m_isDense)
{
diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
index ed031b264..ba2161320 100644
--- a/Eigen/src/Sparse/CholmodSupport.h
+++ b/Eigen/src/Sparse/CholmodSupport.h
@@ -99,10 +99,10 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
}
template<typename MatrixType>
-class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
+class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
{
protected:
- typedef SparseCholesky<MatrixType> Base;
+ typedef SparseLLT<MatrixType> Base;
using Base::Scalar;
using Base::RealScalar;
using Base::MatrixLIsDirty;
@@ -113,14 +113,20 @@ class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
public:
- SparseCholesky(const MatrixType& matrix, int flags = 0)
+ SparseLLT(int flags = 0)
+ : Base(flags), m_cholmodFactor(0)
+ {
+ cholmod_start(&m_cholmod);
+ }
+
+ SparseLLT(const MatrixType& matrix, int flags = 0)
: Base(matrix, flags), m_cholmodFactor(0)
{
cholmod_start(&m_cholmod);
compute(matrix);
}
- ~SparseCholesky()
+ ~SparseLLT()
{
if (m_cholmodFactor)
cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
@@ -140,7 +146,7 @@ class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
};
template<typename MatrixType>
-void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
+void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
{
if (m_cholmodFactor)
{
@@ -169,8 +175,8 @@ void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
}
template<typename MatrixType>
-inline const typename SparseCholesky<MatrixType>::CholMatrixType&
-SparseCholesky<MatrixType,Cholmod>::matrixL() const
+inline const typename SparseLLT<MatrixType>::CholMatrixType&
+SparseLLT<MatrixType,Cholmod>::matrixL() const
{
if (m_status & MatrixLIsDirty)
{
@@ -187,7 +193,7 @@ SparseCholesky<MatrixType,Cholmod>::matrixL() const
template<typename MatrixType>
template<typename Derived>
-void SparseCholesky<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
+void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());
diff --git a/Eigen/src/Sparse/SparseCholesky.h b/Eigen/src/Sparse/SparseCholesky.h
index 839a84bf3..efe005709 100644
--- a/Eigen/src/Sparse/SparseCholesky.h
+++ b/Eigen/src/Sparse/SparseCholesky.h
@@ -38,25 +38,19 @@ enum {
MemoryEfficient = 0x2,
SupernodalMultifrontal = 0x4,
SupernodalLeftLooking = 0x8
-
-
-/*
- CholUseEigen = 0x0, // Eigen's impl is the default
- CholUseTaucs = 0x2,
- CholUseCholmod = 0x4*/
};
/** \ingroup Sparse_Module
*
- * \class SparseCholesky
+ * \class SparseLLT
*
- * \brief Standard Cholesky decomposition of a matrix and associated features
+ * \brief Standard LLT decomposition of a matrix and associated features
*
- * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
+ * \param MatrixType the type of the matrix of which we are computing the LLT decomposition
*
- * \sa class Cholesky, class CholeskyWithoutSquareRoot
+ * \sa class LLT, class LDLT
*/
-template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
+template<typename MatrixType, int Backend = DefaultBackend> class SparseLLT
{
protected:
typedef typename MatrixType::Scalar Scalar;
@@ -70,66 +64,66 @@ template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
public:
- SparseCholesky(const MatrixType& matrix, int flags = 0)
+ SparseLLT(int flags = 0)
+ : m_flags(flags), m_status(0)
+ {
+ m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
+ }
+
+ SparseLLT(const MatrixType& matrix, int flags = 0)
: m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
{
+ m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
compute(matrix);
}
- inline const CholMatrixType& matrixL(void) const { return m_matrix; }
+ void setPrecision(RealScalar v) { m_precision = v; }
+ RealScalar precision() const { return m_precision; }
- /** \returns true if the matrix is positive definite */
- inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
+ void setFlags(int f) { m_flags = f; }
+ int flags() const { return m_flags; }
+
+ void compute(const MatrixType& matrix);
+
+ inline const CholMatrixType& matrixL(void) const { return m_matrix; }
template<typename Derived>
void solveInPlace(MatrixBase<Derived> &b) const;
- void compute(const MatrixType& matrix);
+ /** \returns true if the factorization succeeded */
+ inline bool succeeded(void) const { return m_succeeded; }
protected:
- /** \internal
- * Used to compute and store L
- * The strict upper part is not used and even not initialized.
- */
CholMatrixType m_matrix;
+ RealScalar m_precision;
int m_flags;
mutable int m_status;
- bool m_isPositiveDefinite;
+ bool m_succeeded;
};
-/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
+/** Computes / recomputes the LLT decomposition A = LL^* = U^*U of \a matrix
*/
template<typename MatrixType, int Backend>
-void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
+void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
{
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
- const RealScalar eps = ei_sqrt(precision<Scalar>());
+// const RealScalar eps = ei_sqrt(precision<Scalar>());
// allocate a temporary vector for accumulations
AmbiVector<Scalar> tempVector(size);
+ RealScalar density = a.nonZeros()/RealScalar(size*size);
// TODO estimate the number of nnz
m_matrix.startFill(a.nonZeros()*2);
for (int j = 0; j < size; ++j)
{
-// std::cout << j << "\n";
Scalar x = ei_real(a.coeff(j,j));
int endSize = size-j-1;
- // TODO estimate the number of non zero entries
-// float ratioLhs = float(lhs.nonZeros())/float(lhs.rows()*lhs.cols());
-// float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
-// float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
-
- // let's do a more accurate determination of the nnz ratio for the current column j of res
- //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
- // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
-// float ratioColRes = ratioRes;
-// if (ratioColRes>0.1)
-// tempVector.init(IsSparse);
- tempVector.init(IsDense);
+ // TODO better estimate the density !
+ tempVector.init(density>0.001? IsDense : IsSparse);
tempVector.setBounds(j+1,size);
tempVector.setZero();
// init with current matrix a
@@ -161,7 +155,7 @@ void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
RealScalar rx = ei_sqrt(ei_real(x));
m_matrix.fill(j,j) = rx;
Scalar y = Scalar(1)/rx;
- for (typename AmbiVector<Scalar>::Iterator it(tempVector); it; ++it)
+ for (typename AmbiVector<Scalar>::Iterator it(tempVector, m_precision*rx); it; ++it)
{
m_matrix.fill(it.index(), j) = it.value() * y;
}
@@ -171,7 +165,7 @@ void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
template<typename MatrixType, int Backend>
template<typename Derived>
-void SparseCholesky<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
+void SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());
diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h
index a0a5b4b71..a92a63b57 100644
--- a/Eigen/src/Sparse/TaucsSupport.h
+++ b/Eigen/src/Sparse/TaucsSupport.h
@@ -79,10 +79,10 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(taucs_ccs_matrix& tau
}
template<typename MatrixType>
-class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
+class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
{
protected:
- typedef SparseCholesky<MatrixType> Base;
+ typedef SparseLLT<MatrixType> Base;
using Base::Scalar;
using Base::RealScalar;
using Base::MatrixLIsDirty;
@@ -93,13 +93,18 @@ class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
public:
- SparseCholesky(const MatrixType& matrix, int flags = 0)
+ SparseLLT(int flags = 0)
+ : Base(flags), m_taucsSupernodalFactor(0)
+ {
+ }
+
+ SparseLLT(const MatrixType& matrix, int flags = 0)
: Base(matrix, flags), m_taucsSupernodalFactor(0)
{
compute(matrix);
}
- ~SparseCholesky()
+ ~SparseLLT()
{
if (m_taucsSupernodalFactor)
taucs_supernodal_factor_free(m_taucsSupernodalFactor);
@@ -117,7 +122,7 @@ class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
};
template<typename MatrixType>
-void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
+void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a)
{
if (m_taucsSupernodalFactor)
{
@@ -128,7 +133,7 @@ void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
if (m_flags & IncompleteFactorization)
{
taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
- taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
+ taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0);
m_matrix = Base::CholMatrixType::Map(*taucsRes);
free(taucsRes);
m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
@@ -153,8 +158,8 @@ void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
}
template<typename MatrixType>
-inline const typename SparseCholesky<MatrixType>::CholMatrixType&
-SparseCholesky<MatrixType,Taucs>::matrixL() const
+inline const typename SparseLLT<MatrixType>::CholMatrixType&
+SparseLLT<MatrixType,Taucs>::matrixL() const
{
if (m_status & MatrixLIsDirty)
{
@@ -170,7 +175,7 @@ SparseCholesky<MatrixType,Taucs>::matrixL() const
template<typename MatrixType>
template<typename Derived>
-void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
+void SparseLLT<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());
@@ -179,9 +184,9 @@ void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) cons
{
// ei_assert(!(m_status & SupernodalFactorIsDirty));
// taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b);
- matrixL();
+ //matrixL();
}
-// else
+ else
{
Base::solveInPlace(b);
}
diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h
index 3a1a33b3c..8e71c936d 100644
--- a/Eigen/src/Sparse/TriangularSolver.h
+++ b/Eigen/src/Sparse/TriangularSolver.h
@@ -102,7 +102,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Lower,ColMajor|IsSparse>
typename Lhs::InnerIterator it(lhs, i);
if(!(Lhs::Flags & UnitDiagBit))
{
- std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n";
+ // std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n";
ei_assert(it.index()==i);
other.coeffRef(i,col) /= it.value();
}
diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h
index 2c24c29e6..35c9a5263 100644
--- a/bench/BenchSparseUtil.h
+++ b/bench/BenchSparseUtil.h
@@ -16,7 +16,7 @@ USING_PART_OF_NAMESPACE_EIGEN
#endif
#ifndef SCALAR
-#define SCALAR float
+#define SCALAR double
#endif
typedef SCALAR Scalar;
diff --git a/bench/benchCholesky.cpp b/bench/benchCholesky.cpp
index f64b61b71..e998d8536 100644
--- a/bench/benchCholesky.cpp
+++ b/bench/benchCholesky.cpp
@@ -1,5 +1,5 @@
-// g++ -DNDEBUG -O3 -I.. benchCholesky.cpp -o benchCholesky && ./benchCholesky
+// g++ -DNDEBUG -O3 -I.. benchLLT.cpp -o benchLLT && ./benchLLT
// options:
// -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3
// -DEIGEN_DONT_VECTORIZE
@@ -9,7 +9,7 @@
// -DSCALAR=double
#include <Eigen/Array>
-#include <Eigen/Cholesky>
+#include <Eigen/LLT>
#include <bench/BenchUtil.h>
using namespace Eigen;
@@ -24,7 +24,7 @@ using namespace Eigen;
typedef float Scalar;
template <typename MatrixType>
-__attribute__ ((noinline)) void benchCholesky(const MatrixType& m)
+__attribute__ ((noinline)) void benchLLT(const MatrixType& m)
{
int rows = m.rows();
int cols = m.cols();
@@ -54,7 +54,7 @@ __attribute__ ((noinline)) void benchCholesky(const MatrixType& m)
timerNoSqrt.start();
for (int k=0; k<repeats; ++k)
{
- CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(covMat);
+ LDLT<SquareMatrixType> cholnosqrt(covMat);
acc += cholnosqrt.matrixL().coeff(r,c);
}
timerNoSqrt.stop();
@@ -65,7 +65,7 @@ __attribute__ ((noinline)) void benchCholesky(const MatrixType& m)
timerSqrt.start();
for (int k=0; k<repeats; ++k)
{
- Cholesky<SquareMatrixType> chol(covMat);
+ LLT<SquareMatrixType> chol(covMat);
acc += chol.matrixL().coeff(r,c);
}
timerSqrt.stop();
@@ -124,17 +124,17 @@ int main(int argc, char* argv[])
std::cout << "\n";
for (uint i=0; dynsizes[i]>0; ++i)
- benchCholesky(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
-
-// benchCholesky(Matrix<Scalar,2,2>());
-// benchCholesky(Matrix<Scalar,3,3>());
-// benchCholesky(Matrix<Scalar,4,4>());
-// benchCholesky(Matrix<Scalar,5,5>());
-// benchCholesky(Matrix<Scalar,6,6>());
-// benchCholesky(Matrix<Scalar,7,7>());
-// benchCholesky(Matrix<Scalar,8,8>());
-// benchCholesky(Matrix<Scalar,12,12>());
-// benchCholesky(Matrix<Scalar,16,16>());
+ benchLLT(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
+
+// benchLLT(Matrix<Scalar,2,2>());
+// benchLLT(Matrix<Scalar,3,3>());
+// benchLLT(Matrix<Scalar,4,4>());
+// benchLLT(Matrix<Scalar,5,5>());
+// benchLLT(Matrix<Scalar,6,6>());
+// benchLLT(Matrix<Scalar,7,7>());
+// benchLLT(Matrix<Scalar,8,8>());
+// benchLLT(Matrix<Scalar,12,12>());
+// benchLLT(Matrix<Scalar,16,16>());
return 0;
}
diff --git a/bench/sparse_product.cpp b/bench/sparse_product.cpp
index 846301fa5..edeb08c5d 100644
--- a/bench/sparse_product.cpp
+++ b/bench/sparse_product.cpp
@@ -21,6 +21,18 @@
#define MINDENSITY 0.0004
#endif
+#ifndef NBTRIES
+#define NBTRIES 10
+#endif
+
+#define BENCH(X) \
+ timer.reset(); \
+ for (int _j=0; _j<NBTRIES; ++_j) { \
+ timer.start(); \
+ for (int _k=0; _k<REPEAT; ++_k) { \
+ X \
+ } timer.stop(); }
+
int main(int argc, char *argv[])
{
int rows = SIZE;
@@ -77,32 +89,34 @@ int main(int argc, char *argv[])
{
std::cout << "Eigen sparse\t" << density*100 << "%\n";
- timer.reset();
- timer.start();
- for (int k=0; k<REPEAT; ++k)
- sm3 = sm1 * sm2;
- timer.stop();
+// timer.reset();
+// timer.start();
+ BENCH(for (int k=0; k<REPEAT; ++k) sm3 = sm1 * sm2;)
+// timer.stop();
std::cout << " a * b:\t" << timer.value() << endl;
timer.reset();
timer.start();
- for (int k=0; k<REPEAT; ++k)
- sm3 = sm1.transpose() * sm2;
- timer.stop();
+// std::cerr << "transpose...\n";
+// EigenSparseMatrix sm4 = sm1.transpose();
+// std::cout << sm4.nonZeros() << " == " << sm1.nonZeros() << "\n";
+// exit(1);
+// std::cerr << "transpose OK\n";
+// std::cout << sm1 << "\n\n" << sm1.transpose() << "\n\n" << sm4.transpose() << "\n\n";
+ BENCH(for (int k=0; k<REPEAT; ++k) sm3 = sm1.transpose() * sm2;)
+// timer.stop();
std::cout << " a' * b:\t" << timer.value() << endl;
- timer.reset();
- timer.start();
- for (int k=0; k<REPEAT; ++k)
- sm3 = sm1.transpose() * sm2.transpose();
- timer.stop();
+// timer.reset();
+// timer.start();
+ BENCH( for (int k=0; k<REPEAT; ++k) sm3 = sm1.transpose() * sm2.transpose(); )
+// timer.stop();
std::cout << " a' * b':\t" << timer.value() << endl;
- timer.reset();
- timer.start();
- for (int k=0; k<REPEAT; ++k)
- sm3 = sm1 * sm2.transpose();
- timer.stop();
+// timer.reset();
+// timer.start();
+ BENCH( for (int k=0; k<REPEAT; ++k) sm3 = sm1 * sm2.transpose(); )
+// timer.stop();
std::cout << " a * b' :\t" << timer.value() << endl;
}
diff --git a/doc/snippets/Cholesky_solve.cpp b/doc/snippets/LLT_solve.cpp
index ac743cb55..76ab09ec5 100644
--- a/doc/snippets/Cholesky_solve.cpp
+++ b/doc/snippets/LLT_solve.cpp
@@ -2,5 +2,7 @@ typedef Matrix<float,Dynamic,2> DataMatrix;
// let's generate some samples on the 3D plane of equation z = 2x+3y (with some noise)
DataMatrix samples = DataMatrix::Random(12,2);
VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1;
-// and let's solve samples * x = elevations in least square sense:
-cout << (samples.adjoint() * samples).cholesky().solve((samples.adjoint()*elevations).eval()) << endl;
+// and let's solve samples * [x y]^T = elevations in least square sense:
+Matrix<float,2,1> xy;
+(samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations), &xy);
+cout << xy << endl;
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 80614346c..f7eb7800e 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -33,7 +33,7 @@
template<typename MatrixType> void cholesky(const MatrixType& m)
{
/* this test covers the following files:
- Cholesky.h CholeskyWithoutSquareRoot.h
+ LLT.h LDLT.h
*/
int rows = m.rows();
int cols = m.cols();
@@ -44,8 +44,8 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
MatrixType a0 = MatrixType::Random(rows,cols);
- VectorType vecB = VectorType::Random(rows);
- MatrixType matB = MatrixType::Random(rows,cols);
+ VectorType vecB = VectorType::Random(rows), vecX(rows);
+ MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
SquareMatrixType symm = a0 * a0.adjoint();
// let's make sure the matrix is not singular or near singular
MatrixType a1 = MatrixType::Random(rows,cols);
@@ -80,28 +80,32 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
#endif
{
- CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(symm);
- VERIFY(cholnosqrt.isPositiveDefinite());
- VERIFY_IS_APPROX(symm, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint());
- VERIFY_IS_APPROX(symm * cholnosqrt.solve(vecB), vecB);
- VERIFY_IS_APPROX(symm * cholnosqrt.solve(matB), matB);
+ LDLT<SquareMatrixType> ldlt(symm);
+ VERIFY(ldlt.isPositiveDefinite());
+ VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint());
+ ldlt.solve(vecB, &vecX);
+ VERIFY_IS_APPROX(symm * vecX, vecB);
+ ldlt.solve(matB, &matX);
+ VERIFY_IS_APPROX(symm * matX, matB);
}
{
- Cholesky<SquareMatrixType> chol(symm);
+ LLT<SquareMatrixType> chol(symm);
VERIFY(chol.isPositiveDefinite());
VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint());
- VERIFY_IS_APPROX(symm * chol.solve(vecB), vecB);
- VERIFY_IS_APPROX(symm * chol.solve(matB), matB);
+ chol.solve(vecB, &vecX);
+ VERIFY_IS_APPROX(symm * vecX, vecB);
+ chol.solve(matB, &matX);
+ VERIFY_IS_APPROX(symm * matX, matB);
}
// test isPositiveDefinite on non definite matrix
if (rows>4)
{
SquareMatrixType symm = a0.block(0,0,rows,cols-4) * a0.block(0,0,rows,cols-4).adjoint();
- Cholesky<SquareMatrixType> chol(symm);
+ LLT<SquareMatrixType> chol(symm);
VERIFY(!chol.isPositiveDefinite());
- CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(symm);
+ LDLT<SquareMatrixType> cholnosqrt(symm);
VERIFY(!cholnosqrt.isPositiveDefinite());
}
}
diff --git a/test/sparse.cpp b/test/sparse.cpp
index 39ea05b8b..f54835972 100644
--- a/test/sparse.cpp
+++ b/test/sparse.cpp
@@ -217,7 +217,7 @@ template<typename Scalar> void sparse(int rows, int cols)
// TODO test row major
}
- // test Cholesky
+ // test LLT
{
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-24 11:46:59 +0000