GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-10 02:39:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

remove dead code

Browse Source
This commit is contained in:
Gael Guennebaud 2011-12-08 23:22:28 +01:00
parent e36a4c880a
commit 86bb20c431
8 changed files with 0 additions and 1408 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
unsupported/Eigen/SparseExtra
View File

@ -29,29 +29,6 @@ namespace Eigen {
* \endcode
*/
struct DefaultBackend {};
// solver flags
enum {
CompleteFactorization = 0x0000, // the default
IncompleteFactorization = 0x0001,
MemoryEfficient = 0x0002,
// For LLT Cholesky:
SupernodalMultifrontal = 0x0010,
SupernodalLeftLooking = 0x0020,
// Ordering methods:
NaturalOrdering = 0x0100, // the default
MinimumDegree_AT_PLUS_A = 0x0200,
MinimumDegree_ATA = 0x0300,
ColApproxMinimumDegree = 0x0400,
Metis = 0x0500,
Scotch = 0x0600,
Chaco = 0x0700,
OrderingMask = 0x0f00
};
#include "../../Eigen/src/misc/Solve.h"
#include "../../Eigen/src/misc/SparseSolve.h"
@ -62,10 +39,6 @@ enum {
#include "src/SparseExtra/MarketIO.h"
#include "src/SparseExtra/SparseLLT.h"
#include "src/SparseExtra/SparseLDLTLegacy.h"
#include "src/SparseExtra/SparseLU.h"
} // namespace Eigen
#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"

416
unsupported/Eigen/src/SparseExtra/SparseLDLTLegacy.h
View File

@ -1,416 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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/>.
/*
NOTE: the _symbolic, and _numeric functions has been adapted from
the LDL library:
LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved.
LDL License:
Your use or distribution of LDL or any modified version of
LDL implies that you agree to this License.
This library 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 2.1 of the License, or (at your option) any later version.
This library 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 for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
USA
Permission is hereby granted to use or copy this program under the
terms of the GNU LGPL, provided that the Copyright, this License,
and the Availability of the original version is retained on all copies.
User documentation of any code that uses this code or any modified
version of this code must cite the Copyright, this License, the
Availability note, and "Used by permission." Permission to modify
the code and to distribute modified code is granted, provided the
Copyright, this License, and the Availability note are retained,
and a notice that the code was modified is included.
*/
#ifndef EIGEN_SPARSELDLT_LEGACY_H
#define EIGEN_SPARSELDLT_LEGACY_H
/** \deprecated use class SimplicialLDLT, or class SimplicialLLT, class ConjugateGradient
* \ingroup Sparse_Module
*
* \class SparseLDLT
*
* \brief LDLT Cholesky decomposition of a sparse matrix and associated features
*
* \param MatrixType the type of the matrix of which we are computing the LDLT Cholesky decomposition
*
* \warning the upper triangular part has to be specified. The rest of the matrix is not used. The input matrix must be column major.
*
* \sa class LDLT, class LDLT
*/
template<typename _MatrixType, typename Backend = DefaultBackend>
class SparseLDLT
{
protected:
typedef typename _MatrixType::Scalar Scalar;
typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
typedef Matrix<Scalar,_MatrixType::ColsAtCompileTime,1> VectorType;
enum {
SupernodalFactorIsDirty = 0x10000,
MatrixLIsDirty = 0x20000
};
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Index Index;
typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
/** \deprecated the entire class is deprecated
* Creates a dummy LDLT factorization object with flags \a flags. */
EIGEN_DEPRECATED SparseLDLT(int flags = 0)
: m_flags(flags), m_status(0)
{
eigen_assert((MatrixType::Flags&RowMajorBit)==0);
m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
}
/** \deprecated the entire class is deprecated
* Creates a LDLT object and compute the respective factorization of \a matrix using
* flags \a flags. */
EIGEN_DEPRECATED SparseLDLT(const MatrixType& matrix, int flags = 0)
: m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
{
eigen_assert((MatrixType::Flags&RowMajorBit)==0);
m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
compute(matrix);
}
/** Sets the relative threshold value used to prune zero coefficients during the decomposition.
*
* Setting a value greater than zero speeds up computation, and yields to an imcomplete
* factorization with fewer non zero coefficients. Such approximate factors are especially
* useful to initialize an iterative solver.
*
* \warning if precision is greater that zero, the LDLT factorization is not guaranteed to succeed
* even if the matrix is positive definite.
*
* Note that the exact meaning of this parameter might depends on the actual
* backend. Moreover, not all backends support this feature.
*
* \sa precision() */
void setPrecision(RealScalar v) { m_precision = v; }
/** \returns the current precision.
*
* \sa setPrecision() */
RealScalar precision() const { return m_precision; }
/** Sets the flags. Possible values are:
* - CompleteFactorization
* - IncompleteFactorization
* - MemoryEfficient (hint to use the memory most efficient method offered by the backend)
* - SupernodalMultifrontal (implies a complete factorization if supported by the backend,
* overloads the MemoryEfficient flags)
* - SupernodalLeftLooking (implies a complete factorization if supported by the backend,
* overloads the MemoryEfficient flags)
*
* \sa flags() */
void settags(int f) { m_flags = f; }
/** \returns the current flags */
int flags() const { return m_flags; }
/** Computes/re-computes the LDLT factorization */
void compute(const MatrixType& matrix);
/** Perform a symbolic factorization */
void _symbolic(const MatrixType& matrix);
/** Perform the actual factorization using the previously
* computed symbolic factorization */
bool _numeric(const MatrixType& matrix);
/** \returns the lower triangular matrix L */
inline const CholMatrixType& matrixL(void) const { return m_matrix; }
/** \returns the coefficients of the diagonal matrix D */
inline VectorType vectorD(void) const { return m_diag; }
template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &b) const;
template<typename Rhs>
inline const internal::solve_retval<SparseLDLT<MatrixType>, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(true && "SparseLDLT is not initialized.");
return internal::solve_retval<SparseLDLT<MatrixType>, Rhs>(*this, b.derived());
}
inline Index cols() const { return m_matrix.cols(); }
inline Index rows() const { return m_matrix.rows(); }
inline const VectorType& diag() const { return m_diag; }
/** \returns true if the factorization succeeded */
inline bool succeeded(void) const { return m_succeeded; }
protected:
CholMatrixType m_matrix;
VectorType m_diag;
VectorXi m_parent; // elimination tree
VectorXi m_nonZerosPerCol;
// VectorXi m_w; // workspace
PermutationMatrix<Dynamic,Dynamic,Index> m_P;
PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;
RealScalar m_precision;
int m_flags;
mutable int m_status;
bool m_succeeded;
};
namespace internal {
template<typename _MatrixType, typename Rhs>
struct solve_retval<SparseLDLT<_MatrixType>, Rhs>
: solve_retval_base<SparseLDLT<_MatrixType>, Rhs>
{
typedef SparseLDLT<_MatrixType> SpLDLTDecType;
EIGEN_MAKE_SOLVE_HELPERS(SpLDLTDecType,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
//Index size = dec().matrixL().rows();
eigen_assert(dec().matrixL().rows()==rhs().rows());
Rhs b(rhs().rows(), rhs().cols());
b = rhs();
if (dec().matrixL().nonZeros()>0) // otherwise L==I
dec().matrixL().template triangularView<UnitLower>().solveInPlace(b);
b = b.cwiseQuotient(dec().diag());
if (dec().matrixL().nonZeros()>0) // otherwise L==I
dec().matrixL().adjoint().template triangularView<UnitUpper>().solveInPlace(b);
dst = b;
}
};
} // end namespace internal
/** Computes / recomputes the LDLT decomposition of matrix \a a
* using the default algorithm.
*/
template<typename _MatrixType, typename Backend>
void SparseLDLT<_MatrixType,Backend>::compute(const _MatrixType& a)
{
_symbolic(a);
m_succeeded = _numeric(a);
}
template<typename _MatrixType, typename Backend>
void SparseLDLT<_MatrixType,Backend>::_symbolic(const _MatrixType& a)
{
assert(a.rows()==a.cols());
const Index size = a.rows();
m_matrix.resize(size, size);
m_parent.resize(size);
m_nonZerosPerCol.resize(size);
ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
const Index* Ap = a.outerIndexPtr();
const Index* Ai = a.innerIndexPtr();
Index* Lp = m_matrix.outerIndexPtr();
const Index* P = 0;
Index* Pinv = 0;
if(P)
{
m_P.indices() = Map<const Matrix<Index,Dynamic,1> >(P,size);
m_Pinv = m_P.inverse();
Pinv = m_Pinv.indices().data();
}
else
{
m_P.resize(0);
m_Pinv.resize(0);
}
for (Index 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 */
tags[k] = k; /* mark node k as visited */
m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
Index kk = P ? P[k] : k; /* kth original, or permuted, column */
Index p2 = Ap[kk+1];
for (Index p = Ap[kk]; p < p2; ++p)
{
/* A (i,k) is nonzero (original or permuted A) */
Index i = Pinv ? Pinv[Ai[p]] : Ai[p];
if (i < k)
{
/* follow path from i to root of etree, stop at flagged node */
for (; tags[i] != k; i = m_parent[i])
{
/* find parent of i if not yet determined */
if (m_parent[i] == -1)
m_parent[i] = k;
++m_nonZerosPerCol[i]; /* L (k,i) is nonzero */
tags[i] = k; /* mark i as visited */
}
}
}
}
/* construct Lp index array from m_nonZerosPerCol column counts */
Lp[0] = 0;
for (Index k = 0; k < size; ++k)
Lp[k+1] = Lp[k] + m_nonZerosPerCol[k];
m_matrix.resizeNonZeros(Lp[size]);
}
template<typename _MatrixType, typename Backend>
bool SparseLDLT<_MatrixType,Backend>::_numeric(const _MatrixType& a)
{
assert(a.rows()==a.cols());
const Index size = a.rows();
assert(m_parent.size()==size);
assert(m_nonZerosPerCol.size()==size);
const Index* Ap = a.outerIndexPtr();
const Index* Ai = a.innerIndexPtr();
const Scalar* Ax = a.valuePtr();
const Index* Lp = m_matrix.outerIndexPtr();
Index* Li = m_matrix.innerIndexPtr();
Scalar* Lx = m_matrix.valuePtr();
m_diag.resize(size);
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);
Index* P = 0;
Index* Pinv = 0;
if(m_P.size()==size)
{
P = m_P.indices().data();
Pinv = m_Pinv.indices().data();
}
bool ok = true;
for (Index 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 */
tags[k] = k; /* mark node k as visited */
m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
Index kk = (P) ? (P[k]) : (k); /* kth original, or permuted, column */
Index p2 = Ap[kk+1];
for (Index p = Ap[kk]; p < p2; ++p)
{
Index i = Pinv ? Pinv[Ai[p]] : Ai[p]; /* get A(i,k) */
if (i <= k)
{
y[i] += internal::conj(Ax[p]); /* scatter A(i,k) into Y (sum duplicates) */
Index len;
for (len = 0; tags[i] != k; i = m_parent[i])
{
pattern[len++] = i; /* L(k,i) is nonzero */
tags[i] = k; /* mark i as visited */
}
while (len > 0)
pattern[--top] = pattern[--len];
}
}
/* compute numerical values kth row of L (a sparse triangular solve) */
m_diag[k] = y[k]; /* get D(k,k) and clear Y(k) */
y[k] = 0.0;
for (; top < size; ++top)
{
Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
Scalar yi = (y[i]); /* get and clear Y(i) */
y[i] = 0.0;
Index p2 = Lp[i] + m_nonZerosPerCol[i];
Index p;
for (p = Lp[i]; p < p2; ++p)
y[Li[p]] -= internal::conj(Lx[p]) * (yi);
Scalar l_ki = yi / m_diag[i]; /* the nonzero entry L(k,i) */
m_diag[k] -= l_ki * internal::conj(yi);
Li[p] = k; /* store L(k,i) in column form of L */
Lx[p] = (l_ki);
++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
}
if (m_diag[k] == 0.0)
{
ok = false; /* failure, D(k,k) is zero */
break;
}
}
return ok; /* success, diagonal of D is all nonzero */
}
/** Computes b = L^-T D^-1 L^-1 b */
template<typename _MatrixType, typename Backend>
template<typename Derived>
bool SparseLDLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
{
//Index size = m_matrix.rows();
eigen_assert(m_matrix.rows()==b.rows());
if (!m_succeeded)
return false;
if(m_P.size()>0)
b = m_Pinv * b;
if (m_matrix.nonZeros()>0) // otherwise L==I
m_matrix.template triangularView<UnitLower>().solveInPlace(b);
b = b.cwiseQuotient(m_diag);
if (m_matrix.nonZeros()>0) // otherwise L==I
m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(b);
if(m_P.size()>0)
b = m_P * b;
return true;
}
#endif // EIGEN_SPARSELDLT_LEGACY_H

