diff options
Diffstat (limited to 'eigen/unsupported/Eigen/src')
3 files changed, 55 insertions, 47 deletions
diff --git a/eigen/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/eigen/unsupported/Eigen/src/IterativeSolvers/DGMRES.h index bae04fc..eff2dc8 100644 --- a/eigen/unsupported/Eigen/src/IterativeSolvers/DGMRES.h +++ b/eigen/unsupported/Eigen/src/IterativeSolvers/DGMRES.h @@ -39,7 +39,6 @@ template <typename VectorType, typename IndexType> void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut) { eigen_assert(vec.size() == perm.size()); - typedef typename IndexType::Scalar Index; bool flag; for (Index k = 0; k < ncut; k++) { @@ -112,7 +111,6 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > using Base::_solve_impl; 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; @@ -146,7 +144,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > void _solve_with_guess_impl(const Rhs& b, Dest& x) const { bool failed = false; - for(int j=0; j<b.cols(); ++j) + for(Index j=0; j<b.cols(); ++j) { m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; @@ -170,17 +168,17 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > /** * Get the restart value */ - int restart() { return m_restart; } + Index restart() { return m_restart; } /** * Set the restart value (default is 30) */ - void set_restart(const int restart) { m_restart=restart; } + Index set_restart(const Index restart) { m_restart=restart; } /** * Set the number of eigenvalues to deflate at each restart */ - void setEigenv(const int neig) + void setEigenv(const Index neig) { m_neig = neig; if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates @@ -189,12 +187,12 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > /** * Get the size of the deflation subspace size */ - int deflSize() {return m_r; } + Index deflSize() {return m_r; } /** * Set the maximum size of the deflation subspace */ - void setMaxEigenv(const int maxNeig) { m_maxNeig = maxNeig; } + void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; } protected: // DGMRES algorithm @@ -202,12 +200,12 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const; // Perform one cycle of GMRES template<typename Dest> - int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const; + 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 - int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const; + Index 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; + 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 @@ -221,8 +219,8 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > 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 int m_r; // Current number of deflated eigenvalues, size of m_U - mutable int m_maxNeig; // Maximum number of eigenvalues to deflate + 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; @@ -244,9 +242,9 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh const Preconditioner& precond) const { //Initialization - int n = mat.rows(); + Index n = mat.rows(); DenseVector r0(n); - int nbIts = 0; + 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); @@ -284,7 +282,7 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh */ template< typename _MatrixType, typename _Preconditioner> template<typename Dest> -int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const +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 @@ -293,8 +291,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con m_V.col(0) = r0/beta; m_info = NoConvergence; std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations - int it = 0; // Number of inner iterations - int n = mat.rows(); + 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) { @@ -312,7 +310,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt Scalar coef; - for (int i = 0; i <= it; ++i) + for (Index i = 0; i <= it; ++i) { coef = tv1.dot(m_V.col(i)); tv1 = tv1 - coef * m_V.col(i); @@ -328,7 +326,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con // FIXME Check for happy breakdown // Update Hessenberg matrix with Givens rotations - for (int i = 1; i <= it; ++i) + for (Index i = 1; i <= it; ++i) { m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint()); } @@ -418,7 +416,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, StorageIndex& neig) const +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 typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH; @@ -433,8 +431,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri // Reorder the absolute values of Schur values DenseRealVector modulEig(it); - for (int j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); - perm.setLinSpaced(it,0,it-1); + 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) @@ -442,7 +440,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN); } //Count the real number of extracted eigenvalues (with complex conjugates) - int nbrEig = 0; + Index nbrEig = 0; while (nbrEig < neig) { if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; @@ -451,7 +449,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri // Extract the Schur vectors corresponding to the smallest Ritz values DenseMatrix Sr(it, nbrEig); Sr.setZero(); - for (int j = 0; j < nbrEig; j++) + for (Index j = 0; j < nbrEig; j++) { Sr.col(j) = schurofH.matrixU().col(perm(it-j-1)); } @@ -462,8 +460,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri if (m_r) { // Orthogonalize X against m_U using modified Gram-Schmidt - for (int j = 0; j < nbrEig; j++) - for (int k =0; k < m_r; k++) + 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); } @@ -473,7 +471,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri dgmresInitDeflation(m); DenseMatrix MX(m, nbrEig); DenseVector tv1(m); - for (int j = 0; j < nbrEig; j++) + for (Index j = 0; j < nbrEig; j++) { tv1 = mat * X.col(j); MX.col(j) = precond.solve(tv1); @@ -488,8 +486,8 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri } // Save X into m_U and m_MX in m_MU - for (int j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j); - for (int 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; @@ -502,7 +500,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri } template<typename _MatrixType, typename _Preconditioner> template<typename RhsType, typename DestType> -int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const +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); diff --git a/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index bb6d9e1..85ab3d9 100644 --- a/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -326,6 +326,7 @@ struct matrix_exp_computeUV<MatrixType, long double> } else if (l1norm < 1.125358383453143065081397882891878e+000L) { matrix_exp_pade13(arg, U, V); } else { + const long double maxnorm = 2.884233277829519311757165057717815L; frexp(l1norm / maxnorm, &squarings); if (squarings < 0) squarings = 0; MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings)); @@ -342,6 +343,27 @@ struct matrix_exp_computeUV<MatrixType, long double> } }; +template<typename T> struct is_exp_known_type : false_type {}; +template<> struct is_exp_known_type<float> : true_type {}; +template<> struct is_exp_known_type<double> : true_type {}; +#if LDBL_MANT_DIG <= 112 +template<> struct is_exp_known_type<long double> : true_type {}; +#endif + +template <typename ArgType, typename ResultType> +void matrix_exp_compute(const ArgType& arg, ResultType &result, true_type) // natively supported scalar type +{ + typedef typename ArgType::PlainObject MatrixType; + MatrixType U, V; + int squarings; + matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V) + MatrixType numer = U + V; + MatrixType denom = -U + V; + result = denom.partialPivLu().solve(numer); + for (int i=0; i<squarings; i++) + result *= result; // undo scaling by repeated squaring +} + /* Computes the matrix exponential * @@ -349,26 +371,13 @@ struct matrix_exp_computeUV<MatrixType, long double> * \param result variable in which result will be stored */ template <typename ArgType, typename ResultType> -void matrix_exp_compute(const ArgType& arg, ResultType &result) +void matrix_exp_compute(const ArgType& arg, ResultType &result, false_type) // default { typedef typename ArgType::PlainObject MatrixType; -#if LDBL_MANT_DIG > 112 // rarely happens typedef typename traits<MatrixType>::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename std::complex<RealScalar> ComplexScalar; - if (sizeof(RealScalar) > 14) { - result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>); - return; - } -#endif - MatrixType U, V; - int squarings; - matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V) - MatrixType numer = U + V; - MatrixType denom = -U + V; - result = denom.partialPivLu().solve(numer); - for (int i=0; i<squarings; i++) - result *= result; // undo scaling by repeated squaring + result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>); } } // end namespace Eigen::internal @@ -402,7 +411,7 @@ template<typename Derived> struct MatrixExponentialReturnValue inline void evalTo(ResultType& result) const { const typename internal::nested_eval<Derived, 10>::type tmp(m_src); - internal::matrix_exp_compute(tmp, result); + internal::matrix_exp_compute(tmp, result, internal::is_exp_known_type<typename Derived::Scalar>()); } Index rows() const { return m_src.rows(); } diff --git a/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h b/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h index cdc14f8..41e4af4 100644 --- a/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -109,6 +109,7 @@ namespace internal inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector) { sym = 0; + iscomplex = false; isvector = false; std::ifstream in(filename.c_str(),std::ios::in); if(!in) |