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

merge default branch

Browse Source
This commit is contained in:
Gael Guennebaud
2014-08-29 15:20:31 +02:00
parent f29dbec321 aec3d90ca6
commit 124d12a915
41 changed files with 630 additions and 296 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
unsupported/Eigen/MPRealSupport
View File

@@ -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()));

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

@@ -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;

8
unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h
View File

@@ -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. ;
}

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

@@ -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;

12
unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
View File

@@ -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;

9
unsupported/Eigen/src/SparseExtra/MarketIO.h
View File

@@ -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) )
{

4
unsupported/doc/examples/CMakeLists.txt
View File

@@ -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})

4
unsupported/doc/snippets/CMakeLists.txt
View File

@@ -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})

4
unsupported/test/mpreal/mpreal.h
View File

@@ -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__ */
#endif /* __MPREAL_H__ */