248
unsupported/Eigen/src/SparseExtra/SparseLLT.h
View File

@ -1,248 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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_SPARSELLT_H
#define EIGEN_SPARSELLT_H
/** \deprecated use class SimplicialLDLT, or class SimplicialLLT, class ConjugateGradient
* \ingroup Sparse_Module
*
* \class SparseLLT
*
* \brief LLT Cholesky decomposition of a sparse matrix and associated features
*
* \param MatrixType the type of the matrix of which we are computing the LLT Cholesky decomposition
*
* \sa class LLT, class LDLT
*/
template<typename _MatrixType, typename Backend = DefaultBackend>
class SparseLLT
{
protected:
typedef typename _MatrixType::Scalar Scalar;
typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
enum {
SupernodalFactorIsDirty = 0x10000,
MatrixLIsDirty = 0x20000
};
public:
typedef SparseMatrix<Scalar> CholMatrixType;
typedef _MatrixType MatrixType;
typedef typename MatrixType::Index Index;
/** \deprecated the entire class is deprecated
* Creates a dummy LLT factorization object with flags \a flags. */
EIGEN_DEPRECATED SparseLLT(int flags = 0)
: m_flags(flags), m_status(0)
{
m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
}
/** \deprecated the entire class is deprecated
* Creates a LLT object and compute the respective factorization of \a matrix using
* flags \a flags. */
EIGEN_DEPRECATED 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::NumTraits<RealScalar>::dummy_precision();
compute(matrix);
}
/** Sets the relative threshold value used to prune zero coefficients during the decomposition.
*
* Setting a value greater than zero speeds up computation, and yields to an imcomplete
* factorization with fewer non zero coefficients. Such approximate factors are especially
* useful to initialize an iterative solver.
*
* \warning if precision is greater that zero, the LLT factorization is not guaranteed to succeed
* even if the matrix is positive definite.
*
* Note that the exact meaning of this parameter might depends on the actual
* backend. Moreover, not all backends support this feature.
*
* \sa precision() */
void setPrecision(RealScalar v) { m_precision = v; }
/** \returns the current precision.
*
* \sa setPrecision() */
RealScalar precision() const { return m_precision; }
/** Sets the flags. Possible values are:
* - CompleteFactorization
* - IncompleteFactorization
* - MemoryEfficient (hint to use the memory most efficient method offered by the backend)
* - SupernodalMultifrontal (implies a complete factorization if supported by the backend,
* overloads the MemoryEfficient flags)
* - SupernodalLeftLooking (implies a complete factorization if supported by the backend,
* overloads the MemoryEfficient flags)
*
* \sa flags() */
void setFlags(int f) { m_flags = f; }
/** \returns the current flags */
int flags() const { return m_flags; }
/** Computes/re-computes the LLT factorization */
void compute(const MatrixType& matrix);
/** \returns the lower triangular matrix L */
inline const CholMatrixType& matrixL(void) const { return m_matrix; }
template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &b) const;
template<typename Rhs>
inline const internal::solve_retval<SparseLLT<MatrixType>, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(true && "SparseLLT is not initialized.");
return internal::solve_retval<SparseLLT<MatrixType>, Rhs>(*this, b.derived());
}
inline Index cols() const { return m_matrix.cols(); }
inline Index rows() const { return m_matrix.rows(); }
/** \returns true if the factorization succeeded */
inline bool succeeded(void) const { return m_succeeded; }
protected:
CholMatrixType m_matrix;
RealScalar m_precision;
int m_flags;
mutable int m_status;
bool m_succeeded;
};
namespace internal {
template<typename _MatrixType, typename Rhs>
struct solve_retval<SparseLLT<_MatrixType>, Rhs>
: solve_retval_base<SparseLLT<_MatrixType>, Rhs>
{
typedef SparseLLT<_MatrixType> SpLLTDecType;
EIGEN_MAKE_SOLVE_HELPERS(SpLLTDecType,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
const Index size = dec().matrixL().rows();
eigen_assert(size==rhs().rows());
Rhs b(rhs().rows(), rhs().cols());
b = rhs();
dec().matrixL().template triangularView<Lower>().solveInPlace(b);
dec().matrixL().adjoint().template triangularView<Upper>().solveInPlace(b);
dst = b;
}
};
} // end namespace internal
/** Computes / recomputes the LLT decomposition of matrix \a a
* using the default algorithm.
*/
template<typename _MatrixType, typename Backend>
void SparseLLT<_MatrixType,Backend>::compute(const _MatrixType& a)
{
assert(a.rows()==a.cols());
const Index size = a.rows();
m_matrix.resize(size, size);
// allocate a temporary vector for accumulations
internal::AmbiVector<Scalar,Index> tempVector(size);
RealScalar density = a.nonZeros()/RealScalar(size*size);
// TODO estimate the number of non zeros
m_matrix.setZero();
m_matrix.reserve(a.nonZeros()*10);
for (Index j = 0; j < size; ++j)
{
Scalar x = internal::real(a.coeff(j,j));
// TODO better estimate of the density !
tempVector.init(density>0.001? IsDense : IsSparse);
tempVector.setBounds(j+1,size);
tempVector.setZero();
// init with current matrix a
{
typename _MatrixType::InnerIterator it(a,j);
eigen_assert(it.index()==j &&
"matrix must has non zero diagonal entries and only the lower triangular part must be stored");
++it; // skip diagonal element
for (; it; ++it)
tempVector.coeffRef(it.index()) = it.value();
}
for (Index k=0; k<j+1; ++k)
{
typename CholMatrixType::InnerIterator it(m_matrix, k);
while (it && it.index()<j)
++it;
if (it && it.index()==j)
{
Scalar y = it.value();
x -= internal::abs2(y);
++it; // skip j-th element, and process remaining column coefficients
tempVector.restart();
for (; it; ++it)
{
tempVector.coeffRef(it.index()) -= it.value() * y;
}
}
}
// copy the temporary vector to the respective m_matrix.col()
// while scaling the result by 1/real(x)
RealScalar rx = internal::sqrt(internal::real(x));
m_matrix.insert(j,j) = rx; // FIXME use insertBack
Scalar y = Scalar(1)/rx;
for (typename internal::AmbiVector<Scalar,Index>::Iterator it(tempVector, m_precision*rx); it; ++it)
{
// FIXME use insertBack
m_matrix.insertBack(it.index(), j) = it.value() * y;
}
}
m_matrix.finalize();
}
/** Computes b = L^-T L^-1 b */
template<typename _MatrixType, typename Backend>
template<typename Derived>
bool SparseLLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
{
const Index size = m_matrix.rows();
eigen_assert(size==b.rows());
m_matrix.template triangularView<Lower>().solveInPlace(b);
m_matrix.adjoint().template triangularView<Upper>().solveInPlace(b);
return true;
}
#endif // EIGEN_SPARSELLT_H

