devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Apply clang-format

Browse Source
This commit is contained in:
Tobias Wood
2023-11-29 11:12:48 +00:00
parent 9ea520fc45
commit f38e16c193
534 changed files with 103368 additions and 116934 deletions
Download Patch File Download Diff File Expand all files Collapse all files

627
unsupported/Eigen/src/IterativeSolvers/DGMRES.h
View File

@@ -15,68 +15,62 @@
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
namespace Eigen {
template< typename MatrixType_,
typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
namespace Eigen {
template <typename MatrixType_, typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
class DGMRES;
namespace internal {
template< typename MatrixType_, typename Preconditioner_>
struct traits<DGMRES<MatrixType_,Preconditioner_> >
{
template <typename MatrixType_, typename Preconditioner_>
struct traits<DGMRES<MatrixType_, Preconditioner_> > {
typedef MatrixType_ MatrixType;
typedef Preconditioner_ Preconditioner;
};
/** \brief Computes a permutation vector to have a sorted sequence
* \param vec The vector to reorder.
* \param perm gives the sorted sequence on output. Must be initialized with 0..n-1
* \param ncut Put the ncut smallest elements at the end of the vector
* WARNING This is an expensive sort, so should be used only
* for small size vectors
* TODO Use modified QuickSplit or std::nth_element to get the smallest values
*/
* \param vec The vector to reorder.
* \param perm gives the sorted sequence on output. Must be initialized with 0..n-1
* \param ncut Put the ncut smallest elements at the end of the vector
* WARNING This is an expensive sort, so should be used only
* for small size vectors
* TODO Use modified QuickSplit or std::nth_element to get the smallest values
*/
template <typename VectorType, typename IndexType>
void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut)
{
void sortWithPermutation(VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut) {
eigen_assert(vec.size() == perm.size());
bool flag;
for (Index k = 0; k < ncut; k++)
{
bool flag;
for (Index k = 0; k < ncut; k++) {
flag = false;
for (Index j = 0; j < vec.size()-1; j++)
{
if ( vec(perm(j)) < vec(perm(j+1)) )
{
std::swap(perm(j),perm(j+1));
for (Index j = 0; j < vec.size() - 1; j++) {
if (vec(perm(j)) < vec(perm(j + 1))) {
std::swap(perm(j), perm(j + 1));
flag = true;
}
if (!flag) break; // The vector is in sorted order
if (!flag) break; // The vector is in sorted order
}
}
}
}
} // namespace internal
/**
* \ingroup IterativeLinearSolvers_Module
* \brief A Restarted GMRES with deflation.
* This class implements a modification of the GMRES solver for
* sparse linear systems. The basis is built with modified
* sparse linear systems. The basis is built with modified
* Gram-Schmidt. At each restart, a few approximated eigenvectors
* corresponding to the smallest eigenvalues are used to build a
* preconditioner for the next cycle. This preconditioner
* for deflation can be combined with any other preconditioner,
* the IncompleteLUT for instance. The preconditioner is applied
* preconditioner for the next cycle. This preconditioner
* for deflation can be combined with any other preconditioner,
* the IncompleteLUT for instance. The preconditioner is applied
* at right of the matrix and the combination is multiplicative.
*
*
* \tparam MatrixType_ the type of the sparse matrix A, can be a dense or a sparse matrix.
* \tparam Preconditioner_ the type of the preconditioner. Default is DiagonalPreconditioner
* Typical usage :
* \code
* SparseMatrix<double> A;
* VectorXd x, b;
* VectorXd x, b;
* //Fill A and b ...
* DGMRES<SparseMatrix<double> > solver;
* solver.set_restart(30); // Set restarting value
@@ -84,195 +78,198 @@ void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::
* solver.compute(A);
* x = solver.solve(b);
* \endcode
*
*
* DGMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
*
* References :
* [1] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid
* Algebraic Solvers for Linear Systems Arising from Compressible
* Flows, Computers and Fluids, In Press,
* https://doi.org/10.1016/j.compfluid.2012.03.023
* [2] K. Burrage and J. Erhel, On the performance of various
* https://doi.org/10.1016/j.compfluid.2012.03.023
* [2] K. Burrage and J. Erhel, On the performance of various
* adaptive preconditioned GMRES strategies, 5(1998), 101-121.
* [3] J. Erhel, K. Burrage and B. Pohl, Restarted GMRES
* [3] J. Erhel, K. Burrage and B. Pohl, Restarted GMRES
* preconditioned by deflation,J. Computational and Applied
* Mathematics, 69(1996), 303-318.
* Mathematics, 69(1996), 303-318.
*
*
*/
template< typename MatrixType_, typename Preconditioner_>
class DGMRES : public IterativeSolverBase<DGMRES<MatrixType_,Preconditioner_> >
{
typedef IterativeSolverBase<DGMRES> Base;
using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
using Base::m_isInitialized;
using Base::m_tolerance;
public:
using Base::_solve_impl;
using Base::_solve_with_guess_impl;
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::StorageIndex StorageIndex;
typedef typename MatrixType::RealScalar RealScalar;
typedef Preconditioner_ Preconditioner;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
typedef Matrix<RealScalar,Dynamic,1> DenseRealVector;
typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
template <typename MatrixType_, typename Preconditioner_>
class DGMRES : public IterativeSolverBase<DGMRES<MatrixType_, Preconditioner_> > {
typedef IterativeSolverBase<DGMRES> Base;
using Base::m_error;
using Base::m_info;
using Base::m_isInitialized;
using Base::m_iterations;
using Base::m_tolerance;
using Base::matrix;
public:
using Base::_solve_impl;
using Base::_solve_with_guess_impl;
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::StorageIndex StorageIndex;
typedef typename MatrixType::RealScalar RealScalar;
typedef Preconditioner_ Preconditioner;
typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
typedef Matrix<RealScalar, Dynamic, Dynamic> DenseRealMatrix;
typedef Matrix<Scalar, Dynamic, 1> DenseVector;
typedef Matrix<RealScalar, Dynamic, 1> DenseRealVector;
typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
/** Default constructor. */
DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
DGMRES()
: Base(), m_restart(30), m_neig(0), m_r(0), m_maxNeig(5), m_isDeflAllocated(false), m_isDeflInitialized(false) {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving.
*
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
template<typename MatrixDerived>
explicit DGMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
*
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
template <typename MatrixDerived>
explicit DGMRES(const EigenBase<MatrixDerived>& A)
: Base(A.derived()),
m_restart(30),
m_neig(0),
m_r(0),
m_maxNeig(5),
m_isDeflAllocated(false),
m_isDeflInitialized(false) {}
~DGMRES() {}
/** \internal */
template<typename Rhs,typename Dest>
void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
{
EIGEN_STATIC_ASSERT(Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1, YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
template <typename Rhs, typename Dest>
void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const {
EIGEN_STATIC_ASSERT(Rhs::ColsAtCompileTime == 1 || Dest::ColsAtCompileTime == 1,
YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
dgmres(matrix(), b, x, Base::m_preconditioner);
}
/**
/**
* Get the restart value
*/
*/
Index restart() { return m_restart; }
/**
* Set the restart value (default is 30)
/**
* Set the restart value (default is 30)
*/
void set_restart(const Index restart) { m_restart=restart; }
/**
* Set the number of eigenvalues to deflate at each restart
void set_restart(const Index restart) { m_restart = restart; }
/**
* Set the number of eigenvalues to deflate at each restart
*/
void setEigenv(const Index neig)
{
void setEigenv(const Index neig) {
m_neig = neig;
if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
if (neig + 1 > m_maxNeig) m_maxNeig = neig + 1; // To allow for complex conjugates
}
/**
/**
* Get the size of the deflation subspace size
*/
Index deflSize() {return m_r; }
*/
Index deflSize() { return m_r; }
/**
* Set the maximum size of the deflation subspace
*/
void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; }
protected:
// DGMRES algorithm
template<typename Rhs, typename Dest>
void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
// Perform one cycle of GMRES
template<typename Dest>
Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const;
// Compute data to use for deflation
Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const;
// Apply deflation to a vector
template<typename RhsType, typename DestType>
Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
// Init data for deflation
void dgmresInitDeflation(Index& rows) const;
mutable DenseMatrix m_V; // Krylov basis vectors
mutable DenseMatrix m_H; // Hessenberg matrix
mutable DenseMatrix m_Hes; // Initial hessenberg matrix without Givens rotations applied
mutable Index m_restart; // Maximum size of the Krylov subspace
mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace
mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart
mutable Index m_r; // Current number of deflated eigenvalues, size of m_U
mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate
mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
mutable bool m_isDeflAllocated;
mutable bool m_isDeflInitialized;
//Adaptive strategy
mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
mutable bool m_force; // Force the use of deflation at each restart
};
/**
* \brief Perform several cycles of restarted GMRES with modified Gram Schmidt,
*
protected:
// DGMRES algorithm
template <typename Rhs, typename Dest>
void dgmres(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
// Perform one cycle of GMRES
template <typename Dest>
Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta,
const RealScalar& normRhs, Index& nbIts) const;
// Compute data to use for deflation
Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it,
StorageIndex& neig) const;
// Apply deflation to a vector
template <typename RhsType, typename DestType>
Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
// Init data for deflation
void dgmresInitDeflation(Index& rows) const;
mutable DenseMatrix m_V; // Krylov basis vectors
mutable DenseMatrix m_H; // Hessenberg matrix
mutable DenseMatrix m_Hes; // Initial hessenberg matrix without Givens rotations applied
mutable Index m_restart; // Maximum size of the Krylov subspace
mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace
mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
mutable StorageIndex m_neig; // Number of eigenvalues to extract at each restart
mutable Index m_r; // Current number of deflated eigenvalues, size of m_U
mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate
mutable RealScalar m_lambdaN; // Modulus of the largest eigenvalue of A
mutable bool m_isDeflAllocated;
mutable bool m_isDeflInitialized;
// Adaptive strategy
mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
mutable bool m_force; // Force the use of deflation at each restart
};
/**
* \brief Perform several cycles of restarted GMRES with modified Gram Schmidt,
*
* A right preconditioner is used combined with deflation.
*
*
*/
template< typename MatrixType_, typename Preconditioner_>
template<typename Rhs, typename Dest>
void DGMRES<MatrixType_, Preconditioner_>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
const Preconditioner& precond) const
{
template <typename MatrixType_, typename Preconditioner_>
template <typename Rhs, typename Dest>
void DGMRES<MatrixType_, Preconditioner_>::dgmres(const MatrixType& mat, const Rhs& rhs, Dest& x,
const Preconditioner& precond) const {
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
RealScalar normRhs = rhs.norm();
if(normRhs <= considerAsZero)
{
if (normRhs <= considerAsZero) {
x.setZero();
m_error = 0;
return;
}
//Initialization
// Initialization
m_isDeflInitialized = false;
Index n = mat.rows();
DenseVector r0(n);
Index nbIts = 0;
m_H.resize(m_restart+1, m_restart);
Index n = mat.rows();
DenseVector r0(n);
Index nbIts = 0;
m_H.resize(m_restart + 1, m_restart);
m_Hes.resize(m_restart, m_restart);
m_V.resize(n,m_restart+1);
//Initial residual vector and initial norm
if(x.squaredNorm()==0)
x = precond.solve(rhs);
r0 = rhs - mat * x;
RealScalar beta = r0.norm();
m_error = beta/normRhs;
if(m_error < m_tolerance)
m_info = Success;
m_V.resize(n, m_restart + 1);
// Initial residual vector and initial norm
if (x.squaredNorm() == 0) x = precond.solve(rhs);
r0 = rhs - mat * x;
RealScalar beta = r0.norm();
m_error = beta / normRhs;
if (m_error < m_tolerance)
m_info = Success;
else
m_info = NoConvergence;
// Iterative process
while (nbIts < m_iterations && m_info == NoConvergence)
{
dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
// Compute the new residual vector for the restart
while (nbIts < m_iterations && m_info == NoConvergence) {
dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
// Compute the new residual vector for the restart
if (nbIts < m_iterations && m_info == NoConvergence) {
r0 = rhs - mat * x;
beta = r0.norm();
}
}
}
}
/**
* \brief Perform one restart cycle of DGMRES
@@ -284,231 +281,211 @@ void DGMRES<MatrixType_, Preconditioner_>::dgmres(const MatrixType& mat,const Rh
* \param normRhs The norm of the right hand side vector
* \param nbIts The number of iterations
*/
template< typename MatrixType_, typename Preconditioner_>
template<typename Dest>
Index DGMRES<MatrixType_, Preconditioner_>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const
{
//Initialization
DenseVector g(m_restart+1); // Right hand side of the least square problem
g.setZero();
g(0) = Scalar(beta);
m_V.col(0) = r0/beta;
m_info = NoConvergence;
std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
Index it = 0; // Number of inner iterations
template <typename MatrixType_, typename Preconditioner_>
template <typename Dest>
Index DGMRES<MatrixType_, Preconditioner_>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x,
DenseVector& r0, RealScalar& beta, const RealScalar& normRhs,
Index& nbIts) const {
// Initialization
DenseVector g(m_restart + 1); // Right hand side of the least square problem
g.setZero();
g(0) = Scalar(beta);
m_V.col(0) = r0 / beta;
m_info = NoConvergence;
std::vector<JacobiRotation<Scalar> > gr(m_restart); // Givens rotations
Index it = 0; // Number of inner iterations
Index n = mat.rows();
DenseVector tv1(n), tv2(n); //Temporary vectors
while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
{
DenseVector tv1(n), tv2(n); // Temporary vectors
while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations) {
// Apply preconditioner(s) at right
if (m_isDeflInitialized )
{
dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
tv2 = precond.solve(tv1);
if (m_isDeflInitialized) {
dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
tv2 = precond.solve(tv1);
} else {
tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
}
else
{
tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
}
tv1 = mat * tv2;
tv1 = mat * tv2;
// Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
Scalar coef;
for (Index i = 0; i <= it; ++i)
{
Scalar coef;
for (Index i = 0; i <= it; ++i) {
coef = tv1.dot(m_V.col(i));
tv1 = tv1 - coef * m_V.col(i);
m_H(i,it) = coef;
m_Hes(i,it) = coef;
tv1 = tv1 - coef * m_V.col(i);
m_H(i, it) = coef;
m_Hes(i, it) = coef;
}
// Normalize the vector
coef = tv1.norm();
m_V.col(it+1) = tv1/coef;
m_H(it+1, it) = coef;
// m_Hes(it+1,it) = coef;
// FIXME Check for happy breakdown
// Normalize the vector
coef = tv1.norm();
m_V.col(it + 1) = tv1 / coef;
m_H(it + 1, it) = coef;
// m_Hes(it+1,it) = coef;
// FIXME Check for happy breakdown
// Update Hessenberg matrix with Givens rotations
for (Index i = 1; i <= it; ++i)
{
m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
for (Index i = 1; i <= it; ++i) {
m_H.col(it).applyOnTheLeft(i - 1, i, gr[i - 1].adjoint());
}
// Compute the new plane rotation
gr[it].makeGivens(m_H(it, it), m_H(it+1,it));
// Compute the new plane rotation
gr[it].makeGivens(m_H(it, it), m_H(it + 1, it));
// Apply the new rotation
m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
g.applyOnTheLeft(it,it+1, gr[it].adjoint());
beta = std::abs(g(it+1));
m_error = beta/normRhs;
m_H.col(it).applyOnTheLeft(it, it + 1, gr[it].adjoint());
g.applyOnTheLeft(it, it + 1, gr[it].adjoint());
beta = std::abs(g(it + 1));
m_error = beta / normRhs;
// std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
it++; nbIts++;
if (m_error < m_tolerance)
{
it++;
nbIts++;
if (m_error < m_tolerance) {
// The method has converged
m_info = Success;
break;
}
}
// Compute the new coefficients by solving the least square problem
// it++;
//FIXME Check first if the matrix is singular ... zero diagonal
DenseVector nrs(m_restart);
nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it));
// it++;
// FIXME Check first if the matrix is singular ... zero diagonal
DenseVector nrs(m_restart);
nrs = m_H.topLeftCorner(it, it).template triangularView<Upper>().solve(g.head(it));
// Form the new solution
if (m_isDeflInitialized)
{
tv1 = m_V.leftCols(it) * nrs;
dgmresApplyDeflation(tv1, tv2);
if (m_isDeflInitialized) {
tv1 = m_V.leftCols(it) * nrs;
dgmresApplyDeflation(tv1, tv2);
x = x + precond.solve(tv2);
}
else
x = x + precond.solve(m_V.leftCols(it) * nrs);
} else
x = x + precond.solve(m_V.leftCols(it) * nrs);
// Go for a new cycle and compute data for deflation
if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
dgmresComputeDeflationData(mat, precond, it, m_neig);
return 0;
if (nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r + m_neig) < m_maxNeig)
dgmresComputeDeflationData(mat, precond, it, m_neig);
return 0;
}
template< typename MatrixType_, typename Preconditioner_>
void DGMRES<MatrixType_, Preconditioner_>::dgmresInitDeflation(Index& rows) const
{
template <typename MatrixType_, typename Preconditioner_>
void DGMRES<MatrixType_, Preconditioner_>::dgmresInitDeflation(Index& rows) const {
m_U.resize(rows, m_maxNeig);
m_MU.resize(rows, m_maxNeig);
m_MU.resize(rows, m_maxNeig);
m_T.resize(m_maxNeig, m_maxNeig);
m_lambdaN = 0.0;
m_isDeflAllocated = true;
m_lambdaN = 0.0;
m_isDeflAllocated = true;
}
template< typename MatrixType_, typename Preconditioner_>
inline typename DGMRES<MatrixType_, Preconditioner_>::ComplexVector DGMRES<MatrixType_, Preconditioner_>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const
{
template <typename MatrixType_, typename Preconditioner_>
inline typename DGMRES<MatrixType_, Preconditioner_>::ComplexVector DGMRES<MatrixType_, Preconditioner_>::schurValues(
const ComplexSchur<DenseMatrix>& schurofH) const {
return schurofH.matrixT().diagonal();
}
template< typename MatrixType_, typename Preconditioner_>
inline typename DGMRES<MatrixType_, Preconditioner_>::ComplexVector DGMRES<MatrixType_, Preconditioner_>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
{
template <typename MatrixType_, typename Preconditioner_>
inline typename DGMRES<MatrixType_, Preconditioner_>::ComplexVector DGMRES<MatrixType_, Preconditioner_>::schurValues(
const RealSchur<DenseMatrix>& schurofH) const {
const DenseMatrix& T = schurofH.matrixT();
Index it = T.rows();
ComplexVector eig(it);
Index j = 0;
while (j < it-1)
{
if (T(j+1,j) ==Scalar(0))
{
eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
j++;
}
else
{
eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j));
eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1));
while (j < it - 1) {
if (T(j + 1, j) == Scalar(0)) {
eig(j) = std::complex<RealScalar>(T(j, j), RealScalar(0));
j++;
} else {
eig(j) = std::complex<RealScalar>(T(j, j), T(j + 1, j));
eig(j + 1) = std::complex<RealScalar>(T(j, j + 1), T(j + 1, j + 1));
j++;
}
}
if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
if (j < it - 1) eig(j) = std::complex<RealScalar>(T(j, j), RealScalar(0));
return eig;
}
template< typename MatrixType_, typename Preconditioner_>
Index DGMRES<MatrixType_, Preconditioner_>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const
{
template <typename MatrixType_, typename Preconditioner_>
Index DGMRES<MatrixType_, Preconditioner_>::dgmresComputeDeflationData(const MatrixType& mat,
const Preconditioner& precond, const Index& it,
StorageIndex& neig) const {
// First, find the Schur form of the Hessenberg matrix H
std::conditional_t<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> > schurofH;
bool computeU = true;
DenseMatrix matrixQ(it,it);
DenseMatrix matrixQ(it, it);
matrixQ.setIdentity();
schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU);
schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it, it), matrixQ, computeU);
ComplexVector eig(it);
Matrix<StorageIndex,Dynamic,1>perm(it);
Matrix<StorageIndex, Dynamic, 1> perm(it);
eig = this->schurValues(schurofH);
// Reorder the absolute values of Schur values
DenseRealVector modulEig(it);
for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j));
perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(it-1));
DenseRealVector modulEig(it);
for (Index j = 0; j < it; ++j) modulEig(j) = std::abs(eig(j));
perm.setLinSpaced(it, 0, internal::convert_index<StorageIndex>(it - 1));
internal::sortWithPermutation(modulEig, perm, neig);
if (!m_lambdaN)
{
if (!m_lambdaN) {
m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
}
//Count the real number of extracted eigenvalues (with complex conjugates)
Index nbrEig = 0;
while (nbrEig < neig)
{
if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++;
else nbrEig += 2;
// Count the real number of extracted eigenvalues (with complex conjugates)
Index nbrEig = 0;
while (nbrEig < neig) {
if (eig(perm(it - nbrEig - 1)).imag() == RealScalar(0))
nbrEig++;
else
nbrEig += 2;
}
// Extract the Schur vectors corresponding to the smallest Ritz values
DenseMatrix Sr(it, nbrEig);
DenseMatrix Sr(it, nbrEig);
Sr.setZero();
for (Index j = 0; j < nbrEig; j++)
{
Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
for (Index j = 0; j < nbrEig; j++) {
Sr.col(j) = schurofH.matrixU().col(perm(it - j - 1));
}
// Form the Schur vectors of the initial matrix using the Krylov basis
DenseMatrix X;
DenseMatrix X;
X = m_V.leftCols(it) * Sr;
if (m_r)
{
// Orthogonalize X against m_U using modified Gram-Schmidt
for (Index j = 0; j < nbrEig; j++)
for (Index k =0; k < m_r; k++)
X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k);
if (m_r) {
// Orthogonalize X against m_U using modified Gram-Schmidt
for (Index j = 0; j < nbrEig; j++)
for (Index k = 0; k < m_r; k++) X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j))) * m_U.col(k);
}
// Compute m_MX = A * M^-1 * X
Index m = m_V.rows();
if (!m_isDeflAllocated)
dgmresInitDeflation(m);
if (!m_isDeflAllocated) dgmresInitDeflation(m);
DenseMatrix MX(m, nbrEig);
DenseVector tv1(m);
for (Index j = 0; j < nbrEig; j++)
{
for (Index j = 0; j < nbrEig; j++) {
tv1 = mat * X.col(j);
MX.col(j) = precond.solve(tv1);
}
//Update m_T = [U'MU U'MX; X'MU X'MX]
m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX;
if(m_r)
{
m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX;
// Update m_T = [U'MU U'MX; X'MU X'MX]
m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX;
if (m_r) {
m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX;
m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r);
}
// Save X into m_U and m_MX in m_MU
for (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
for (Index j = 0; j < nbrEig; j++) m_U.col(m_r + j) = X.col(j);
for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r + j) = MX.col(j);
// Increase the size of the invariant subspace
m_r += nbrEig;
m_r += nbrEig;
// Factorize m_T into m_luT
m_luT.compute(m_T.topLeftCorner(m_r, m_r));
//FIXME CHeck if the factorization was correctly done (nonsingular matrix)
// FIXME CHeck if the factorization was correctly done (nonsingular matrix)
m_isDeflInitialized = true;
return 0;
return 0;
}
template<typename MatrixType_, typename Preconditioner_>
template<typename RhsType, typename DestType>
Index DGMRES<MatrixType_, Preconditioner_>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
{
DenseVector x1 = m_U.leftCols(m_r).transpose() * x;
y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
return 0;
template <typename MatrixType_, typename Preconditioner_>
template <typename RhsType, typename DestType>
Index DGMRES<MatrixType_, Preconditioner_>::dgmresApplyDeflation(const RhsType& x, DestType& y) const {
DenseVector x1 = m_U.leftCols(m_r).transpose() * x;
y = x + m_U.leftCols(m_r) * (m_lambdaN * m_luT.solve(x1) - x1);
return 0;
}
} // end namespace Eigen
#endif
} // end namespace Eigen
#endif

285
unsupported/Eigen/src/IterativeSolvers/GMRES.h
View File

@@ -19,58 +19,56 @@ namespace Eigen {
namespace internal {
/**
* Generalized Minimal Residual Algorithm based on the
* Arnoldi algorithm implemented with Householder reflections.
*
* Parameters:
* \param mat matrix of linear system of equations
* \param rhs right hand side vector of linear system of equations
* \param x on input: initial guess, on output: solution
* \param precond preconditioner used
* \param iters on input: maximum number of iterations to perform
* on output: number of iterations performed
* \param restart number of iterations for a restart
* \param tol_error on input: relative residual tolerance
* on output: residuum achieved
*
* \sa IterativeMethods::bicgstab()
*
*
* For references, please see:
*
* Saad, Y. and Schultz, M. H.
* GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems.
* SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869.
*
* Saad, Y.
* Iterative Methods for Sparse Linear Systems.
* Society for Industrial and Applied Mathematics, Philadelphia, 2003.
*
* Walker, H. F.
* Implementations of the GMRES method.
* Comput.Phys.Comm. 53, 1989, pp. 311 - 320.
*
* Walker, H. F.
* Implementation of the GMRES Method using Householder Transformations.
* SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163.
*
*/
template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond,
Index &iters, const Index &restart, typename Dest::RealScalar & tol_error) {
using std::sqrt;
* Generalized Minimal Residual Algorithm based on the
* Arnoldi algorithm implemented with Householder reflections.
*
* Parameters:
* \param mat matrix of linear system of equations
* \param rhs right hand side vector of linear system of equations
* \param x on input: initial guess, on output: solution
* \param precond preconditioner used
* \param iters on input: maximum number of iterations to perform
* on output: number of iterations performed
* \param restart number of iterations for a restart
* \param tol_error on input: relative residual tolerance
* on output: residuum achieved
*
* \sa IterativeMethods::bicgstab()
*
*
* For references, please see:
*
* Saad, Y. and Schultz, M. H.
* GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems.
* SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869.
*
* Saad, Y.
* Iterative Methods for Sparse Linear Systems.
* Society for Industrial and Applied Mathematics, Philadelphia, 2003.
*
* Walker, H. F.
* Implementations of the GMRES method.
* Comput.Phys.Comm. 53, 1989, pp. 311 - 320.
*
* Walker, H. F.
* Implementation of the GMRES Method using Householder Transformations.
* SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163.
*
*/
template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
bool gmres(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond, Index& iters,
const Index& restart, typename Dest::RealScalar& tol_error) {
using std::abs;
using std::sqrt;
typedef typename Dest::RealScalar RealScalar;
typedef typename Dest::Scalar Scalar;
typedef Matrix < Scalar, Dynamic, 1 > VectorType;
typedef Matrix < Scalar, Dynamic, Dynamic, ColMajor> FMatrixType;
typedef Matrix<Scalar, Dynamic, 1> VectorType;
typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> FMatrixType;
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
if(rhs.norm() <= considerAsZero)
{
if (rhs.norm() <= considerAsZero) {
x.setZero();
tol_error = 0;
return true;
@@ -83,26 +81,25 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
const Index m = mat.rows();
// residual and preconditioned residual
VectorType p0 = rhs - mat*x;
VectorType p0 = rhs - mat * x;
VectorType r0 = precond.solve(p0);
const RealScalar r0Norm = r0.norm();
// is initial guess already good enough?
if(r0Norm == 0)
{
if (r0Norm == 0) {
tol_error = 0;
return true;
}
// storage for Hessenberg matrix and Householder data
FMatrixType H = FMatrixType::Zero(m, restart + 1);
VectorType w = VectorType::Zero(restart + 1);
VectorType tau = VectorType::Zero(restart + 1);
FMatrixType H = FMatrixType::Zero(m, restart + 1);
VectorType w = VectorType::Zero(restart + 1);
VectorType tau = VectorType::Zero(restart + 1);
// storage for Jacobi rotations
std::vector < JacobiRotation < Scalar > > G(restart);
std::vector<JacobiRotation<Scalar> > G(restart);
// storage for temporaries
VectorType t(m), v(m), workspace(m), x_new(m);
@@ -111,9 +108,8 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
RealScalar beta;
r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
w(0) = Scalar(beta);
for (Index k = 1; k <= restart; ++k)
{
for (Index k = 1; k <= restart; ++k) {
++iters;
v = VectorType::Unit(m, k - 1);
@@ -134,10 +130,8 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
}
if (v.tail(m - k).norm() != 0.0)
{
if (k <= restart)
{
if (v.tail(m - k).norm() != 0.0) {
if (k <= restart) {
// generate new Householder vector
Ref<VectorType> Hk_tail = H.col(k).tail(m - k - 1);
v.tail(m - k).makeHouseholder(Hk_tail, tau.coeffRef(k), beta);
@@ -147,17 +141,14 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
}
}
if (k > 1)
{
for (Index i = 0; i < k - 1; ++i)
{
if (k > 1) {
for (Index i = 0; i < k - 1; ++i) {
// apply old Givens rotations to v
v.applyOnTheLeft(i, i + 1, G[i].adjoint());
}
}
if (k<m && v(k) != (Scalar) 0)
{
if (k < m && v(k) != (Scalar)0) {
// determine next Givens rotation
G[k - 1].makeGivens(v(k - 1), v(k));
@@ -167,21 +158,19 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
}
// insert coefficients into upper matrix triangle
H.col(k-1).head(k) = v.head(k);
H.col(k - 1).head(k) = v.head(k);
tol_error = abs(w(k)) / r0Norm;
bool stop = (k==m || tol_error < tol || iters == maxIters);
bool stop = (k == m || tol_error < tol || iters == maxIters);
if (stop || k == restart)
{
if (stop || k == restart) {
// solve upper triangular system
Ref<VectorType> y = w.head(k);
H.topLeftCorner(k, k).template triangularView <Upper>().solveInPlace(y);
H.topLeftCorner(k, k).template triangularView<Upper>().solveInPlace(y);
// use Horner-like scheme to calculate solution vector
x_new.setZero();
for (Index i = k - 1; i >= 0; --i)
{
for (Index i = k - 1; i >= 0; --i) {
x_new(i) += y(i);
// apply Householder reflection H_{i} to x_new
x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
@@ -189,16 +178,13 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
x += x_new;
if(stop)
{
if (stop) {
return true;
}
else
{
k=0;
} else {
k = 0;
// reset data for restart
p0.noalias() = rhs - mat*x;
p0.noalias() = rhs - mat * x;
r0 = precond.solve(p0);
// clear Hessenberg matrix and Householder data
@@ -214,125 +200,116 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
}
return false;
}
}
} // namespace internal
template< typename MatrixType_,
typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
template <typename MatrixType_, typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
class GMRES;
namespace internal {
template< typename MatrixType_, typename Preconditioner_>
struct traits<GMRES<MatrixType_,Preconditioner_> >
{
template <typename MatrixType_, typename Preconditioner_>
struct traits<GMRES<MatrixType_, Preconditioner_> > {
typedef MatrixType_ MatrixType;
typedef Preconditioner_ Preconditioner;
};
}
} // namespace internal
/** \ingroup IterativeLinearSolvers_Module
* \brief A GMRES solver for sparse square problems
*
* This class allows to solve for A.x = b sparse linear problems using a generalized minimal
* residual method. The vectors x and b can be either dense or sparse.
*
* \tparam MatrixType_ the type of the sparse matrix A, can be a dense or a sparse matrix.
* \tparam Preconditioner_ the type of the preconditioner. Default is DiagonalPreconditioner
*
* The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
* and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
* and NumTraits<Scalar>::epsilon() for the tolerance.
*
* This class can be used as the direct solver classes. Here is a typical usage example:
* \code
* int n = 10000;
* VectorXd x(n), b(n);
* SparseMatrix<double> A(n,n);
* // fill A and b
* GMRES<SparseMatrix<double> > solver(A);
* x = solver.solve(b);
* std::cout << "#iterations: " << solver.iterations() << std::endl;
* std::cout << "estimated error: " << solver.error() << std::endl;
* // update b, and solve again
* x = solver.solve(b);
* \endcode
*
* By default the iterations start with x=0 as an initial guess of the solution.
* One can control the start using the solveWithGuess() method.
*
* GMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
*
* \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/
template< typename MatrixType_, typename Preconditioner_>
class GMRES : public IterativeSolverBase<GMRES<MatrixType_,Preconditioner_> >
{
* \brief A GMRES solver for sparse square problems
*
* This class allows to solve for A.x = b sparse linear problems using a generalized minimal
* residual method. The vectors x and b can be either dense or sparse.
*
* \tparam MatrixType_ the type of the sparse matrix A, can be a dense or a sparse matrix.
* \tparam Preconditioner_ the type of the preconditioner. Default is DiagonalPreconditioner
*
* The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
* and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
* and NumTraits<Scalar>::epsilon() for the tolerance.
*
* This class can be used as the direct solver classes. Here is a typical usage example:
* \code
* int n = 10000;
* VectorXd x(n), b(n);
* SparseMatrix<double> A(n,n);
* // fill A and b
* GMRES<SparseMatrix<double> > solver(A);
* x = solver.solve(b);
* std::cout << "#iterations: " << solver.iterations() << std::endl;
* std::cout << "estimated error: " << solver.error() << std::endl;
* // update b, and solve again
* x = solver.solve(b);
* \endcode
*
* By default the iterations start with x=0 as an initial guess of the solution.
* One can control the start using the solveWithGuess() method.
*
* GMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
*
* \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/
template <typename MatrixType_, typename Preconditioner_>
class GMRES : public IterativeSolverBase<GMRES<MatrixType_, Preconditioner_> > {
typedef IterativeSolverBase<GMRES> Base;
using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
using Base::m_isInitialized;
using Base::m_iterations;
using Base::matrix;
private:
private:
Index m_restart;
public:
public:
using Base::_solve_impl;
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Preconditioner_ Preconditioner;
public:
public:
/** Default constructor. */
GMRES() : Base(), m_restart(30) {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving.
*
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
template<typename MatrixDerived>
*
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
template <typename MatrixDerived>
explicit GMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30) {}
~GMRES() {}
/** Get the number of iterations after that a restart is performed.
*/
*/
Index get_restart() { return m_restart; }
/** Set the number of iterations after that a restart is performed.
* \param restart number of iterations for a restarti, default is 30.
*/
void set_restart(const Index restart) { m_restart=restart; }
* \param restart number of iterations for a restarti, default is 30.
*/
void set_restart(const Index restart) { m_restart = restart; }
/** \internal */
template<typename Rhs,typename Dest>
void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
{
template <typename Rhs, typename Dest>
void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const {
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
bool ret = internal::gmres(matrix(), b, x, Base::m_preconditioner, m_iterations, m_restart, m_error);
m_info = (!ret) ? NumericalIssue
: m_error <= Base::m_tolerance ? Success
: NoConvergence;
m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence;
}
protected:
protected:
};
} // end namespace Eigen
} // end namespace Eigen
#endif // EIGEN_GMRES_H
#endif // EIGEN_GMRES_H

119
unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h
View File

@@ -13,81 +13,76 @@
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
namespace Eigen {
namespace Eigen {
template <typename Scalar_>
class IncompleteLU : public SparseSolverBase<IncompleteLU<Scalar_> >
{
protected:
typedef SparseSolverBase<IncompleteLU<Scalar_> > Base;
using Base::m_isInitialized;
typedef Scalar_ Scalar;
typedef Matrix<Scalar,Dynamic,1> Vector;
typedef typename Vector::Index Index;
typedef SparseMatrix<Scalar,RowMajor> FactorType;
class IncompleteLU : public SparseSolverBase<IncompleteLU<Scalar_> > {
protected:
typedef SparseSolverBase<IncompleteLU<Scalar_> > Base;
using Base::m_isInitialized;
public:
typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
typedef Scalar_ Scalar;
typedef Matrix<Scalar, Dynamic, 1> Vector;
typedef typename Vector::Index Index;
typedef SparseMatrix<Scalar, RowMajor> FactorType;
IncompleteLU() {}
public:
typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
template<typename MatrixType>
IncompleteLU(const MatrixType& mat)
{
compute(mat);
}
IncompleteLU() {}
Index rows() const { return m_lu.rows(); }
Index cols() const { return m_lu.cols(); }
template <typename MatrixType>
IncompleteLU(const MatrixType& mat) {
compute(mat);
}
template<typename MatrixType>
IncompleteLU& compute(const MatrixType& mat)
{
m_lu = mat;
int size = mat.cols();
Vector diag(size);
for(int i=0; i<size; ++i)
{
typename FactorType::InnerIterator k_it(m_lu,i);
for(; k_it && k_it.index()<i; ++k_it)
{
int k = k_it.index();
k_it.valueRef() /= diag(k);
Index rows() const { return m_lu.rows(); }
Index cols() const { return m_lu.cols(); }
typename FactorType::InnerIterator j_it(k_it);
typename FactorType::InnerIterator kj_it(m_lu, k);
while(kj_it && kj_it.index()<=k) ++kj_it;
for(++j_it; j_it; )
{
if(kj_it.index()==j_it.index())
{
j_it.valueRef() -= k_it.value() * kj_it.value();
++j_it;
++kj_it;
}
else if(kj_it.index()<j_it.index()) ++kj_it;
else ++j_it;
}
template <typename MatrixType>
IncompleteLU& compute(const MatrixType& mat) {
m_lu = mat;
int size = mat.cols();
Vector diag(size);
for (int i = 0; i < size; ++i) {
typename FactorType::InnerIterator k_it(m_lu, i);
for (; k_it && k_it.index() < i; ++k_it) {
int k = k_it.index();
k_it.valueRef() /= diag(k);
typename FactorType::InnerIterator j_it(k_it);
typename FactorType::InnerIterator kj_it(m_lu, k);
while (kj_it && kj_it.index() <= k) ++kj_it;
for (++j_it; j_it;) {
if (kj_it.index() == j_it.index()) {
j_it.valueRef() -= k_it.value() * kj_it.value();
++j_it;
++kj_it;
} else if (kj_it.index() < j_it.index())
++kj_it;
else
++j_it;
}
if(k_it && k_it.index()==i) diag(i) = k_it.value();
else diag(i) = 1;
}
m_isInitialized = true;
return *this;
if (k_it && k_it.index() == i)
diag(i) = k_it.value();
else
diag(i) = 1;
}
m_isInitialized = true;
return *this;
}
template<typename Rhs, typename Dest>
void _solve_impl(const Rhs& b, Dest& x) const
{
x = m_lu.template triangularView<UnitLower>().solve(b);
x = m_lu.template triangularView<Upper>().solve(x);
}
template <typename Rhs, typename Dest>
void _solve_impl(const Rhs& b, Dest& x) const {
x = m_lu.template triangularView<UnitLower>().solve(b);
x = m_lu.template triangularView<Upper>().solve(x);
}
protected:
FactorType m_lu;
protected:
FactorType m_lu;
};
} // end namespace Eigen
} // end namespace Eigen
#endif // EIGEN_INCOMPLETE_LU_H
#endif // EIGEN_INCOMPLETE_LU_H

3
unsupported/Eigen/src/IterativeSolvers/InternalHeaderCheck.h
View File

@@ -1,3 +1,4 @@
#ifndef EIGEN_ITERATIVE_SOLVERS_MODULE_H
#error "Please include unsupported/Eigen/IterativeSolvers instead of including headers inside the src directory directly."
#error \
"Please include unsupported/Eigen/IterativeSolvers instead of including headers inside the src directory directly."
#endif

474
unsupported/Eigen/src/IterativeSolvers/MINRES.h
View File

@@ -9,262 +9,246 @@
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MINRES_H
#define EIGEN_MINRES_H
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
namespace Eigen {
namespace internal {
/** \internal Low-level MINRES algorithm
* \param mat The matrix A
* \param rhs The right hand side vector b
* \param x On input and initial solution, on output the computed solution.
* \param precond A right preconditioner being able to efficiently solve for an
* approximation of Ax=b (regardless of b)
* \param iters On input the max number of iteration, on output the number of performed iterations.
* \param tol_error On input the tolerance error, on output an estimation of the relative error.
*/
template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
EIGEN_DONT_INLINE
void minres(const MatrixType& mat, const Rhs& rhs, Dest& x,
const Preconditioner& precond, Index& iters,
typename Dest::RealScalar& tol_error)
{
using std::sqrt;
typedef typename Dest::RealScalar RealScalar;
typedef typename Dest::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,1> VectorType;
// Check for zero rhs
const RealScalar rhsNorm2(rhs.squaredNorm());
if(rhsNorm2 == 0)
{
x.setZero();
iters = 0;
tol_error = 0;
return;
}
// initialize
const Index maxIters(iters); // initialize maxIters to iters
const Index N(mat.cols()); // the size of the matrix
const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
// Initialize preconditioned Lanczos
VectorType v_old(N); // will be initialized inside loop
VectorType v( VectorType::Zero(N) ); //initialize v
VectorType v_new(rhs-mat*x); //initialize v_new
RealScalar residualNorm2(v_new.squaredNorm());
VectorType w(N); // will be initialized inside loop
VectorType w_new(precond.solve(v_new)); // initialize w_new
// RealScalar beta; // will be initialized inside loop
RealScalar beta_new2(v_new.dot(w_new));
eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
RealScalar beta_new(sqrt(beta_new2));
const RealScalar beta_one(beta_new);
// Initialize other variables
RealScalar c(1.0); // the cosine of the Givens rotation
RealScalar c_old(1.0);
RealScalar s(0.0); // the sine of the Givens rotation
RealScalar s_old(0.0); // the sine of the Givens rotation
VectorType p_oold(N); // will be initialized in loop
VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
VectorType p(p_old); // initialize p=0
RealScalar eta(1.0);
iters = 0; // reset iters
while ( iters < maxIters )
{
// Preconditioned Lanczos
/* Note that there are 4 variants on the Lanczos algorithm. These are
* described in Paige, C. C. (1972). Computational variants of
* the Lanczos method for the eigenproblem. IMA Journal of Applied
* Mathematics, 10(3), 373-381. The current implementation corresponds
* to the case A(2,7) in the paper. It also corresponds to
* algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
* Systems, 2003 p.173. For the preconditioned version see
* A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
*/
const RealScalar beta(beta_new);
v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
v_new /= beta_new; // overwrite v_new for next iteration
w_new /= beta_new; // overwrite w_new for next iteration
v = v_new; // update
w = w_new; // update
v_new.noalias() = mat*w - beta*v_old; // compute v_new
const RealScalar alpha = v_new.dot(w);
v_new -= alpha*v; // overwrite v_new
w_new = precond.solve(v_new); // overwrite w_new
beta_new2 = v_new.dot(w_new); // compute beta_new
eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
beta_new = sqrt(beta_new2); // compute beta_new
// Givens rotation
const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration
const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration
const RealScalar r1_hat=c*alpha-c_old*s*beta;
const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) );
c_old = c; // store for next iteration
s_old = s; // store for next iteration
c=r1_hat/r1; // new cosine
s=beta_new/r1; // new sine
// Update solution
p_oold = p_old;
p_old = p;
p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
x += beta_one*c*eta*p;
/* Update the squared residual. Note that this is the estimated residual.
The real residual |Ax-b|^2 may be slightly larger */
residualNorm2 *= s*s;
if ( residualNorm2 < threshold2)
{
break;
}
eta=-s*eta; // update eta
iters++; // increment iteration number (for output purposes)
}
/* Compute error. Note that this is the estimated error. The real
error |Ax-b|/|b| may be slightly larger */
tol_error = std::sqrt(residualNorm2 / rhsNorm2);
}
}
template< typename MatrixType_, int UpLo_=Lower,
typename Preconditioner_ = IdentityPreconditioner>
class MINRES;
namespace internal {
template< typename MatrixType_, int UpLo_, typename Preconditioner_>
struct traits<MINRES<MatrixType_,UpLo_,Preconditioner_> >
{
typedef MatrixType_ MatrixType;
typedef Preconditioner_ Preconditioner;
};
}
/** \ingroup IterativeLinearSolvers_Module
* \brief A minimal residual solver for sparse symmetric problems
*
* This class allows to solve for A.x = b sparse linear problems using the MINRES algorithm
* of Paige and Saunders (1975). The sparse matrix A must be symmetric (possibly indefinite).
* The vectors x and b can be either dense or sparse.
*
* \tparam MatrixType_ the type of the sparse matrix A, can be a dense or a sparse matrix.
* \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower,
* Upper, or Lower|Upper in which the full matrix entries will be considered. Default is Lower.
* \tparam Preconditioner_ the type of the preconditioner. Default is DiagonalPreconditioner
*
* The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
* and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
* and NumTraits<Scalar>::epsilon() for the tolerance.
*
* This class can be used as the direct solver classes. Here is a typical usage example:
* \code
* int n = 10000;
* VectorXd x(n), b(n);
* SparseMatrix<double> A(n,n);
* // fill A and b
* MINRES<SparseMatrix<double> > mr;
* mr.compute(A);
* x = mr.solve(b);
* std::cout << "#iterations: " << mr.iterations() << std::endl;
* std::cout << "estimated error: " << mr.error() << std::endl;
* // update b, and solve again
* x = mr.solve(b);
* \endcode
*
* By default the iterations start with x=0 as an initial guess of the solution.
* One can control the start using the solveWithGuess() method.
*
* MINRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
*
* \sa class ConjugateGradient, BiCGSTAB, SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
namespace internal {
/** \internal Low-level MINRES algorithm
* \param mat The matrix A
* \param rhs The right hand side vector b
* \param x On input and initial solution, on output the computed solution.
* \param precond A right preconditioner being able to efficiently solve for an
* approximation of Ax=b (regardless of b)
* \param iters On input the max number of iteration, on output the number of performed iterations.
* \param tol_error On input the tolerance error, on output an estimation of the relative error.
*/
template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
EIGEN_DONT_INLINE void minres(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond,
Index& iters, typename Dest::RealScalar& tol_error) {
using std::sqrt;
typedef typename Dest::RealScalar RealScalar;
typedef typename Dest::Scalar Scalar;
typedef Matrix<Scalar, Dynamic, 1> VectorType;
// Check for zero rhs
const RealScalar rhsNorm2(rhs.squaredNorm());
if (rhsNorm2 == 0) {
x.setZero();
iters = 0;
tol_error = 0;
return;
}
// initialize
const Index maxIters(iters); // initialize maxIters to iters
const Index N(mat.cols()); // the size of the matrix
const RealScalar threshold2(tol_error * tol_error * rhsNorm2); // convergence threshold (compared to residualNorm2)
// Initialize preconditioned Lanczos
VectorType v_old(N); // will be initialized inside loop
VectorType v(VectorType::Zero(N)); // initialize v
VectorType v_new(rhs - mat * x); // initialize v_new
RealScalar residualNorm2(v_new.squaredNorm());
VectorType w(N); // will be initialized inside loop
VectorType w_new(precond.solve(v_new)); // initialize w_new
// RealScalar beta; // will be initialized inside loop
RealScalar beta_new2(v_new.dot(w_new));
eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
RealScalar beta_new(sqrt(beta_new2));
const RealScalar beta_one(beta_new);
// Initialize other variables
RealScalar c(1.0); // the cosine of the Givens rotation
RealScalar c_old(1.0);
RealScalar s(0.0); // the sine of the Givens rotation
RealScalar s_old(0.0); // the sine of the Givens rotation
VectorType p_oold(N); // will be initialized in loop
VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
VectorType p(p_old); // initialize p=0
RealScalar eta(1.0);
iters = 0; // reset iters
while (iters < maxIters) {
// Preconditioned Lanczos
/* Note that there are 4 variants on the Lanczos algorithm. These are
* described in Paige, C. C. (1972). Computational variants of
* the Lanczos method for the eigenproblem. IMA Journal of Applied
* Mathematics, 10(3), 373-381. The current implementation corresponds
* to the case A(2,7) in the paper. It also corresponds to
* algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
* Systems, 2003 p.173. For the preconditioned version see
* A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
*/
template< typename MatrixType_, int UpLo_, typename Preconditioner_>
class MINRES : public IterativeSolverBase<MINRES<MatrixType_,UpLo_,Preconditioner_> >
{
typedef IterativeSolverBase<MINRES> Base;
using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
using Base::m_isInitialized;
public:
using Base::_solve_impl;
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Preconditioner_ Preconditioner;
enum {UpLo = UpLo_};
public:
/** Default constructor. */
MINRES() : Base() {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving.
*
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
template<typename MatrixDerived>
explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
/** Destructor. */
~MINRES(){}
const RealScalar beta(beta_new);
v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
v_new /= beta_new; // overwrite v_new for next iteration
w_new /= beta_new; // overwrite w_new for next iteration
v = v_new; // update
w = w_new; // update
v_new.noalias() = mat * w - beta * v_old; // compute v_new
const RealScalar alpha = v_new.dot(w);
v_new -= alpha * v; // overwrite v_new
w_new = precond.solve(v_new); // overwrite w_new
beta_new2 = v_new.dot(w_new); // compute beta_new
eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
beta_new = sqrt(beta_new2); // compute beta_new
/** \internal */
template<typename Rhs,typename Dest>
void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
{
typedef typename Base::MatrixWrapper MatrixWrapper;
typedef typename Base::ActualMatrixType ActualMatrixType;
enum {
TransposeInput = (!MatrixWrapper::MatrixFree)
&& (UpLo==(Lower|Upper))
&& (!MatrixType::IsRowMajor)
&& (!NumTraits<Scalar>::IsComplex)
};
typedef std::conditional_t<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&> RowMajorWrapper;
EIGEN_STATIC_ASSERT(internal::check_implication(MatrixWrapper::MatrixFree, UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
typedef std::conditional_t<UpLo==(Lower|Upper),
RowMajorWrapper,
typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
> SelfAdjointWrapper;
// Givens rotation
const RealScalar r2 = s * alpha + c * c_old * beta; // s, s_old, c and c_old are still from previous iteration
const RealScalar r3 = s_old * beta; // s, s_old, c and c_old are still from previous iteration
const RealScalar r1_hat = c * alpha - c_old * s * beta;
const RealScalar r1 = sqrt(std::pow(r1_hat, 2) + std::pow(beta_new, 2));
c_old = c; // store for next iteration
s_old = s; // store for next iteration
c = r1_hat / r1; // new cosine
s = beta_new / r1; // new sine
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
RowMajorWrapper row_mat(matrix());
internal::minres(SelfAdjointWrapper(row_mat), b, x,
Base::m_preconditioner, m_iterations, m_error);
m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
}
protected:
// Update solution
p_oold = p_old;
p_old = p;
p.noalias() = (w - r2 * p_old - r3 * p_oold) / r1; // IS NOALIAS REQUIRED?
x += beta_one * c * eta * p;
/* Update the squared residual. Note that this is the estimated residual.
The real residual |Ax-b|^2 may be slightly larger */
residualNorm2 *= s * s;
if (residualNorm2 < threshold2) {
break;
}
eta = -s * eta; // update eta
iters++; // increment iteration number (for output purposes)
}
/* Compute error. Note that this is the estimated error. The real
error |Ax-b|/|b| may be slightly larger */
tol_error = std::sqrt(residualNorm2 / rhsNorm2);
}
} // namespace internal
template <typename MatrixType_, int UpLo_ = Lower, typename Preconditioner_ = IdentityPreconditioner>
class MINRES;
namespace internal {
template <typename MatrixType_, int UpLo_, typename Preconditioner_>
struct traits<MINRES<MatrixType_, UpLo_, Preconditioner_> > {
typedef MatrixType_ MatrixType;
typedef Preconditioner_ Preconditioner;
};
} // namespace internal
/** \ingroup IterativeLinearSolvers_Module
* \brief A minimal residual solver for sparse symmetric problems
*
* This class allows to solve for A.x = b sparse linear problems using the MINRES algorithm
* of Paige and Saunders (1975). The sparse matrix A must be symmetric (possibly indefinite).
* The vectors x and b can be either dense or sparse.
*
* \tparam MatrixType_ the type of the sparse matrix A, can be a dense or a sparse matrix.
* \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower,
* Upper, or Lower|Upper in which the full matrix entries will be considered. Default is Lower.
* \tparam Preconditioner_ the type of the preconditioner. Default is DiagonalPreconditioner
*
* The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
* and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
* and NumTraits<Scalar>::epsilon() for the tolerance.
*
* This class can be used as the direct solver classes. Here is a typical usage example:
* \code
* int n = 10000;
* VectorXd x(n), b(n);
* SparseMatrix<double> A(n,n);
* // fill A and b
* MINRES<SparseMatrix<double> > mr;
* mr.compute(A);
* x = mr.solve(b);
* std::cout << "#iterations: " << mr.iterations() << std::endl;
* std::cout << "estimated error: " << mr.error() << std::endl;
* // update b, and solve again
* x = mr.solve(b);
* \endcode
*
* By default the iterations start with x=0 as an initial guess of the solution.
* One can control the start using the solveWithGuess() method.
*
* MINRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
*
* \sa class ConjugateGradient, BiCGSTAB, SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/
template <typename MatrixType_, int UpLo_, typename Preconditioner_>
class MINRES : public IterativeSolverBase<MINRES<MatrixType_, UpLo_, Preconditioner_> > {
typedef IterativeSolverBase<MINRES> Base;
using Base::m_error;
using Base::m_info;
using Base::m_isInitialized;
using Base::m_iterations;
using Base::matrix;
public:
using Base::_solve_impl;
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Preconditioner_ Preconditioner;
enum { UpLo = UpLo_ };
public:
/** Default constructor. */
MINRES() : Base() {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving.
*
* This constructor is a shortcut for the default constructor followed
* by a call to compute().
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
template <typename MatrixDerived>
explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
/** Destructor. */
~MINRES() {}
/** \internal */
template <typename Rhs, typename Dest>
void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const {
typedef typename Base::MatrixWrapper MatrixWrapper;
typedef typename Base::ActualMatrixType ActualMatrixType;
enum {
TransposeInput = (!MatrixWrapper::MatrixFree) && (UpLo == (Lower | Upper)) && (!MatrixType::IsRowMajor) &&
(!NumTraits<Scalar>::IsComplex)
};
typedef std::conditional_t<TransposeInput, Transpose<const ActualMatrixType>, ActualMatrixType const&>
RowMajorWrapper;
EIGEN_STATIC_ASSERT(internal::check_implication(MatrixWrapper::MatrixFree, UpLo == (Lower | Upper)),
MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
typedef std::conditional_t<UpLo == (Lower | Upper), RowMajorWrapper,
typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type>
SelfAdjointWrapper;
} // end namespace Eigen
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
RowMajorWrapper row_mat(matrix());
internal::minres(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, m_error);
m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
}
#endif // EIGEN_MINRES_H
protected:
};
} // end namespace Eigen
#endif // EIGEN_MINRES_H

324
unsupported/Eigen/src/IterativeSolvers/Scaling.h
View File

@@ -16,181 +16,157 @@
namespace Eigen {
/**
* \ingroup IterativeSolvers_Module
* \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
*
* This class can be used as a preprocessing tool to accelerate the convergence of iterative methods
*
* This feature is useful to limit the pivoting amount during LU/ILU factorization
* The scaling strategy as presented here preserves the symmetry of the problem
* NOTE It is assumed that the matrix does not have empty row or column,
*
* Example with key steps
* \code
* VectorXd x(n), b(n);
* SparseMatrix<double> A;
* // fill A and b;
* IterScaling<SparseMatrix<double> > scal;
* // Compute the left and right scaling vectors. The matrix is equilibrated at output
* scal.computeRef(A);
* // Scale the right hand side
* b = scal.LeftScaling().cwiseProduct(b);
* // Now, solve the equilibrated linear system with any available solver
*
* // Scale back the computed solution
* x = scal.RightScaling().cwiseProduct(x);
* \endcode
*
* \tparam MatrixType_ the type of the matrix. It should be a real square sparsematrix
*
* References : D. Ruiz and B. Ucar, A Symmetry Preserving Algorithm for Matrix Scaling, INRIA Research report RR-7552
*
* \sa \ref IncompleteLUT
*/
template<typename MatrixType_>
class IterScaling
{
public:
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index;
public:
IterScaling() { init(); }
IterScaling(const MatrixType& matrix)
{
init();
compute(matrix);
}
~IterScaling() { }
/**
* Compute the left and right diagonal matrices to scale the input matrix @p mat
*
* FIXME This algorithm will be modified such that the diagonal elements are permuted on the diagonal.
*
* \sa LeftScaling() RightScaling()
*/
void compute (const MatrixType& mat)
{
using std::abs;
int m = mat.rows();
int n = mat.cols();
eigen_assert((m>0 && m == n) && "Please give a non - empty matrix");
m_left.resize(m);
m_right.resize(n);
m_left.setOnes();
m_right.setOnes();
m_matrix = mat;
VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
Dr.resize(m); Dc.resize(n);
DrRes.resize(m); DcRes.resize(n);
double EpsRow = 1.0, EpsCol = 1.0;
int its = 0;
do
{ // Iterate until the infinite norm of each row and column is approximately 1
// Get the maximum value in each row and column
Dr.setZero(); Dc.setZero();
for (int k=0; k<m_matrix.outerSize(); ++k)
{
for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
{
if ( Dr(it.row()) < abs(it.value()) )
Dr(it.row()) = abs(it.value());
if ( Dc(it.col()) < abs(it.value()) )
Dc(it.col()) = abs(it.value());
}
* \ingroup IterativeSolvers_Module
* \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
*
* This class can be used as a preprocessing tool to accelerate the convergence of iterative methods
*
* This feature is useful to limit the pivoting amount during LU/ILU factorization
* The scaling strategy as presented here preserves the symmetry of the problem
* NOTE It is assumed that the matrix does not have empty row or column,
*
* Example with key steps
* \code
* VectorXd x(n), b(n);
* SparseMatrix<double> A;
* // fill A and b;
* IterScaling<SparseMatrix<double> > scal;
* // Compute the left and right scaling vectors. The matrix is equilibrated at output
* scal.computeRef(A);
* // Scale the right hand side
* b = scal.LeftScaling().cwiseProduct(b);
* // Now, solve the equilibrated linear system with any available solver
*
* // Scale back the computed solution
* x = scal.RightScaling().cwiseProduct(x);
* \endcode
*
* \tparam MatrixType_ the type of the matrix. It should be a real square sparsematrix
*
* References : D. Ruiz and B. Ucar, A Symmetry Preserving Algorithm for Matrix Scaling, INRIA Research report RR-7552
*
* \sa \ref IncompleteLUT
*/
template <typename MatrixType_>
class IterScaling {
public:
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index;
public:
IterScaling() { init(); }
IterScaling(const MatrixType& matrix) {
init();
compute(matrix);
}
~IterScaling() {}
/**
* Compute the left and right diagonal matrices to scale the input matrix @p mat
*
* FIXME This algorithm will be modified such that the diagonal elements are permuted on the diagonal.
*
* \sa LeftScaling() RightScaling()
*/
void compute(const MatrixType& mat) {
using std::abs;
int m = mat.rows();
int n = mat.cols();
eigen_assert((m > 0 && m == n) && "Please give a non - empty matrix");
m_left.resize(m);
m_right.resize(n);
m_left.setOnes();
m_right.setOnes();
m_matrix = mat;
VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
Dr.resize(m);
Dc.resize(n);
DrRes.resize(m);
DcRes.resize(n);
double EpsRow = 1.0, EpsCol = 1.0;
int its = 0;
do { // Iterate until the infinite norm of each row and column is approximately 1
// Get the maximum value in each row and column
Dr.setZero();
Dc.setZero();
for (int k = 0; k < m_matrix.outerSize(); ++k) {
for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) {
if (Dr(it.row()) < abs(it.value())) Dr(it.row()) = abs(it.value());
if (Dc(it.col()) < abs(it.value())) Dc(it.col()) = abs(it.value());
}
for (int i = 0; i < m; ++i)
{
Dr(i) = std::sqrt(Dr(i));
}
for (int i = 0; i < m; ++i) {
Dr(i) = std::sqrt(Dr(i));
}
for (int i = 0; i < n; ++i) {
Dc(i) = std::sqrt(Dc(i));
}
// Save the scaling factors
for (int i = 0; i < m; ++i) {
m_left(i) /= Dr(i);
}
for (int i = 0; i < n; ++i) {
m_right(i) /= Dc(i);
}
// Scale the rows and the columns of the matrix
DrRes.setZero();
DcRes.setZero();
for (int k = 0; k < m_matrix.outerSize(); ++k) {
for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) {
it.valueRef() = it.value() / (Dr(it.row()) * Dc(it.col()));
// Accumulate the norms of the row and column vectors
if (DrRes(it.row()) < abs(it.value())) DrRes(it.row()) = abs(it.value());
if (DcRes(it.col()) < abs(it.value())) DcRes(it.col()) = abs(it.value());
}
for (int i = 0; i < n; ++i)
{
Dc(i) = std::sqrt(Dc(i));
}
// Save the scaling factors
for (int i = 0; i < m; ++i)
{
m_left(i) /= Dr(i);
}
for (int i = 0; i < n; ++i)
{
m_right(i) /= Dc(i);
}
// Scale the rows and the columns of the matrix
DrRes.setZero(); DcRes.setZero();
for (int k=0; k<m_matrix.outerSize(); ++k)
{
for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
{
it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
// Accumulate the norms of the row and column vectors
if ( DrRes(it.row()) < abs(it.value()) )
DrRes(it.row()) = abs(it.value());
if ( DcRes(it.col()) < abs(it.value()) )
DcRes(it.col()) = abs(it.value());
}
}
DrRes.array() = (1-DrRes.array()).abs();
EpsRow = DrRes.maxCoeff();
DcRes.array() = (1-DcRes.array()).abs();
EpsCol = DcRes.maxCoeff();
its++;
}while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
m_isInitialized = true;
}
/** Compute the left and right vectors to scale the vectors
* the input matrix is scaled with the computed vectors at output
*
* \sa compute()
*/
void computeRef (MatrixType& mat)
{
compute (mat);
mat = m_matrix;
}
/** Get the vector to scale the rows of the matrix
*/
VectorXd& LeftScaling()
{
return m_left;
}
/** Get the vector to scale the columns of the matrix
*/
VectorXd& RightScaling()
{
return m_right;
}
/** Set the tolerance for the convergence of the iterative scaling algorithm
*/
void setTolerance(double tol)
{
m_tol = tol;
}
protected:
void init()
{
m_tol = 1e-10;
m_maxits = 5;
m_isInitialized = false;
}
MatrixType m_matrix;
mutable ComputationInfo m_info;
bool m_isInitialized;
VectorXd m_left; // Left scaling vector
VectorXd m_right; // m_right scaling vector
double m_tol;
int m_maxits; // Maximum number of iterations allowed
}
DrRes.array() = (1 - DrRes.array()).abs();
EpsRow = DrRes.maxCoeff();
DcRes.array() = (1 - DcRes.array()).abs();
EpsCol = DcRes.maxCoeff();
its++;
} while ((EpsRow > m_tol || EpsCol > m_tol) && (its < m_maxits));
m_isInitialized = true;
}
/** Compute the left and right vectors to scale the vectors
* the input matrix is scaled with the computed vectors at output
*
* \sa compute()
*/
void computeRef(MatrixType& mat) {
compute(mat);
mat = m_matrix;
}
/** Get the vector to scale the rows of the matrix
*/
VectorXd& LeftScaling() { return m_left; }
/** Get the vector to scale the columns of the matrix
*/
VectorXd& RightScaling() { return m_right; }
/** Set the tolerance for the convergence of the iterative scaling algorithm
*/
void setTolerance(double tol) { m_tol = tol; }
protected:
void init() {
m_tol = 1e-10;
m_maxits = 5;
m_isInitialized = false;
}
MatrixType m_matrix;
mutable ComputationInfo m_info;
bool m_isInitialized;
VectorXd m_left; // Left scaling vector
VectorXd m_right; // m_right scaling vector
double m_tol;
int m_maxits; // Maximum number of iterations allowed
};
}
} // namespace Eigen
#endif