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

Merged eigen/eigen into default

Browse Source
This commit is contained in:
Tal Hadad
2015-12-20 12:50:07 +02:00
parent 6752a69aa5 6d777e1bc7
commit fabd8474ff
296 changed files with 10487 additions and 4148 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h Normal file → Executable file
View File

@@ -99,7 +99,11 @@ class AutoDiffScalar
{}
template<typename OtherDerType>
AutoDiffScalar(const AutoDiffScalar<OtherDerType>& other)
AutoDiffScalar(const AutoDiffScalar<OtherDerType>& other
#ifndef EIGEN_PARSED_BY_DOXYGEN
, typename internal::enable_if<internal::is_same<Scalar,typename OtherDerType::Scalar>::value,void*>::type = 0
#endif
)
: m_value(other.value()), m_derivatives(other.derivatives())
{}
@@ -127,6 +131,14 @@ class AutoDiffScalar
return *this;
}
inline AutoDiffScalar& operator=(const Scalar& other)
{
m_value = other;
if(m_derivatives.size()>0)
m_derivatives.setZero();
return *this;
}
// inline operator const Scalar& () const { return m_value; }
// inline operator Scalar& () { return m_value; }
@@ -626,9 +638,10 @@ EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(acos,
template<typename DerType> struct NumTraits<AutoDiffScalar<DerType> >
: NumTraits< typename NumTraits<typename DerType::Scalar>::Real >
{
typedef AutoDiffScalar<Matrix<typename NumTraits<typename DerType::Scalar>::Real,DerType::RowsAtCompileTime,DerType::ColsAtCompileTime> > Real;
typedef AutoDiffScalar<Matrix<typename NumTraits<typename DerType::Scalar>::Real,DerType::RowsAtCompileTime,DerType::ColsAtCompileTime,
DerType::Options, DerType::MaxRowsAtCompileTime, DerType::MaxColsAtCompileTime> > Real;
typedef AutoDiffScalar<DerType> NonInteger;
typedef AutoDiffScalar<DerType>& Nested;
typedef AutoDiffScalar<DerType> Nested;
enum{
RequireInitialization = 1
};

1
unsupported/Eigen/src/CMakeLists.txt
View File

@@ -1,5 +1,6 @@
ADD_SUBDIRECTORY(AutoDiff)
ADD_SUBDIRECTORY(BVH)
ADD_SUBDIRECTORY(Eigenvalues)
ADD_SUBDIRECTORY(FFT)
ADD_SUBDIRECTORY(IterativeSolvers)
ADD_SUBDIRECTORY(LevenbergMarquardt)

6
unsupported/Eigen/src/Eigenvalues/CMakeLists.txt Normal file
View File

@@ -0,0 +1,6 @@
FILE(GLOB Eigen_Eigenvalues_SRCS "*.h")
INSTALL(FILES
${Eigen_Eigenvalues_SRCS}
DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Eigenvalues COMPONENT Devel
)

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

@@ -40,7 +40,6 @@ void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::
{
eigen_assert(vec.size() == perm.size());
typedef typename IndexType::Scalar Index;
typedef typename VectorType::Scalar Scalar;
bool flag;
for (Index k = 0; k < ncut; k++)
{
@@ -84,6 +83,8 @@ void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::
* 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
@@ -101,7 +102,7 @@ template< typename _MatrixType, typename _Preconditioner>
class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
{
typedef IterativeSolverBase<DGMRES> Base;
using Base::mp_matrix;
using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
@@ -112,6 +113,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index;
typedef typename MatrixType::StorageIndex StorageIndex;
typedef typename MatrixType::RealScalar RealScalar;
typedef _Preconditioner Preconditioner;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
@@ -134,8 +136,8 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
DGMRES(const MatrixType& A) : Base(A),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false)
{}
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() {}
@@ -150,7 +152,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
m_error = Base::m_tolerance;
typename Dest::ColXpr xj(x,j);
dgmres(mp_matrix, b.col(j), xj, Base::m_preconditioner);
dgmres(matrix(), b.col(j), xj, Base::m_preconditioner);
}
m_info = failed ? NumericalIssue
: m_error <= Base::m_tolerance ? Success
@@ -202,7 +204,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
template<typename Dest>
int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const;
// Compute data to use for deflation
int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const;
int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const;
// Apply deflation to a vector
template<typename RhsType, typename DestType>
int dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
@@ -218,7 +220,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
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 int m_neig; //Number of eigenvalues to extract at each restart
mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart
mutable int m_r; // Current number of deflated eigenvalues, size of m_U
mutable int m_maxNeig; // Maximum number of eigenvalues to deflate
mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
@@ -338,7 +340,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con
beta = std::abs(g(it+1));
m_error = beta/normRhs;
std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
// std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
it++; nbIts++;
if (m_error < m_tolerance)
@@ -416,7 +418,7 @@ inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_Matr
}
template< typename _MatrixType, typename _Preconditioner>
int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const
int 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
typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH;
@@ -426,7 +428,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri
schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU);
ComplexVector eig(it);
Matrix<Index,Dynamic,1>perm(it);
Matrix<StorageIndex,Dynamic,1>perm(it);
eig = this->schurValues(schurofH);
// Reorder the absolute values of Schur values

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

@@ -251,13 +251,15 @@ struct traits<GMRES<_MatrixType,_Preconditioner> >
* 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::mp_matrix;
using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
@@ -288,7 +290,8 @@ public:
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
GMRES(const MatrixType& A) : Base(A), m_restart(30) {}
template<typename MatrixDerived>
explicit GMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30) {}
~GMRES() {}
@@ -312,7 +315,7 @@ public:
m_error = Base::m_tolerance;
typename Dest::ColXpr xj(x,j);
if(!internal::gmres(mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
if(!internal::gmres(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
failed = true;
}
m_info = failed ? NumericalIssue

305
unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
View File

@@ -1,305 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// 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_INCOMPLETE_CHOlESKY_H
#define EIGEN_INCOMPLETE_CHOlESKY_H
#include "Eigen/src/IterativeLinearSolvers/IncompleteLUT.h"
#include <Eigen/OrderingMethods>
#include <list>
namespace Eigen {
/**
* \brief Modified Incomplete Cholesky with dual threshold
*
* References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
* Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999
*
* \tparam _MatrixType The type of the sparse matrix. It should be a symmetric
* matrix. It is advised to give a row-oriented sparse matrix
* \tparam _UpLo The triangular part of the matrix to reference.
* \tparam _OrderingType
*/
template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> >
class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
{
protected:
typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
using Base::m_isInitialized;
public:
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef _OrderingType OrderingType;
typedef typename OrderingType::PermutationType PermutationType;
typedef typename PermutationType::StorageIndex StorageIndex;
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
typedef FactorType MatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorSx;
typedef Matrix<RealScalar,Dynamic,1> VectorRx;
typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
typedef std::vector<std::list<StorageIndex> > VectorList;
enum { UpLo = _UpLo };
public:
IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {}
template<typename MatrixType>
IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false)
{
compute(matrix);
}
Index rows() const { return m_L.rows(); }
Index cols() const { return m_L.cols(); }
/** \brief Reports whether previous computation was successful.
*
* \returns \c Success if computation was succesful,
* \c NumericalIssue if the matrix appears to be negative.
*/
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
return m_info;
}
/**
* \brief Set the initial shift parameter
*/
void setInitialShift(RealScalar shift) { m_initialShift = shift; }
/**
* \brief Computes the fill reducing permutation vector.
*/
template<typename MatrixType>
void analyzePattern(const MatrixType& mat)
{
OrderingType ord;
PermutationType pinv;
ord(mat.template selfadjointView<UpLo>(), pinv);
if(pinv.size()>0) m_perm = pinv.inverse();
else m_perm.resize(0);
m_analysisIsOk = true;
}
template<typename MatrixType>
void factorize(const MatrixType& amat);
template<typename MatrixType>
void compute(const MatrixType& matrix)
{
analyzePattern(matrix);
factorize(matrix);
}
template<typename Rhs, typename Dest>
void _solve_impl(const Rhs& b, Dest& x) const
{
eigen_assert(m_factorizationIsOk && "factorize() should be called first");
if (m_perm.rows() == b.rows()) x = m_perm * b;
else x = b;
x = m_scale.asDiagonal() * x;
x = m_L.template triangularView<Lower>().solve(x);
x = m_L.adjoint().template triangularView<Upper>().solve(x);
x = m_scale.asDiagonal() * x;
if (m_perm.rows() == b.rows())
x = m_perm.inverse() * x;
}
protected:
FactorType m_L; // The lower part stored in CSC
VectorRx m_scale; // The vector for scaling the matrix
RealScalar m_initialShift; // The initial shift parameter
bool m_analysisIsOk;
bool m_factorizationIsOk;
ComputationInfo m_info;
PermutationType m_perm;
private:
inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol);
};
template<typename Scalar, int _UpLo, typename OrderingType>
template<typename _MatrixType>
void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
{
using std::sqrt;
eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
// Dropping strategy : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
m_L.resize(mat.rows(), mat.cols());
// Apply the fill-reducing permutation computed in analyzePattern()
if (m_perm.rows() == mat.rows() ) // To detect the null permutation
{
// The temporary is needed to make sure that the diagonal entry is properly sorted
FactorType tmp(mat.rows(), mat.cols());
tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
}
else
{
m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
}
Index n = m_L.cols();
Index nnz = m_L.nonZeros();
Map<VectorSx> vals(m_L.valuePtr(), nnz); //values
Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices
Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
VectorSx col_vals(n); // Store a nonzero values in each column
VectorIx col_irow(n); // Row indices of nonzero elements in each column
VectorIx col_pattern(n);
col_pattern.fill(-1);
StorageIndex col_nnz;
// Computes the scaling factors
m_scale.resize(n);
m_scale.setZero();
for (Index j = 0; j < n; j++)
for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
{
m_scale(j) += numext::abs2(vals(k));
if(rowIdx[k]!=j)
m_scale(rowIdx[k]) += numext::abs2(vals(k));
}
m_scale = m_scale.cwiseSqrt().cwiseSqrt();
// Scale and compute the shift for the matrix
RealScalar mindiag = NumTraits<RealScalar>::highest();
for (Index j = 0; j < n; j++)
{
for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
vals[k] /= (m_scale(j)*m_scale(rowIdx[k]));
eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
}
RealScalar shift = 0;
if(mindiag <= RealScalar(0.))
shift = m_initialShift - mindiag;
// Apply the shift to the diagonal elements of the matrix
for (Index j = 0; j < n; j++)
vals[colPtr[j]] += shift;
// jki version of the Cholesky factorization
for (Index j=0; j < n; ++j)
{
// Left-looking factorization of the j-th column
// First, load the j-th column into col_vals
Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
col_nnz = 0;
for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
{
StorageIndex l = rowIdx[i];
col_vals(col_nnz) = vals[i];
col_irow(col_nnz) = l;
col_pattern(l) = col_nnz;
col_nnz++;
}
{
typename std::list<StorageIndex>::iterator k;
// Browse all previous columns that will update column j
for(k = listCol[j].begin(); k != listCol[j].end(); k++)
{
Index jk = firstElt(*k); // First element to use in the column
eigen_internal_assert(rowIdx[jk]==j);
Scalar v_j_jk = numext::conj(vals[jk]);
jk += 1;
for (Index i = jk; i < colPtr[*k+1]; i++)
{
StorageIndex l = rowIdx[i];
if(col_pattern[l]<0)
{
col_vals(col_nnz) = vals[i] * v_j_jk;
col_irow[col_nnz] = l;
col_pattern(l) = col_nnz;
col_nnz++;
}
else
col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
}
updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
}
}
// Scale the current column
if(numext::real(diag) <= 0)
{
std::cerr << "\nNegative diagonal during Incomplete factorization at position " << j << " (value = " << diag << ")\n";
m_info = NumericalIssue;
return;
}
RealScalar rdiag = sqrt(numext::real(diag));
vals[colPtr[j]] = rdiag;
for (Index k = 0; k<col_nnz; ++k)
{
Index i = col_irow[k];
//Scale
col_vals(k) /= rdiag;
//Update the remaining diagonals with col_vals
vals[colPtr[i]] -= numext::abs2(col_vals(k));
}
// Select the largest p elements
// p is the original number of elements in the column (without the diagonal)
Index p = colPtr[j+1] - colPtr[j] - 1 ;
Ref<VectorSx> cvals = col_vals.head(col_nnz);
Ref<VectorIx> cirow = col_irow.head(col_nnz);
internal::QuickSplit(cvals,cirow, p);
// Insert the largest p elements in the matrix
Index cpt = 0;
for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
{
vals[i] = col_vals(cpt);
rowIdx[i] = col_irow(cpt);
// restore col_pattern:
col_pattern(col_irow(cpt)) = -1;
cpt++;
}
// Get the first smallest row index and put it after the diagonal element
Index jk = colPtr(j)+1;
updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
}
m_factorizationIsOk = true;
m_isInitialized = true;
m_info = Success;
}
template<typename Scalar, int _UpLo, typename OrderingType>
inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol)
{
if (jk < colPtr(col+1) )
{
Index p = colPtr(col+1) - jk;
Index minpos;
rowIdx.segment(jk,p).minCoeff(&minpos);
minpos += jk;
if (rowIdx(minpos) != rowIdx(jk))
{
//Swap
std::swap(rowIdx(jk),rowIdx(minpos));
std::swap(vals(jk),vals(minpos));
}
firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
}
}
} // end namespace Eigen
#endif

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