166
unsupported/Eigen/src/SparseExtra/SparseLU.h
View File

@ -1,166 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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_SPARSELU_H
#define EIGEN_SPARSELU_H
enum {
SvNoTrans = 0,
SvTranspose = 1,
SvAdjoint = 2
};
/** \deprecated use class BiCGSTAB, class SuperLU, or class UmfPackLU
* \ingroup Sparse_Module
*
* \class SparseLU
*
* \brief LU decomposition of a sparse matrix and associated features
*
* \param _MatrixType the type of the matrix of which we are computing the LU factorization
*
* \sa class FullPivLU, class SparseLLT
*/
template<typename _MatrixType, typename Backend = DefaultBackend>
class SparseLU
{
protected:
typedef typename _MatrixType::Scalar Scalar;
typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
typedef SparseMatrix<Scalar> LUMatrixType;
enum {
MatrixLUIsDirty = 0x10000
};
public:
typedef _MatrixType MatrixType;
/** \deprecated the entire class is deprecated
* Creates a dummy LU factorization object with flags \a flags. */
EIGEN_DEPRECATED SparseLU(int flags = 0)
: m_flags(flags), m_status(0)
{
m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
}
/** \deprecated the entire class is deprecated
* Creates a LU object and compute the respective factorization of \a matrix using
* flags \a flags. */
EIGEN_DEPRECATED SparseLU(const _MatrixType& matrix, int flags = 0)
: /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0)
{
m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
compute(matrix);
}
/** Sets the relative threshold value used to prune zero coefficients during the decomposition.
*
* Setting a value greater than zero speeds up computation, and yields to an imcomplete
* factorization with fewer non zero coefficients. Such approximate factors are especially
* useful to initialize an iterative solver.
*
* Note that the exact meaning of this parameter might depends on the actual
* backend. Moreover, not all backends support this feature.
*
* \sa precision() */
void setPrecision(RealScalar v) { m_precision = v; }
/** \returns the current precision.
*
* \sa setPrecision() */
RealScalar precision() const { return m_precision; }
/** Sets the flags. Possible values are:
* - CompleteFactorization
* - IncompleteFactorization
* - MemoryEfficient
* - one of the ordering methods
* - etc...
*
* \sa flags() */
void setFlags(int f) { m_flags = f; }
/** \returns the current flags */
int flags() const { return m_flags; }
void setOrderingMethod(int m)
{
eigen_assert( (m&~OrderingMask) == 0 && m!=0 && "invalid ordering method");
m_flags = m_flags&~OrderingMask | m&OrderingMask;
}
int orderingMethod() const
{
return m_flags&OrderingMask;
}
/** Computes/re-computes the LU factorization */
void compute(const _MatrixType& matrix);
/** \returns the lower triangular matrix L */
//inline const _MatrixType& matrixL() const { return m_matrixL; }
/** \returns the upper triangular matrix U */
//inline const _MatrixType& matrixU() const { return m_matrixU; }
template<typename BDerived, typename XDerived>
bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
const int transposed = SvNoTrans) const;
/** \returns true if the factorization succeeded */
inline bool succeeded(void) const { return m_succeeded; }
protected:
RealScalar m_precision;
int m_flags;
mutable int m_status;
bool m_succeeded;
};
/** Computes / recomputes the LU decomposition of matrix \a a
* using the default algorithm.
*/
template<typename _MatrixType, typename Backend>
void SparseLU<_MatrixType,Backend>::compute(const _MatrixType& )
{
eigen_assert(false && "not implemented yet");
}
/** Computes *x = U^-1 L^-1 b
*
* If \a transpose is set to SvTranspose or SvAdjoint, the solution
* of the transposed/adjoint system is computed instead.
*
* Not all backends implement the solution of the transposed or
* adjoint system.
*/
template<typename _MatrixType, typename Backend>
template<typename BDerived, typename XDerived>
bool SparseLU<_MatrixType,Backend>::solve(const MatrixBase<BDerived> &, MatrixBase<XDerived>* , const int ) const
{
eigen_assert(false && "not implemented yet");
return false;
}
#endif // EIGEN_SPARSELU_H

