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

Pulled in the latest changes from the Eigen trunk

Browse Source
This commit is contained in:
Benoit Steiner
2014-08-13 22:25:29 -07:00
parent 916ef48846 e51da9c3a8
commit 16047c8d4a
211 changed files with 2321 additions and 13262 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
unsupported/Eigen/CXX11/Core
View File

@@ -34,6 +34,7 @@
#if __cplusplus <= 199711L
#include "src/Core/util/EmulateCXX11Meta.h"
#else
#include <array>
#include "src/Core/util/CXX11Workarounds.h"
#include "src/Core/util/CXX11Meta.h"
#endif

57
unsupported/Eigen/MPRealSupport
View File

@@ -27,6 +27,8 @@ namespace Eigen {
* via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a>
* library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>.
*
* \warning MPFR C++ is licensed under the GPL.
*
* You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder.
*
* Here is an example:
@@ -86,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);
}
};
@@ -139,64 +141,53 @@ int main()
public:
typedef mpfr::mpreal ResScalar;
enum {
nr = 2, // must be 2 for proper packing...
nr = 1,
mr = 1,
WorkSpaceFactor = nr,
LhsProgress = 1,
RhsProgress = 1
};
};
template<typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,mr,nr,ConjugateLhs,ConjugateRhs>
template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,1,1,ConjugateLhs,ConjugateRhs>
{
typedef mpfr::mpreal mpreal;
EIGEN_DONT_INLINE
void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha,
Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0, mpreal* /*unpackedB*/ = 0)
Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
{
mpreal acc1, acc2, tmp;
if(rows==0 || cols==0 || depth==0)
return;
mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())),
tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr()));
if(strideA==-1) strideA = depth;
if(strideB==-1) strideB = depth;
for(Index j=0; j<cols; j+=nr)
for(Index i=0; i<rows; ++i)
{
Index actual_nr = (std::min<Index>)(nr,cols-j);
mpreal *C1 = res + j*resStride;
mpreal *C2 = res + (j+1)*resStride;
for(Index i=0; i<rows; i++)
for(Index j=0; j<cols; ++j)
{
mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr;
mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA;
mpreal *C1 = res + j*resStride;
const mpreal *A = blockA + i*strideA + offsetA;
const mpreal *B = blockB + j*strideB + offsetB;
acc1 = 0;
acc2 = 0;
for(Index k=0; k<depth; k++)
{
mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd());
mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_srcptr(), B[k].mpfr_srcptr(), mpreal::get_default_rnd());
mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd());
if(actual_nr==2) {
mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[1].mpfr_ptr(), mpreal::get_default_rnd());
mpfr_add(acc2.mpfr_ptr(), acc2.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd());
}
B+=actual_nr;
}
mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(), mpreal::get_default_rnd());
if(actual_nr==2) {
mpfr_mul(acc2.mpfr_ptr(), acc2.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
mpfr_add(C2[i].mpfr_ptr(), C2[i].mpfr_ptr(), acc2.mpfr_ptr(), mpreal::get_default_rnd());
}
mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd());
mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_srcptr(), acc1.mpfr_srcptr(), mpreal::get_default_rnd());
}
}
}
};
} // end namespace internal
}

4
unsupported/Eigen/SVD
View File

@@ -29,10 +29,6 @@
#include "../../Eigen/src/SVD/JacobiSVD_MKL.h"
#endif
#ifdef EIGEN2_SUPPORT
#include "../../Eigen/src/Eigen2Support/SVD.h"
#endif
#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
#endif // EIGEN_SVD_MODULE_H

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

@@ -2,7 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2012 Kolja Brix <brix@igpm.rwth-aaachen.de>
// Copyright (C) 2012, 2014 Kolja Brix <brix@igpm.rwth-aaachen.de>
//
// 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
@@ -72,16 +72,20 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
VectorType p0 = rhs - mat*x;
VectorType r0 = precond.solve(p0);
// RealScalar r0_sqnorm = r0.squaredNorm();
// is initial guess already good enough?
if(abs(r0.norm()) < tol) {
return true;
}
VectorType w = VectorType::Zero(restart + 1);
FMatrixType H = FMatrixType::Zero(m, restart + 1);
FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix
VectorType tau = VectorType::Zero(restart + 1);
std::vector < JacobiRotation < Scalar > > G(restart);
// generate first Householder vector
VectorType e;
VectorType e(m-1);
RealScalar beta;
r0.makeHouseholder(e, tau.coeffRef(0), beta);
w(0)=(Scalar) beta;

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

