diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-08-29 15:20:31 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-08-29 15:20:31 +0200 |
commit | 124d12a915129bc36ebe87f483712505a11dc91f (patch) | |
tree | 5c0b12148e55cfbfa2c2e69368d982774d96193f /unsupported | |
parent | f29dbec321617d46287c4415889c4485ad70bea3 (diff) | |
parent | aec3d90ca65528fdface6013ccbcc33b04ada867 (diff) |
merge default branch
Diffstat (limited to 'unsupported')
-rw-r--r-- | unsupported/Eigen/MPRealSupport | 10 | ||||
-rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/GMRES.h | 156 | ||||
-rw-r--r-- | unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h | 8 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h | 4 | ||||
-rw-r--r-- | unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h | 12 | ||||
-rw-r--r-- | unsupported/Eigen/src/SparseExtra/MarketIO.h | 9 | ||||
-rw-r--r-- | unsupported/doc/examples/CMakeLists.txt | 4 | ||||
-rw-r--r-- | unsupported/doc/snippets/CMakeLists.txt | 4 | ||||
-rw-r--r-- | unsupported/test/mpreal/mpreal.h | 4 |
9 files changed, 110 insertions, 101 deletions
diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport index 0584e470e..630e0aa77 100644 --- a/unsupported/Eigen/MPRealSupport +++ b/unsupported/Eigen/MPRealSupport @@ -88,9 +88,9 @@ int main() inline static Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); } inline static Real dummy_precision() - { - unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; - return mpfr::machine_epsilon(weak_prec); + { + mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; + return mpfr::machine_epsilon(weak_prec); } }; @@ -159,7 +159,11 @@ int main() { if(rows==0 || cols==0 || depth==0) return; +<<<<<<< local +======= + +>>>>>>> other mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())), tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr())); diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index c8c84069e..67498705b 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -11,7 +11,7 @@ #ifndef EIGEN_GMRES_H #define EIGEN_GMRES_H -namespace Eigen { +namespace Eigen { namespace internal { @@ -27,11 +27,11 @@ namespace internal { * \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: residual tolerance + * \param tol_error on input: relative residual tolerance * on output: residuum achieved * - * \sa IterativeMethods::bicgstab() - * + * \sa IterativeMethods::bicgstab() + * * * For references, please see: * @@ -70,18 +70,24 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition const int m = mat.rows(); - VectorType p0 = rhs - mat*x; + // residual and preconditioned residual + const VectorType p0 = rhs - mat*x; VectorType r0 = precond.solve(p0); - + + const RealScalar r0Norm = r0.norm(); + // is initial guess already good enough? - if(abs(r0.norm()) < tol) { - return true; + 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); - - FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix VectorType tau = VectorType::Zero(restart + 1); + + // storage for Jacobi rotations std::vector < JacobiRotation < Scalar > > G(restart); // generate first Householder vector @@ -112,11 +118,10 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition } if (v.tail(m - k).norm() != 0.0) { - if (k <= restart) { // generate new Householder vector - VectorType e(m - k - 1); + VectorType e(m - k - 1); RealScalar beta; v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta); H.col(k).tail(m - k - 1) = e; @@ -125,78 +130,77 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data()); } - } + } - if (k > 1) { - for (int i = 0; i < k - 1; ++i) { - // apply old Givens rotations to v - v.applyOnTheLeft(i, i + 1, G[i].adjoint()); - } - } + if (k > 1) { + for (int 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) { - // determine next Givens rotation - G[k - 1].makeGivens(v(k - 1), v(k)); + if (k<m && v(k) != (Scalar) 0) { - // apply Givens rotation to v and w - v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); - w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); + // determine next Givens rotation + G[k - 1].makeGivens(v(k - 1), v(k)); - } + // apply Givens rotation to v and w + v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); + w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); + } - // insert coefficients into upper matrix triangle - H.col(k - 1).head(k) = v.head(k); + // insert coefficients into upper matrix triangle + H.col(k - 1).head(k) = v.head(k); - bool stop=(k==m || abs(w(k)) < tol || iters == maxIters); + bool stop=(k==m || abs(w(k)) < tol * r0Norm || iters == maxIters); - if (stop || k == restart) { + if (stop || k == restart) { - // solve upper triangular system - VectorType y = w.head(k); - H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y); + // solve upper triangular system + VectorType y = w.head(k); + H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y); - // use Horner-like scheme to calculate solution vector - VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); + // use Horner-like scheme to calculate solution vector + VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); - // apply Householder reflection H_{k} to x_new - x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); + // apply Householder reflection H_{k} to x_new + x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); - for (int i = k - 2; i >= 0; --i) { - x_new += y(i) * VectorType::Unit(m, 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()); - } + for (int i = k - 2; i >= 0; --i) { + x_new += y(i) * VectorType::Unit(m, 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()); + } - x += x_new; + x += x_new; - if (stop) { - return true; - } else { - k=0; + if (stop) { + return true; + } else { + k=0; - // reset data for a restart r0 = rhs - mat * x; - VectorType p0=mat*x; - VectorType p1=precond.solve(p0); - r0 = rhs - p1; -// r0_sqnorm = r0.squaredNorm(); - w = VectorType::Zero(restart + 1); - H = FMatrixType::Zero(m, restart + 1); - tau = VectorType::Zero(restart + 1); + // reset data for restart + const VectorType p0 = rhs - mat*x; + r0 = precond.solve(p0); - // generate first Householder vector - RealScalar beta; - r0.makeHouseholder(e, tau.coeffRef(0), beta); - w(0)=(Scalar) beta; - H.bottomLeftCorner(m - 1, 1) = e; + // clear Hessenberg matrix and Householder data + H = FMatrixType::Zero(m, restart + 1); + w = VectorType::Zero(restart + 1); + tau = VectorType::Zero(restart + 1); - } + // generate first Householder vector + RealScalar beta; + r0.makeHouseholder(e, tau.coeffRef(0), beta); + w(0)=(Scalar) beta; + H.bottomLeftCorner(m - 1, 1) = e; - } + } + } } - + return false; } @@ -230,7 +234,7 @@ struct traits<GMRES<_MatrixType,_Preconditioner> > * 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; @@ -244,7 +248,7 @@ struct traits<GMRES<_MatrixType,_Preconditioner> > * // 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. Here is a step by * step execution example starting with a random guess and printing the evolution @@ -260,7 +264,7 @@ struct traits<GMRES<_MatrixType,_Preconditioner> > * } while (solver.info()!=Success && i<100); * \endcode * Note that such a step by step excution is slightly slower. - * + * * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner */ template< typename _MatrixType, typename _Preconditioner> @@ -272,10 +276,10 @@ class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> > using Base::m_iterations; using Base::m_info; using Base::m_isInitialized; - + private: int m_restart; - + public: typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; @@ -289,10 +293,10 @@ public: 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 @@ -301,16 +305,16 @@ public: GMRES(const MatrixType& A) : Base(A), m_restart(30) {} ~GMRES() {} - + /** Get the number of iterations after that a restart is performed. */ int 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 int restart) { m_restart=restart; } - + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A * \a x0 as an initial solution. * @@ -326,17 +330,17 @@ public: return internal::solve_retval_with_guess <GMRES, Rhs, Guess>(*this, b.derived(), x0); } - + /** \internal */ template<typename Rhs,typename Dest> void _solveWithGuess(const Rhs& b, Dest& x) const - { + { bool failed = false; for(int j=0; j<b.cols(); ++j) { m_iterations = Base::maxIterations(); 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)) failed = true; diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h index 51dd1d3c4..7cebe4e06 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h @@ -144,11 +144,13 @@ class LevenbergMarquardt : internal::no_assignment_operator /** Sets the default parameters */ void resetParameters() - { + { + using std::sqrt; + m_factor = 100.; m_maxfev = 400; - m_ftol = std::sqrt(NumTraits<RealScalar>::epsilon()); - m_xtol = std::sqrt(NumTraits<RealScalar>::epsilon()); + m_ftol = sqrt(NumTraits<RealScalar>::epsilon()); + m_xtol = sqrt(NumTraits<RealScalar>::epsilon()); m_gtol = 0. ; m_epsfcn = 0. ; } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 05df5a102..160120d03 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -176,8 +176,8 @@ void matrix_exp_pade17(const MatrixType &A, MatrixType &U, MatrixType &V) const MatrixType A4 = A2 * A2; const MatrixType A6 = A4 * A2; const MatrixType A8 = A4 * A4; - V = b[17] * m_tmp1 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage - matrixType tmp = A8 * V; + V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage + MatrixType tmp = A8 * V; tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); U.noalias() = A * tmp; diff --git a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h index bfeb26fc9..ecb8dccf4 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h @@ -55,8 +55,8 @@ public: Parameters() : factor(Scalar(100.)) , maxfev(400) - , ftol(std::sqrt(NumTraits<Scalar>::epsilon())) - , xtol(std::sqrt(NumTraits<Scalar>::epsilon())) + , ftol(sqrt_(NumTraits<Scalar>::epsilon())) + , xtol(sqrt_(NumTraits<Scalar>::epsilon())) , gtol(Scalar(0.)) , epsfcn(Scalar(0.)) {} Scalar factor; @@ -72,7 +72,7 @@ public: LevenbergMarquardtSpace::Status lmder1( FVectorType &x, - const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) + const Scalar tol = sqrt_(NumTraits<Scalar>::epsilon()) ); LevenbergMarquardtSpace::Status minimize(FVectorType &x); @@ -83,12 +83,12 @@ public: FunctorType &functor, FVectorType &x, Index *nfev, - const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) + const Scalar tol = sqrt_(NumTraits<Scalar>::epsilon()) ); LevenbergMarquardtSpace::Status lmstr1( FVectorType &x, - const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) + const Scalar tol = sqrt_(NumTraits<Scalar>::epsilon()) ); LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x); @@ -109,6 +109,8 @@ public: Scalar lm_param(void) { return par; } private: + static Scalar sqrt_(const Scalar& x) { using std::sqrt; return sqrt(x); } + FunctorType &functor; Index n; Index m; diff --git a/unsupported/Eigen/src/SparseExtra/MarketIO.h b/unsupported/Eigen/src/SparseExtra/MarketIO.h index 1c40d3f7c..25ff4228d 100644 --- a/unsupported/Eigen/src/SparseExtra/MarketIO.h +++ b/unsupported/Eigen/src/SparseExtra/MarketIO.h @@ -133,6 +133,7 @@ template<typename SparseMatrixType> bool loadMarket(SparseMatrixType& mat, const std::string& filename) { typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::Index Index; std::ifstream input(filename.c_str(),std::ios::in); if(!input) return false; @@ -142,11 +143,11 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename) bool readsizes = false; - typedef Triplet<Scalar,int> T; + typedef Triplet<Scalar,Index> T; std::vector<T> elements; - int M(-1), N(-1), NNZ(-1); - int count = 0; + Index M(-1), N(-1), NNZ(-1); + Index count = 0; while(input.getline(buffer, maxBuffersize)) { // skip comments @@ -169,7 +170,7 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename) } else { - int i(-1), j(-1); + Index i(-1), j(-1); Scalar value; if( internal::GetMarketLine(line, M, N, i, j, value) ) { diff --git a/unsupported/doc/examples/CMakeLists.txt b/unsupported/doc/examples/CMakeLists.txt index 978f9afd8..c47646dfc 100644 --- a/unsupported/doc/examples/CMakeLists.txt +++ b/unsupported/doc/examples/CMakeLists.txt @@ -10,12 +10,10 @@ FOREACH(example_src ${examples_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(example_${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - GET_TARGET_PROPERTY(example_executable - example_${example} LOCATION) ADD_CUSTOM_COMMAND( TARGET example_${example} POST_BUILD - COMMAND ${example_executable} + COMMAND example_${example} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out ) ADD_DEPENDENCIES(unsupported_examples example_${example}) diff --git a/unsupported/doc/snippets/CMakeLists.txt b/unsupported/doc/snippets/CMakeLists.txt index 4a4157933..f0c5cc2a8 100644 --- a/unsupported/doc/snippets/CMakeLists.txt +++ b/unsupported/doc/snippets/CMakeLists.txt @@ -14,12 +14,10 @@ FOREACH(snippet_src ${snippets_SRCS}) if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO) target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}) endif() - GET_TARGET_PROPERTY(compile_snippet_executable - ${compile_snippet_target} LOCATION) ADD_CUSTOM_COMMAND( TARGET ${compile_snippet_target} POST_BUILD - COMMAND ${compile_snippet_executable} + COMMAND ${compile_snippet_target} ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out ) ADD_DEPENDENCIES(unsupported_snippets ${compile_snippet_target}) diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index 4c29d2ac2..dddda7dd9 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -1698,7 +1698,7 @@ inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcpt //////////////////////////////////////////////////////////////////////////
// Type Converters
-inline bool mpreal::toBool (mp_rnd_t mode) const { return mpfr_zero_p (mpfr_srcptr()) == 0; }
+inline bool mpreal::toBool (mp_rnd_t /*mode*/) const { return mpfr_zero_p (mpfr_srcptr()) == 0; }
inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si (mpfr_srcptr(), mode); }
inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui (mpfr_srcptr(), mode); }
inline float mpreal::toFloat (mp_rnd_t mode) const { return mpfr_get_flt(mpfr_srcptr(), mode); }
@@ -3070,4 +3070,4 @@ namespace std }
-#endif /* __MPREAL_H__ */
\ No newline at end of file +#endif /* __MPREAL_H__ */
|