100
unsupported/test/cg.cpp
View File

@ -1,100 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.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/>.
#include "sparse.h"
#include <Eigen/IterativeSolvers>
template<typename Scalar,typename Index> void cg(int size)
{
double density = (std::max)(8./(size*size), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
typedef SparseMatrix<Scalar,ColMajor,Index> SparseMatrixType;
SparseMatrixType m2(size,size);
DenseMatrix refMat2(size,size);
DenseVector b = DenseVector::Random(size);
DenseVector ref_x(size), x(size);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, 0, 0);
// for(int i=0; i<rows; ++i)
// m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
SparseMatrixType m3 = m2 * m2.adjoint(), m3_lo(size,size), m3_up(size,size);
DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
m3_lo.template selfadjointView<Lower>().rankUpdate(m2,0);
m3_up.template selfadjointView<Upper>().rankUpdate(m2,0);
ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
x = ConjugateGradient<SparseMatrixType, Lower>().compute(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, full storage, lower");
x.setRandom();
x = ConjugateGradient<SparseMatrixType, Lower>().compute(m3).solveWithGuess(b,x);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solveWithGuess, full storage, lower");
x = ConjugateGradient<SparseMatrixType, Upper>().compute(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, full storage, upper, single dense rhs");
x = ConjugateGradient<SparseMatrixType, Lower>(m3_lo).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, lower only, single dense rhs");
x = ConjugateGradient<SparseMatrixType, Upper>(m3_up).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, upper only, single dense rhs");
x = ConjugateGradient<SparseMatrixType, Lower, IdentityPreconditioner>().compute(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, full storage, lower");
x = ConjugateGradient<SparseMatrixType, Upper, IdentityPreconditioner>().compute(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, full storage, upper, single dense rhs");
x = ConjugateGradient<SparseMatrixType, Lower, IdentityPreconditioner>(m3_lo).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, lower only, single dense rhs");
x = ConjugateGradient<SparseMatrixType, Upper, IdentityPreconditioner>(m3_up).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "ConjugateGradient: solve, upper only, single dense rhs");
ref_x = refMat2.lu().solve(b);
x = BiCGSTAB<SparseMatrixType, IdentityPreconditioner>(m2).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "BiCGSTAB: solve, I, single dense rhs");
x = BiCGSTAB<SparseMatrixType>(m2).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "BiCGSTAB: solve, diag, single dense rhs");
}
void test_cg()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( (cg<double,int>(8)) );
CALL_SUBTEST_1( (cg<double,long int>(8)) );
CALL_SUBTEST_2( (cg<std::complex<double>,int>(internal::random<int>(1,300))) );
CALL_SUBTEST_1( (cg<double,int>(internal::random<int>(1,300))) );
}
}