@@ -9,6 +9,9 @@
#ifndef EIGEN_ITERSCALING_H
#define EIGEN_ITERSCALING_H
namespace Eigen {
/**
* \ingroup IterativeSolvers_Module
* \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
@@ -41,8 +44,6 @@
*
* \sa \ref IncompleteLUT
*/
namespace Eigen {
using std::abs;
template<typename _MatrixType>
class IterScaling
{
@@ -71,6 +72,7 @@ class IterScaling
*/
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");

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;

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

64
unsupported/Eigen/src/Splines/Spline.h
View File

@@ -44,9 +44,15 @@ namespace Eigen
/** \brief The data type used to store knot vectors. */
typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
/** \brief The data type used to store parameter vectors. */
typedef typename SplineTraits<Spline>::ParameterVectorType ParameterVectorType;
/** \brief The data type used to store non-zero basis functions. */
typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
/** \brief The data type used to store the values of the basis function derivatives. */
typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType;
/** \brief The data type representing the spline's control points. */
typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
@@ -203,10 +209,25 @@ namespace Eigen
**/
static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
/**
* \copydoc Spline::basisFunctionDerivatives
* \param degree The degree of the underlying spline
* \param knots The underlying spline's knot vector.
**/
static BasisDerivativeType BasisFunctionDerivatives(
const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots);
private:
KnotVectorType m_knots; /*!< Knot vector. */
ControlPointVectorType m_ctrls; /*!< Control points. */
template <typename DerivativeType>
static void BasisFunctionDerivativesImpl(
const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
const DenseIndex order,
const DenseIndex p,
const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
DerivativeType& N_);
};
template <typename _Scalar, int _Dim, int _Degree>
@@ -345,20 +366,24 @@ namespace Eigen
}
/* --------------------------------------------------------------------------------------------- */
template <typename SplineType, typename DerivativeType>
void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_)
template <typename _Scalar, int _Dim, int _Degree>
template <typename DerivativeType>
void Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivativesImpl(
const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
const DenseIndex order,
const DenseIndex p,
const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
DerivativeType& N_)
{
typedef Spline<_Scalar, _Dim, _Degree> SplineType;
enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
typedef typename SplineTraits<SplineType>::Scalar Scalar;
typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType;
const KnotVectorType& U = spline.knots();
const DenseIndex p = spline.degree();
const DenseIndex span = spline.span(u);
const DenseIndex span = SplineType::Span(u, p, U);
const DenseIndex n = (std::min)(p, order);
@@ -455,8 +480,8 @@ namespace Eigen
typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
{
typename SplineTraits< Spline >::BasisDerivativeType der;
basisFunctionDerivativesImpl(*this, u, order, der);
typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType der;
BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
return der;
}
@@ -465,8 +490,21 @@ namespace Eigen
typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
{
typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der;
basisFunctionDerivativesImpl(*this, u, order, der);
typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType der;
BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
return der;
}
template <typename _Scalar, int _Dim, int _Degree>
typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivatives(
const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
const DenseIndex order,
const DenseIndex degree,
const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
{
typename SplineTraits<Spline>::BasisDerivativeType der;
BasisFunctionDerivativesImpl(u, order, degree, knots, der);
return der;
}
}

274
unsupported/Eigen/src/Splines/SplineFitting.h
View File

@@ -10,10 +10,14 @@
#ifndef EIGEN_SPLINE_FITTING_H
#define EIGEN_SPLINE_FITTING_H
#include <algorithm>
#include <functional>
#include <numeric>
#include <vector>
#include "SplineFwd.h"
#include <Eigen/LU>
#include <Eigen/QR>
namespace Eigen
@@ -49,6 +53,129 @@ namespace Eigen
knots.segment(knots.size()-degree-1,degree+1) = KnotVectorType::Ones(degree+1);
}
/**
* \brief Computes knot averages when derivative constraints are present.
* Note that this is a technical interpretation of the referenced article
* since the algorithm contained therein is incorrect as written.
* \ingroup Splines_Module
*
* \param[in] parameters The parameters at which the interpolation B-Spline
* will intersect the given interpolation points. The parameters
* are assumed to be a non-decreasing sequence.
* \param[in] degree The degree of the interpolating B-Spline. This must be
* greater than zero.
* \param[in] derivativeIndices The indices corresponding to parameters at
* which there are derivative constraints. The indices are assumed
* to be a non-decreasing sequence.
* \param[out] knots The calculated knot vector. These will be returned as a
* non-decreasing sequence
*
* \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
* Curve interpolation with directional constraints for engineering design.
* Engineering with Computers
**/
template <typename KnotVectorType, typename ParameterVectorType, typename IndexArray>
void KnotAveragingWithDerivatives(const ParameterVectorType& parameters,
const unsigned int degree,
const IndexArray& derivativeIndices,
KnotVectorType& knots)
{
typedef typename ParameterVectorType::Scalar Scalar;
DenseIndex numParameters = parameters.size();
DenseIndex numDerivatives = derivativeIndices.size();
if (numDerivatives < 1)
{
KnotAveraging(parameters, degree, knots);
return;
}
DenseIndex startIndex;
DenseIndex endIndex;
DenseIndex numInternalDerivatives = numDerivatives;
if (derivativeIndices[0] == 0)
{
startIndex = 0;
--numInternalDerivatives;
}
else
{
startIndex = 1;
}
if (derivativeIndices[numDerivatives - 1] == numParameters - 1)
{
endIndex = numParameters - degree;
--numInternalDerivatives;
}
else
{
endIndex = numParameters - degree - 1;
}
// There are (endIndex - startIndex + 1) knots obtained from the averaging
// and 2 for the first and last parameters.
DenseIndex numAverageKnots = endIndex - startIndex + 3;
KnotVectorType averageKnots(numAverageKnots);
averageKnots[0] = parameters[0];
int newKnotIndex = 0;
for (DenseIndex i = startIndex; i <= endIndex; ++i)
averageKnots[++newKnotIndex] = parameters.segment(i, degree).mean();
averageKnots[++newKnotIndex] = parameters[numParameters - 1];
newKnotIndex = -1;
ParameterVectorType temporaryParameters(numParameters + 1);
KnotVectorType derivativeKnots(numInternalDerivatives);
for (unsigned int i = 0; i < numAverageKnots - 1; ++i)
{
temporaryParameters[0] = averageKnots[i];
ParameterVectorType parameterIndices(numParameters);
int temporaryParameterIndex = 1;
for (int j = 0; j < numParameters; ++j)
{
Scalar parameter = parameters[j];
if (parameter >= averageKnots[i] && parameter < averageKnots[i + 1])
{
parameterIndices[temporaryParameterIndex] = j;
temporaryParameters[temporaryParameterIndex++] = parameter;
}
}
temporaryParameters[temporaryParameterIndex] = averageKnots[i + 1];
for (int j = 0; j <= temporaryParameterIndex - 2; ++j)
{
for (DenseIndex k = 0; k < derivativeIndices.size(); ++k)
{
if (parameterIndices[j + 1] == derivativeIndices[k]
&& parameterIndices[j + 1] != 0
&& parameterIndices[j + 1] != numParameters - 1)
{
derivativeKnots[++newKnotIndex] = temporaryParameters.segment(j, 3).mean();
break;
}
}
}
}
KnotVectorType temporaryKnots(averageKnots.size() + derivativeKnots.size());
std::merge(averageKnots.data(), averageKnots.data() + averageKnots.size(),
derivativeKnots.data(), derivativeKnots.data() + derivativeKnots.size(),
temporaryKnots.data());
// Number of control points (one for each point and derivative) plus spline order.
DenseIndex numKnots = numParameters + numDerivatives + degree + 1;
knots.resize(numKnots);
knots.head(degree).fill(temporaryKnots[0]);
knots.tail(degree).fill(temporaryKnots.template tail<1>()[0]);
knots.segment(degree, temporaryKnots.size()) = temporaryKnots;
}
/**
* \brief Computes chord length parameters which are required for spline interpolation.
* \ingroup Splines_Module
@@ -86,6 +213,7 @@ namespace Eigen
struct SplineFitting
{
typedef typename SplineType::KnotVectorType KnotVectorType;
typedef typename SplineType::ParameterVectorType ParameterVectorType;
/**
* \brief Fits an interpolating Spline to the given data points.
@@ -109,6 +237,52 @@ namespace Eigen
**/
template <typename PointArrayType>
static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters);
/**
* \brief Fits an interpolating spline to the given data points and
* derivatives.
*
* \param points The points for which an interpolating spline will be computed.
* \param derivatives The desired derivatives of the interpolating spline at interpolation
* points.
* \param derivativeIndices An array indicating which point each derivative belongs to. This
* must be the same size as @a derivatives.
* \param degree The degree of the interpolating spline.
*
* \returns A spline interpolating @a points with @a derivatives at those points.
*
* \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
* Curve interpolation with directional constraints for engineering design.
* Engineering with Computers
**/
template <typename PointArrayType, typename IndexArray>
static SplineType InterpolateWithDerivatives(const PointArrayType& points,
const PointArrayType& derivatives,
const IndexArray& derivativeIndices,
const unsigned int degree);
/**
* \brief Fits an interpolating spline to the given data points and derivatives.
*
* \param points The points for which an interpolating spline will be computed.
* \param derivatives The desired derivatives of the interpolating spline at interpolation points.
* \param derivativeIndices An array indicating which point each derivative belongs to. This
* must be the same size as @a derivatives.
* \param degree The degree of the interpolating spline.
* \param parameters The parameters corresponding to the interpolation points.
*
* \returns A spline interpolating @a points with @a derivatives at those points.
*
* \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
* Curve interpolation with directional constraints for engineering design.
* Engineering with Computers
*/
template <typename PointArrayType, typename IndexArray>
static SplineType InterpolateWithDerivatives(const PointArrayType& points,
const PointArrayType& derivatives,
const IndexArray& derivativeIndices,
const unsigned int degree,
const ParameterVectorType& parameters);
};
template <typename SplineType>
@@ -151,6 +325,106 @@ namespace Eigen
ChordLengths(pts, chord_lengths);
return Interpolate(pts, degree, chord_lengths);
}
template <typename SplineType>
template <typename PointArrayType, typename IndexArray>
SplineType
SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points,
const PointArrayType& derivatives,
const IndexArray& derivativeIndices,
const unsigned int degree,
const ParameterVectorType& parameters)
{
typedef typename SplineType::KnotVectorType::Scalar Scalar;
typedef typename SplineType::ControlPointVectorType ControlPointVectorType;
typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;
const DenseIndex n = points.cols() + derivatives.cols();
KnotVectorType knots;
KnotAveragingWithDerivatives(parameters, degree, derivativeIndices, knots);
// fill matrix
MatrixType A = MatrixType::Zero(n, n);
// Use these dimensions for quicker populating, then transpose for solving.
MatrixType b(points.rows(), n);
DenseIndex startRow;
DenseIndex derivativeStart;
// End derivatives.
if (derivativeIndices[0] == 0)
{
A.template block<1, 2>(1, 0) << -1, 1;
Scalar y = (knots(degree + 1) - knots(0)) / degree;
b.col(1) = y*derivatives.col(0);
startRow = 2;
derivativeStart = 1;
}
else
{
startRow = 1;
derivativeStart = 0;
}
if (derivativeIndices[derivatives.cols() - 1] == points.cols() - 1)
{
A.template block<1, 2>(n - 2, n - 2) << -1, 1;
Scalar y = (knots(knots.size() - 1) - knots(knots.size() - (degree + 2))) / degree;
b.col(b.cols() - 2) = y*derivatives.col(derivatives.cols() - 1);
}
DenseIndex row = startRow;
DenseIndex derivativeIndex = derivativeStart;
for (DenseIndex i = 1; i < parameters.size() - 1; ++i)
{
const DenseIndex span = SplineType::Span(parameters[i], degree, knots);
if (derivativeIndices[derivativeIndex] == i)
{
A.block(row, span - degree, 2, degree + 1)
= SplineType::BasisFunctionDerivatives(parameters[i], 1, degree, knots);
b.col(row++) = points.col(i);
b.col(row++) = derivatives.col(derivativeIndex++);
}
else
{
A.row(row++).segment(span - degree, degree + 1)
= SplineType::BasisFunctions(parameters[i], degree, knots);
}
}
b.col(0) = points.col(0);
b.col(b.cols() - 1) = points.col(points.cols() - 1);
A(0,0) = 1;
A(n - 1, n - 1) = 1;
// Solve
FullPivLU<MatrixType> lu(A);
ControlPointVectorType controlPoints = lu.solve(MatrixType(b.transpose())).transpose();
SplineType spline(knots, controlPoints);
return spline;
}
template <typename SplineType>
template <typename PointArrayType, typename IndexArray>
SplineType
SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points,
const PointArrayType& derivatives,
const IndexArray& derivativeIndices,
const unsigned int degree)
{
ParameterVectorType parameters;
ChordLengths(points, parameters);
return InterpolateWithDerivatives(points, derivatives, derivativeIndices, degree, parameters);
}
}
#endif // EIGEN_SPLINE_FITTING_H