@@ -191,6 +191,8 @@ namespace Eigen {
* 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>
@@ -198,7 +200,7 @@ namespace Eigen {
{
typedef IterativeSolverBase<MINRES> Base;
using Base::mp_matrix;
using Base::matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
@@ -227,7 +229,8 @@ namespace Eigen {
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
MINRES(const MatrixType& A) : Base(A) {}
template<typename MatrixDerived>
explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
/** Destructor. */
~MINRES(){}
@@ -236,21 +239,31 @@ namespace Eigen {
template<typename Rhs,typename Dest>
void _solve_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 typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper;
EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
typedef typename internal::conditional<UpLo==(Lower|Upper),
Ref<const MatrixType>&,
SparseSelfAdjointView<const Ref<const MatrixType>, UpLo>
>::type MatrixWrapperType;
RowMajorWrapper,
typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
>::type SelfAdjointWrapper;
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
RowMajorWrapper row_mat(matrix());
for(int j=0; j<b.cols(); ++j)
{
m_iterations = Base::maxIterations();
m_error = Base::m_tolerance;
typename Dest::ColXpr xj(x,j);
internal::minres(MatrixWrapperType(mp_matrix), b.col(j), xj,
internal::minres(SelfAdjointWrapper(row_mat), b.col(j), xj,
Base::m_preconditioner, m_iterations, m_error);
}

2
unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h
View File

@@ -240,7 +240,7 @@ struct traits<KroneckerProductSparse<_Lhs,_Rhs> >
Flags = ((LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
| EvalBeforeNestingBit | EvalBeforeAssigningBit,
CoeffReadCost = Dynamic
CoeffReadCost = HugeCost
};
typedef SparseMatrix<Scalar, 0, StorageIndex> ReturnType;

2
unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
View File

@@ -348,7 +348,7 @@ void matrix_exp_compute(const MatrixType& arg, ResultType &result)
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename std::complex<RealScalar> ComplexScalar;
if (sizeof(RealScalar) > 14) {
result = arg.matrixFunction(StdStemFunctions<ComplexScalar>::exp);
result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
return;
}
#endif

2
unsupported/Eigen/src/Skyline/SkylineProduct.h
View File

@@ -49,7 +49,7 @@ struct internal::traits<SkylineProduct<LhsNested, RhsNested, ProductMode> > {
| EvalBeforeAssigningBit
| EvalBeforeNestingBit,
CoeffReadCost = Dynamic
CoeffReadCost = HugeCost
};
typedef typename internal::conditional<ResultIsSkyline,

2
unsupported/Eigen/src/SparseExtra/DynamicSparseMatrix.h
View File

@@ -56,6 +56,8 @@ template<typename _Scalar, int _Options, typename _StorageIndex>
class DynamicSparseMatrix
: public SparseMatrixBase<DynamicSparseMatrix<_Scalar, _Options, _StorageIndex> >
{
typedef SparseMatrixBase<DynamicSparseMatrix> Base;
using Base::convert_index;
public:
EIGEN_SPARSE_PUBLIC_INTERFACE(DynamicSparseMatrix)
// FIXME: why are these operator already alvailable ???