179
unsupported/test/sparse_ldlt.cpp
View File

@ -1,179 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 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/>.
#define EIGEN_NO_DEPRECATED_WARNING
#include "sparse.h"
#include <Eigen/SparseExtra>
#ifdef EIGEN_CHOLMOD_SUPPORT
#include <Eigen/CholmodSupport>
#endif
template<typename Scalar,typename Index> void sparse_ldlt(int rows, int cols)
{
static bool odd = true;
odd = !odd;
double density = (std::max)(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
typedef SparseMatrix<Scalar,ColMajor,Index> SparseMatrixType;
SparseMatrixType m2(rows, cols);
DenseMatrix refMat2(rows, cols);
DenseVector b = DenseVector::Random(cols);
DenseVector refX(cols), x(cols);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
SparseMatrixType m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
refX = refMat3.template selfadjointView<Upper>().ldlt().solve(b);
typedef SparseMatrix<Scalar,Upper|SelfAdjoint,Index> SparseSelfAdjointMatrix;
x = b;
SparseLDLT<SparseSelfAdjointMatrix> ldlt(m3);
if (ldlt.succeeded())
ldlt.solveInPlace(x);
else
std::cerr << "warning LDLT failed\n";
VERIFY_IS_APPROX(refMat3.template selfadjointView<Upper>() * x, b);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
#ifdef EIGEN_CHOLMOD_SUPPORT
{
x = b;
SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt2(m3);
if (ldlt2.succeeded())
{
ldlt2.solveInPlace(x);
VERIFY_IS_APPROX(refMat3.template selfadjointView<Upper>() * x, b);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solveInPlace");
x = ldlt2.solve(b);
VERIFY_IS_APPROX(refMat3.template selfadjointView<Upper>() * x, b);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solve");
}
else
std::cerr << "warning LDLT failed\n";
}
#endif
// new Simplicial LLT
// new API
{
SparseMatrixType m2(rows, cols);
DenseMatrix refMat2(rows, cols);
DenseVector b = DenseVector::Random(cols);
DenseVector ref_x(cols), x(cols);
DenseMatrix B = DenseMatrix::Random(rows,cols);
DenseMatrix ref_X(rows,cols), X(rows,cols);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0);
for(int i=0; i<rows; ++i)
m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
SparseMatrixType m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
m3_lo.template selfadjointView<Lower>().rankUpdate(m2,0);
m3_up.template selfadjointView<Upper>().rankUpdate(m2,0);
// with a single vector as the rhs
ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
x = SimplicialCholesky<SparseMatrixType, Lower>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, lower, single dense rhs");
x = SimplicialCholesky<SparseMatrixType, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, single dense rhs");
x = SimplicialCholesky<SparseMatrixType, Lower>(m3_lo).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, lower only, single dense rhs");
x = SimplicialCholesky<SparseMatrixType, Upper>(m3_up).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, upper only, single dense rhs");
// with multiple rhs
ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B);
X = SimplicialCholesky<SparseMatrixType, Lower>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, lower, multiple dense rhs");
X = SimplicialCholesky<SparseMatrixType, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, multiple dense rhs");
// with a sparse rhs
SparseMatrixType spB(rows,cols), spX(rows,cols);
B.diagonal().array() += 1;
spB = B.sparseView(0.5,1);
ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB));
spX = SimplicialCholesky<SparseMatrixType, Lower>(m3).solve(spB);
VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: SimplicialCholesky solve, multiple sparse rhs");
//
spX = SimplicialCholesky<SparseMatrixType, Upper>(m3).solve(spB);
VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: SimplicialCholesky solve, multiple sparse rhs");
}
// for(int i=0; i<rows; ++i)
// m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
//
// refX = refMat2.template selfadjointView<Upper>().ldlt().solve(b);
// typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix;
// x = b;
// SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
// if (ldlt.succeeded())
// ldlt.solveInPlace(x);
// else
// std::cerr << "warning LDLT failed\n";
//
// VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
// VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
}
void test_sparse_ldlt()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( (sparse_ldlt<double,int>(8, 8)) );
CALL_SUBTEST_1( (sparse_ldlt<double,long int>(8, 8)) );
int s = internal::random<int>(1,300);
CALL_SUBTEST_2( (sparse_ldlt<std::complex<double>,int>(s,s)) );
CALL_SUBTEST_1( (sparse_ldlt<double,int>(s,s)) );
}
}