3
unsupported/Eigen/src/Splines/SplineFwd.h
View File

@@ -48,6 +48,9 @@ namespace Eigen
/** \brief The data type used to store knot vectors. */
typedef Array<Scalar,1,Dynamic> KnotVectorType;
/** \brief The data type used to store parameter vectors. */
typedef Array<Scalar,1,Dynamic> ParameterVectorType;
/** \brief The data type representing the spline's control points. */
typedef Array<Scalar,Dimension,Dynamic> ControlPointVectorType;

11
unsupported/doc/Overview.dox
View File

@@ -1,14 +1,15 @@
/// \brief Namespace containing all symbols from the %Eigen library.
namespace Eigen {
/** \mainpage Eigen's unsupported modules
/** \mainpage %Eigen's unsupported modules
This is the API documentation for Eigen's unsupported modules.
This is the API documentation for %Eigen's unsupported modules.
These modules are contributions from various users. They are provided "as is", without any support.
Click on the \e Modules tab at the top of this page to get a list of all unsupported modules.
Don't miss the <a href="..//index.html">official Eigen documentation</a>.
Don't miss the <a href="../index.html">official Eigen documentation</a>.
*/
@@ -18,8 +19,10 @@ Don't miss the <a href="..//index.html">official Eigen documentation</a>.
The unsupported modules are contributions from various users. They are
provided "as is", without any support. Nevertheless, some of them are
subject to be included in Eigen in the future.
subject to be included in %Eigen in the future.
*/
/// \internal \brief Namespace containing low-level routines from the %Eigen library.
namespace internal {}
}