144
unsupported/test/sparse_llt.cpp
View File

@ -1,144 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 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/>.
#define EIGEN_NO_DEPRECATED_WARNING
#include "sparse.h"
#include <Eigen/SparseExtra>
#ifdef EIGEN_CHOLMOD_SUPPORT
#include <Eigen/CholmodSupport>
#endif
template<typename Scalar,typename Index> void sparse_llt(int rows, int cols)
{
double density = (std::max)(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
typedef SparseMatrix<Scalar,ColMajor,Index> SparseMatrixType;
// TODO fix the issue with complex (see SparseLLT::solveInPlace)
SparseMatrixType m2(rows, cols);
DenseMatrix refMat2(rows, cols);
DenseVector b = DenseVector::Random(cols);
DenseVector ref_x(cols), x(cols);
DenseMatrix B = DenseMatrix::Random(rows,cols);
DenseMatrix ref_X(rows,cols), X(rows,cols);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0);
for(int i=0; i<rows; ++i)
m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
ref_x = refMat2.template selfadjointView<Lower>().llt().solve(b);
if (!NumTraits<Scalar>::IsComplex)
{
x = b;
SparseLLT<SparseMatrixType > (m2).solveInPlace(x);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: default");
}
#ifdef EIGEN_CHOLMOD_SUPPORT
// legacy API
{
// Cholmod, as configured in CholmodSupport.h, only supports self-adjoint matrices
SparseMatrixType m3 = m2.adjoint()*m2;
DenseMatrix refMat3 = refMat2.adjoint()*refMat2;
ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
x = b;
SparseLLT<SparseMatrixType, Cholmod>(m3).solveInPlace(x);
VERIFY((m3*x).isApprox(b,test_precision<Scalar>()) && "LLT legacy: cholmod solveInPlace");
x = SparseLLT<SparseMatrixType, Cholmod>(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT legacy: cholmod solve");
}
// new API
{
// Cholmod, as configured in CholmodSupport.h, only supports self-adjoint matrices
SparseMatrixType m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
m3_lo.template selfadjointView<Lower>().rankUpdate(m2,0);
m3_up.template selfadjointView<Upper>().rankUpdate(m2,0);
// with a single vector as the rhs
ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
x = CholmodDecomposition<SparseMatrixType, Lower>(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
x = CholmodDecomposition<SparseMatrixType, Upper>(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
x = CholmodDecomposition<SparseMatrixType, Lower>(m3_lo).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
x = CholmodDecomposition<SparseMatrixType, Upper>(m3_up).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
// with multiple rhs
ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B);
#ifndef EIGEN_DEFAULT_TO_ROW_MAJOR
// TODO make sure the API is properly documented about this fact
X = CholmodDecomposition<SparseMatrixType, Lower>(m3).solve(B);
VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "LLT: cholmod solve, multiple dense rhs");
X = CholmodDecomposition<SparseMatrixType, Upper>(m3).solve(B);
VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "LLT: cholmod solve, multiple dense rhs");
#endif
// with a sparse rhs
SparseMatrixType spB(rows,cols), spX(rows,cols);
B.diagonal().array() += 1;
spB = B.sparseView(0.5,1);
ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB));
spX = CholmodDecomposition<SparseMatrixType, Lower>(m3).solve(spB);
VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
spX = CholmodDecomposition<SparseMatrixType, Upper>(m3).solve(spB);
VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
}
#endif
}
void test_sparse_llt()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( (sparse_llt<double,int>(8, 8)) );
int s = internal::random<int>(1,300);
CALL_SUBTEST_2( (sparse_llt<std::complex<double>,int>(s,s)) );
CALL_SUBTEST_1( (sparse_llt<double,int>(s,s)) );
CALL_SUBTEST_1( (sparse_llt<double,long int>(s,s)) );
}
}

128
unsupported/test/sparse_lu_legacy.cpp
View File

@ -1,128 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 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/>.
#define EIGEN_NO_DEPRECATED_WARNING
#include "sparse.h"
#include <Eigen/SparseExtra>
#ifdef EIGEN_UMFPACK_SUPPORT
#include <Eigen/UmfPackSupport>
#endif
#ifdef EIGEN_SUPERLU_SUPPORT
#include <Eigen/SuperLUSupport>
#endif
template<typename Scalar> void sparse_lu_legacy(int rows, int cols)
{
double density = (std::max)(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
DenseVector vec1 = DenseVector::Random(rows);
std::vector<Vector2i> zeroCoords;
std::vector<Vector2i> nonzeroCoords;
SparseMatrix<Scalar> m2(rows, cols);
DenseMatrix refMat2(rows, cols);
DenseVector b = DenseVector::Random(cols);
DenseVector refX(cols), x(cols);
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
FullPivLU<DenseMatrix> refLu(refMat2);
refX = refLu.solve(b);
#if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT)
Scalar refDet = refLu.determinant();
#endif
x.setZero();
// // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x);
// // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default");
#ifdef EIGEN_UMFPACK_SUPPORT
{
// check solve
x.setZero();
SparseLU<SparseMatrix<Scalar>,UmfPack> lu(m2);
VERIFY(lu.succeeded() && "umfpack LU decomposition failed");
VERIFY(lu.solve(b,&x) && "umfpack LU solving failed");
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: umfpack");
VERIFY_IS_APPROX(refDet,lu.determinant());
// TODO check the extracted data
//std::cerr << slu.matrixL() << "\n";
}
#endif
#ifdef EIGEN_SUPERLU_SUPPORT
// legacy, deprecated API
{
x.setZero();
SparseLU<SparseMatrix<Scalar>,SuperLULegacy> slu(m2);
if (slu.succeeded())
{
DenseVector oldb = b;
if (slu.solve(b,&x)) {
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU");
}
else
std::cerr << "super lu solving failed\n";
VERIFY(oldb.isApprox(b) && "the rhs should not be modified!");
// std::cerr << refDet << " == " << slu.determinant() << "\n";
if (slu.solve(b, &x, SvTranspose)) {
VERIFY(b.isApprox(m2.transpose() * x, test_precision<Scalar>()));
}
else
std::cerr << "super lu solving failed\n";
if (slu.solve(b, &x, SvAdjoint)) {
VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>()));
}
else
std::cerr << "super lu solving failed\n";
if (!NumTraits<Scalar>::IsComplex) {
VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
}
}
else
std::cerr << "super lu factorize failed\n";
}
#endif
}
void test_sparse_lu_legacy()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(sparse_lu_legacy<double>(8, 8) );
int s = internal::random<int>(1,300);
CALL_SUBTEST_1(sparse_lu_legacy<std::complex<double> >(s,s) );
CALL_SUBTEST_1(sparse_lu_legacy<double>(s,s) );
}
}