4
unsupported/test/CMakeLists.txt
View File

@@ -100,7 +100,7 @@ if(EIGEN_TEST_CXX11)
ei_add_test(cxx11_meta "-std=c++0x")
ei_add_test(cxx11_tensor_simple "-std=c++0x")
ei_add_test(cxx11_tensor_symmetry "-std=c++0x")
ei_add_test(cxx11_tensor_assign "-std=c++0x")
# ei_add_test(cxx11_tensor_assign "-std=c++0x")
ei_add_test(cxx11_tensor_comparisons "-std=c++0x")
ei_add_test(cxx11_tensor_contraction "-std=c++0x")
ei_add_test(cxx11_tensor_convolution "-std=c++0x")
@@ -108,7 +108,7 @@ if(EIGEN_TEST_CXX11)
# ei_add_test(cxx11_tensor_fixed_size "-std=c++0x")
ei_add_test(cxx11_tensor_lvalue "-std=c++0x")
ei_add_test(cxx11_tensor_map "-std=c++0x")
ei_add_test(cxx11_tensor_morphing "-std=c++0x")
# ei_add_test(cxx11_tensor_morphing "-std=c++0x")
# ei_add_test(cxx11_tensor_device "-std=c++0x")
# ei_add_test(cxx11_tensor_fixed_size "-std=c++0x")
# ei_add_test(cxx11_tensor_thread_pool "-std=c++0x")

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

File diff suppressed because it is too large Load Diff

3
unsupported/test/mpreal_support.cpp
View File

@@ -32,12 +32,11 @@ void test_mpreal_support()
VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs());
VERIFY_IS_APPROX(A.array().sin(), sin(A.array()));
VERIFY_IS_APPROX(A.array().cos(), cos(A.array()));
// Cholesky
X = S.selfadjointView<Lower>().llt().solve(B);
VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
// partial LU
X = A.lu().solve(B);
VERIFY_IS_APPROX((A*X).eval(),B);

73
unsupported/test/splines.cpp
View File

@@ -13,23 +13,23 @@
namespace Eigen {
// lets do some explicit instantiations and thus
// force the compilation of all spline functions...
template class Spline<double, 2, Dynamic>;
template class Spline<double, 3, Dynamic>;
// lets do some explicit instantiations and thus
// force the compilation of all spline functions...
template class Spline<double, 2, Dynamic>;
template class Spline<double, 3, Dynamic>;
template class Spline<double, 2, 2>;
template class Spline<double, 2, 3>;
template class Spline<double, 2, 4>;
template class Spline<double, 2, 5>;
template class Spline<double, 2, 2>;
template class Spline<double, 2, 3>;
template class Spline<double, 2, 4>;
template class Spline<double, 2, 5>;
template class Spline<float, 2, Dynamic>;
template class Spline<float, 3, Dynamic>;
template class Spline<float, 2, Dynamic>;
template class Spline<float, 3, Dynamic>;
template class Spline<float, 3, 2>;
template class Spline<float, 3, 3>;
template class Spline<float, 3, 4>;
template class Spline<float, 3, 5>;
template class Spline<float, 3, 2>;
template class Spline<float, 3, 3>;
template class Spline<float, 3, 4>;
template class Spline<float, 3, 5>;
}
@@ -234,11 +234,48 @@ void check_global_interpolation2d()
}
}
void check_global_interpolation_with_derivatives2d()
{
typedef Spline2d::PointType PointType;
typedef Spline2d::KnotVectorType KnotVectorType;
const unsigned int numPoints = 100;
const unsigned int dimension = 2;
const unsigned int degree = 3;
ArrayXXd points = ArrayXXd::Random(dimension, numPoints);
KnotVectorType knots;
Eigen::ChordLengths(points, knots);
ArrayXXd derivatives = ArrayXXd::Random(dimension, numPoints);
VectorXd derivativeIndices(numPoints);
for (Eigen::DenseIndex i = 0; i < numPoints; ++i)
derivativeIndices(i) = static_cast<double>(i);
const Spline2d spline = SplineFitting<Spline2d>::InterpolateWithDerivatives(
points, derivatives, derivativeIndices, degree);
for (Eigen::DenseIndex i = 0; i < points.cols(); ++i)
{
PointType point = spline(knots(i));
PointType referencePoint = points.col(i);
VERIFY_IS_APPROX(point, referencePoint);
PointType derivative = spline.derivatives(knots(i), 1).col(1);
PointType referenceDerivative = derivatives.col(i);
VERIFY_IS_APPROX(derivative, referenceDerivative);
}
}
void test_splines()
{
CALL_SUBTEST( eval_spline3d() );
CALL_SUBTEST( eval_spline3d_onbrks() );
CALL_SUBTEST( eval_closed_spline2d() );
CALL_SUBTEST( check_global_interpolation2d() );
for (int i = 0; i < g_repeat; ++i)
{
CALL_SUBTEST( eval_spline3d() );
CALL_SUBTEST( eval_spline3d_onbrks() );
CALL_SUBTEST( eval_closed_spline2d() );
CALL_SUBTEST( check_global_interpolation2d() );
CALL_SUBTEST( check_global_interpolation_with_derivatives2d() );
}
}