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

syncing this fork with upstream

Browse Source
This commit is contained in:
Deven Desai
2018-06-13 12:09:52 -04:00
parent 8fbd47052b d3a380af4d
commit d1d22ef0f4
88 changed files with 1266 additions and 395 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
unsupported/Eigen/src/BVH/KdBVH.h
View File

@@ -170,7 +170,7 @@ private:
typedef internal::vector_int_pair<Scalar, Dim> VIPair;
typedef std::vector<VIPair, aligned_allocator<VIPair> > VIPairList;
typedef Matrix<Scalar, Dim, 1> VectorType;
struct VectorComparator //compares vectors, or, more specificall, VIPairs along a particular dimension
struct VectorComparator //compares vectors, or more specifically, VIPairs along a particular dimension
{
VectorComparator(int inDim) : dim(inDim) {}
inline bool operator()(const VIPair &v1, const VIPair &v2) const { return v1.first[dim] < v2.first[dim]; }

2
unsupported/Eigen/src/Eigenvalues/ArpackSelfAdjointEigenSolver.h
View File

@@ -300,7 +300,7 @@ public:
/** \brief Reports whether previous computation was successful.
*
* \returns \c Success if computation was succesful, \c NoConvergence otherwise.
* \returns \c Success if computation was successful, \c NoConvergence otherwise.
*/
ComputationInfo info() const
{

2
unsupported/Eigen/src/EulerAngles/EulerSystem.h
View File

@@ -12,7 +12,7 @@
namespace Eigen
{
// Forward declerations
// Forward declarations
template <typename _Scalar, class _System>
class EulerAngles;

2
unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h
View File

@@ -99,7 +99,7 @@ void pseudo_inverse(const CMatrix &C, CINVMatrix &CINV)
/** \ingroup IterativeSolvers_Module
* Constrained conjugate gradient
*
* Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the contraint \f$ Cx \le f \f$
* Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the constraint \f$ Cx \le f \f$
*/
template<typename TMatrix, typename CMatrix,
typename VectorX, typename VectorB, typename VectorF>

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

@@ -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; }
void 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,27 +200,27 @@ 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
void dgmresInitDeflation(Index& rows) const;
mutable DenseMatrix m_V; // Krylov basis vectors
mutable DenseMatrix m_H; // Hessenberg matrix
mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied
mutable DenseMatrix m_Hes; // Initial hessenberg matrix without Givens rotations applied
mutable Index m_restart; // Maximum size of the Krylov subspace
mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace
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 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,13 +242,13 @@ 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);
//Initial residual vector and intial norm
//Initial residual vector and initial norm
x = precond.solve(x);
r0 = rhs - mat * x;
RealScalar beta = r0.norm();
@@ -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);

2
unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h
View File

@@ -73,7 +73,7 @@ void lmqrsolv(
qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj;
wa[k] = temp;
/* accumulate the tranformation in the row of s. */
/* accumulate the transformation in the row of s. */
for (i = k+1; i<n; ++i) {
temp = givens.c() * s(i,k) + givens.s() * sdiag[i];
sdiag[i] = -givens.s() * s(i,k) + givens.c() * sdiag[i];

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

@@ -233,9 +233,9 @@ class LevenbergMarquardt : internal::no_assignment_operator
/**
* \brief Reports whether the minimization was successful
* \returns \c Success if the minimization was succesful,
* \returns \c Success if the minimization was successful,
* \c NumericalIssue if a numerical problem arises during the
* minimization process, for exemple during the QR factorization
* minimization process, for example during the QR factorization
* \c NoConvergence if the minimization did not converge after
* the maximum number of function evaluation allowed
* \c InvalidInput if the input matrix is invalid

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

@@ -313,7 +313,7 @@ struct matrix_exp_computeUV<MatrixType, long double>
matrix_exp_pade17(A, U, V);
}
#elif LDBL_MANT_DIG <= 112 // quadruple precison
#elif LDBL_MANT_DIG <= 112 // quadruple precision
if (l1norm < 1.639394610288918690547467954466970e-005L) {
matrix_exp_pade3(arg, U, V);

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

@@ -81,7 +81,7 @@ class MatrixPowerParenthesesReturnValue : public ReturnByValue< MatrixPowerParen
*
* \note Currently this class is only used by MatrixPower. One may
* insist that this be nested into MatrixPower. This class is here to
* faciliate future development of triangular matrix functions.
* facilitate future development of triangular matrix functions.
*/
template<typename MatrixType>
class MatrixPowerAtomic : internal::noncopyable

2
unsupported/Eigen/src/NonLinearOptimization/qrsolv.h
View File

@@ -61,7 +61,7 @@ void qrsolv(
qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj;
wa[k] = temp;
/* accumulate the tranformation in the row of s. */
/* accumulate the transformation in the row of s. */
for (i = k+1; i<n; ++i) {
temp = givens.c() * s(i,k) + givens.s() * sdiag[i];
sdiag[i] = -givens.s() * s(i,k) + givens.c() * sdiag[i];

2
unsupported/Eigen/src/NonLinearOptimization/r1updt.h
View File

@@ -22,7 +22,7 @@ void r1updt(
Scalar temp;
JacobiRotation<Scalar> givens;
// r1updt had a broader usecase, but we dont use it here. And, more
// r1updt had a broader usecase, but we don't use it here. And, more
// importantly, we can not test it.
eigen_assert(m==n);
eigen_assert(u.size()==m);

4
unsupported/Eigen/src/Polynomials/Companion.h
View File

@@ -104,7 +104,7 @@ class companion
/** Helper function for the balancing algorithm.
* \returns true if the row and the column, having colNorm and rowNorm
* as norms, are balanced, false otherwise.
* colB and rowB are repectively the multipliers for
* colB and rowB are respectively the multipliers for
* the column and the row in order to balance them.
* */
bool balanced( RealScalar colNorm, RealScalar rowNorm,
@@ -113,7 +113,7 @@ class companion
/** Helper function for the balancing algorithm.
* \returns true if the row and the column, having colNorm and rowNorm
* as norms, are balanced, false otherwise.
* colB and rowB are repectively the multipliers for
* colB and rowB are respectively the multipliers for
* the column and the row in order to balance them.
* */
bool balancedR( RealScalar colNorm, RealScalar rowNorm,

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

@@ -41,7 +41,7 @@ public:
/** Sets the relative threshold value used to prune zero coefficients during the decomposition.
*
* Setting a value greater than zero speeds up computation, and yields to an imcomplete
* Setting a value greater than zero speeds up computation, and yields to an incomplete
* factorization with fewer non zero coefficients. Such approximate factors are especially
* useful to initialize an iterative solver.
*

16
unsupported/Eigen/src/Skyline/SkylineMatrix.h
View File

@@ -206,26 +206,26 @@ public:
if (col > row) //upper matrix
{
const Index minOuterIndex = inner - m_data.upperProfile(inner);
eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(outer >= minOuterIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
}
if (col < row) //lower matrix
{
const Index minInnerIndex = outer - m_data.lowerProfile(outer);
eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(inner >= minInnerIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
}
} else {
if (outer > inner) //upper matrix
{
const Index maxOuterIndex = inner + m_data.upperProfile(inner);
eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(outer <= maxOuterIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
}
if (outer < inner) //lower matrix
{
const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(inner <= maxInnerIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
}
}
@@ -300,11 +300,11 @@ public:
if (IsRowMajor) {
const Index minInnerIndex = outer - m_data.lowerProfile(outer);
eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(inner >= minInnerIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
} else {
const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(inner <= maxInnerIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
}
}
@@ -336,11 +336,11 @@ public:
if (IsRowMajor) {
const Index minOuterIndex = inner - m_data.upperProfile(inner);
eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(outer >= minOuterIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
} else {
const Index maxOuterIndex = inner + m_data.upperProfile(inner);
eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage");
eigen_assert(outer <= maxOuterIndex && "You tried to access a coeff that does not exist in the storage");
return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
}
}

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

@@ -187,7 +187,7 @@ template<typename _Scalar, int _Options, typename _StorageIndex>
/** Does nothing: provided for compatibility with SparseMatrix */
inline void finalize() {}
/** Suppress all nonzeros which are smaller than \a reference under the tolerence \a epsilon */
/** Suppress all nonzeros which are smaller than \a reference under the tolerance \a epsilon */
void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
{
for (Index j=0; j<outerSize(); ++j)
@@ -224,21 +224,21 @@ template<typename _Scalar, int _Options, typename _StorageIndex>
}
}
/** The class DynamicSparseMatrix is deprectaed */
/** The class DynamicSparseMatrix is deprecated */
EIGEN_DEPRECATED inline DynamicSparseMatrix()
: m_innerSize(0), m_data(0)
{
eigen_assert(innerSize()==0 && outerSize()==0);
}
/** The class DynamicSparseMatrix is deprectaed */
/** The class DynamicSparseMatrix is deprecated */
EIGEN_DEPRECATED inline DynamicSparseMatrix(Index rows, Index cols)
: m_innerSize(0)
{
resize(rows, cols);
}
/** The class DynamicSparseMatrix is deprectaed */
/** The class DynamicSparseMatrix is deprecated */
template<typename OtherDerived>
EIGEN_DEPRECATED explicit inline DynamicSparseMatrix(const SparseMatrixBase<OtherDerived>& other)
: m_innerSize(0)

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

@@ -104,7 +104,7 @@ namespace internal
out << value.real << " " << value.imag()<< "\n";
}
} // end namepsace internal
} // end namespace internal
inline bool getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
{

42
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h
View File

@@ -33,6 +33,48 @@ igamma(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerive
);
}
/** \cpp11 \returns an expression of the coefficient-wise igamma_der_a(\a a, \a x) to the given arrays.
*
* This function computes the coefficient-wise derivative of the incomplete
* gamma function with respect to the parameter a.
*
* \note This function supports only float and double scalar types in c++11
* mode. To support other scalar types,
* or float/double in non c++11 mode, the user has to provide implementations
* of igamma_der_a(T,T) for any scalar
* type T to be supported.
*
* \sa Eigen::igamma(), Eigen::lgamma()
*/
template <typename Derived, typename ExponentDerived>
inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_der_a_op<typename Derived::Scalar>, const Derived, const ExponentDerived>
igamma_der_a(const Eigen::ArrayBase<Derived>& a, const Eigen::ArrayBase<ExponentDerived>& x) {
return Eigen::CwiseBinaryOp<Eigen::internal::scalar_igamma_der_a_op<typename Derived::Scalar>, const Derived, const ExponentDerived>(
a.derived(),
x.derived());
}
/** \cpp11 \returns an expression of the coefficient-wise gamma_sample_der_alpha(\a alpha, \a sample) to the given arrays.
*
* This function computes the coefficient-wise derivative of the sample
* of a Gamma(alpha, 1) random variable with respect to the parameter alpha.
*
* \note This function supports only float and double scalar types in c++11
* mode. To support other scalar types,
* or float/double in non c++11 mode, the user has to provide implementations
* of gamma_sample_der_alpha(T,T) for any scalar
* type T to be supported.
*
* \sa Eigen::igamma(), Eigen::lgamma()
*/
template <typename AlphaDerived, typename SampleDerived>
inline const Eigen::CwiseBinaryOp<Eigen::internal::scalar_gamma_sample_der_alpha_op<typename AlphaDerived::Scalar>, const AlphaDerived, const SampleDerived>
gamma_sample_der_alpha(const Eigen::ArrayBase<AlphaDerived>& alpha, const Eigen::ArrayBase<SampleDerived>& sample) {
return Eigen::CwiseBinaryOp<Eigen::internal::scalar_gamma_sample_der_alpha_op<typename AlphaDerived::Scalar>, const AlphaDerived, const SampleDerived>(
alpha.derived(),
sample.derived());
}
/** \cpp11 \returns an expression of the coefficient-wise igammac(\a a, \a x) to the given arrays.
*
* This function computes the coefficient-wise complementary incomplete gamma function.

54
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h
View File

@@ -41,6 +41,60 @@ struct functor_traits<scalar_igamma_op<Scalar> > {
};
};
/** \internal
* \brief Template functor to compute the derivative of the incomplete gamma
* function igamma_der_a(a, x)
*
* \sa class CwiseBinaryOp, Cwise::igamma_der_a
*/
template <typename Scalar>
struct scalar_igamma_der_a_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_der_a_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& a, const Scalar& x) const {
using numext::igamma_der_a;
return igamma_der_a(a, x);
}
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const {
return internal::pigamma_der_a(a, x);
}
};
template <typename Scalar>
struct functor_traits<scalar_igamma_der_a_op<Scalar> > {
enum {
// 2x the cost of igamma
Cost = 40 * NumTraits<Scalar>::MulCost + 20 * NumTraits<Scalar>::AddCost,
PacketAccess = packet_traits<Scalar>::HasIGammaDerA
};
};
/** \internal
* \brief Template functor to compute the derivative of the sample
* of a Gamma(alpha, 1) random variable with respect to the parameter alpha
* gamma_sample_der_alpha(alpha, sample)
*
* \sa class CwiseBinaryOp, Cwise::gamma_sample_der_alpha
*/
template <typename Scalar>
struct scalar_gamma_sample_der_alpha_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_gamma_sample_der_alpha_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& alpha, const Scalar& sample) const {
using numext::gamma_sample_der_alpha;
return gamma_sample_der_alpha(alpha, sample);
}
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& alpha, const Packet& sample) const {
return internal::pgamma_sample_der_alpha(alpha, sample);
}
};
template <typename Scalar>
struct functor_traits<scalar_gamma_sample_der_alpha_op<Scalar> > {
enum {
// 2x the cost of igamma, minus the lgamma cost (the lgamma cancels out)
Cost = 30 * NumTraits<Scalar>::MulCost + 15 * NumTraits<Scalar>::AddCost,
PacketAccess = packet_traits<Scalar>::HasGammaSampleDerAlpha
};
};
/** \internal
* \brief Template functor to compute the complementary incomplete gamma function igammac(a, x)

8
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h
View File

@@ -33,6 +33,14 @@ template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half erfc(const Eigen::h
template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igamma(const Eigen::half& a, const Eigen::half& x) {
return Eigen::half(Eigen::numext::igamma(static_cast<float>(a), static_cast<float>(x)));
}
template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igamma_der_a(const Eigen::half& a, const Eigen::half& x) {
return Eigen::half(Eigen::numext::igamma_der_a(static_cast<float>(a), static_cast<float>(x)));
}
template <>
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half gamma_sample_der_alpha(const Eigen::half& alpha, const Eigen::half& sample) {
return Eigen::half(Eigen::numext::gamma_sample_der_alpha(static_cast<float>(alpha), static_cast<float>(sample)));
}
template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igammac(const Eigen::half& a, const Eigen::half& x) {
return Eigen::half(Eigen::numext::igammac(static_cast<float>(a), static_cast<float>(x)));
}

576
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
View File

@@ -521,6 +521,197 @@ struct cephes_helper<double> {
}
};
enum IgammaComputationMode { VALUE, DERIVATIVE, SAMPLE_DERIVATIVE };
template <typename Scalar, IgammaComputationMode mode>
EIGEN_DEVICE_FUNC
int igamma_num_iterations() {
/* Returns the maximum number of internal iterations for igamma computation.
*/
if (mode == VALUE) {
return 2000;
}
if (internal::is_same<Scalar, float>::value) {
return 200;
} else if (internal::is_same<Scalar, double>::value) {
return 500;
} else {
return 2000;
}
}
template <typename Scalar, IgammaComputationMode mode>
struct igammac_cf_impl {
/* Computes igamc(a, x) or derivative (depending on the mode)
* using the continued fraction expansion of the complementary
* incomplete Gamma function.
*
* Preconditions:
* a > 0
* x >= 1
* x >= a
*/
EIGEN_DEVICE_FUNC
static Scalar run(Scalar a, Scalar x) {
const Scalar zero = 0;
const Scalar one = 1;
const Scalar two = 2;
const Scalar machep = cephes_helper<Scalar>::machep();
const Scalar big = cephes_helper<Scalar>::big();
const Scalar biginv = cephes_helper<Scalar>::biginv();
if ((numext::isinf)(x)) {
return zero;
}
// continued fraction
Scalar y = one - a;
Scalar z = x + y + one;
Scalar c = zero;
Scalar pkm2 = one;
Scalar qkm2 = x;
Scalar pkm1 = x + one;
Scalar qkm1 = z * x;
Scalar ans = pkm1 / qkm1;
Scalar dpkm2_da = zero;
Scalar dqkm2_da = zero;
Scalar dpkm1_da = zero;
Scalar dqkm1_da = -x;
Scalar dans_da = (dpkm1_da - ans * dqkm1_da) / qkm1;
for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) {
c += one;
y += one;
z += two;
Scalar yc = y * c;
Scalar pk = pkm1 * z - pkm2 * yc;
Scalar qk = qkm1 * z - qkm2 * yc;
Scalar dpk_da = dpkm1_da * z - pkm1 - dpkm2_da * yc + pkm2 * c;
Scalar dqk_da = dqkm1_da * z - qkm1 - dqkm2_da * yc + qkm2 * c;
if (qk != zero) {
Scalar ans_prev = ans;
ans = pk / qk;
Scalar dans_da_prev = dans_da;
dans_da = (dpk_da - ans * dqk_da) / qk;
if (mode == VALUE) {
if (numext::abs(ans_prev - ans) <= machep * numext::abs(ans)) {
break;
}
} else {
if (numext::abs(dans_da - dans_da_prev) <= machep) {
break;
}
}
}
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
dpkm2_da = dpkm1_da;
dpkm1_da = dpk_da;
dqkm2_da = dqkm1_da;
dqkm1_da = dqk_da;
if (numext::abs(pk) > big) {
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
dpkm2_da *= biginv;
dpkm1_da *= biginv;
dqkm2_da *= biginv;
dqkm1_da *= biginv;
}
}
/* Compute x**a * exp(-x) / gamma(a) */
Scalar logax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a);
Scalar ax = numext::exp(logax);
Scalar dax_da = ax * dlogax_da;
switch (mode) {
case VALUE:
return ans * ax;
case DERIVATIVE:
return ans * dax_da + dans_da * ax;
case SAMPLE_DERIVATIVE:
return -(dans_da + ans * dlogax_da) * x;
}
}
};
template <typename Scalar, IgammaComputationMode mode>
struct igamma_series_impl {
/* Computes igam(a, x) or its derivative (depending on the mode)
* using the series expansion of the incomplete Gamma function.
*
* Preconditions:
* x > 0
* a > 0
* !(x > 1 && x > a)
*/
EIGEN_DEVICE_FUNC
static Scalar run(Scalar a, Scalar x) {
const Scalar zero = 0;
const Scalar one = 1;
const Scalar machep = cephes_helper<Scalar>::machep();
/* power series */
Scalar r = a;
Scalar c = one;
Scalar ans = one;
Scalar dc_da = zero;
Scalar dans_da = zero;
for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) {
r += one;
Scalar term = x / r;
Scalar dterm_da = -x / (r * r);
dc_da = term * dc_da + dterm_da * c;
dans_da += dc_da;
c *= term;
ans += c;
if (mode == VALUE) {
if (c <= machep * ans) {
break;
}
} else {
if (numext::abs(dc_da) <= machep * numext::abs(dans_da)) {
break;
}
}
}
/* Compute x**a * exp(-x) / gamma(a + 1) */
Scalar logax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a + one);
Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a + one);
Scalar ax = numext::exp(logax);
Scalar dax_da = ax * dlogax_da;
switch (mode) {
case VALUE:
return ans * ax;
case DERIVATIVE:
return ans * dax_da + dans_da * ax;
case SAMPLE_DERIVATIVE:
return -(dans_da + ans * dlogax_da) * x / a;
}
}
};
#if !EIGEN_HAS_C99_MATH
template <typename Scalar>
@@ -535,8 +726,6 @@ struct igammac_impl {
#else
template <typename Scalar> struct igamma_impl; // predeclare igamma_impl
template <typename Scalar>
struct igammac_impl {
EIGEN_DEVICE_FUNC
@@ -604,97 +793,15 @@ struct igammac_impl {
return nan;
}
if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans
if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans
return nan;
}
if ((x < one) || (x < a)) {
/* The checks above ensure that we meet the preconditions for
* igamma_impl::Impl(), so call it, rather than igamma_impl::Run().
* Calling Run() would also work, but in that case the compiler may not be
* able to prove that igammac_impl::Run and igamma_impl::Run are not
* mutually recursive. This leads to worse code, particularly on
* platforms like nvptx, where recursion is allowed only begrudgingly.
*/
return (one - igamma_impl<Scalar>::Impl(a, x));
return (one - igamma_series_impl<Scalar, VALUE>::run(a, x));
}
return Impl(a, x);
}
private:
/* igamma_impl calls igammac_impl::Impl. */
friend struct igamma_impl<Scalar>;
/* Actually computes igamc(a, x).
*
* Preconditions:
* a > 0
* x >= 1
* x >= a
*/
EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) {
const Scalar zero = 0;
const Scalar one = 1;
const Scalar two = 2;
const Scalar machep = cephes_helper<Scalar>::machep();
const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
const Scalar big = cephes_helper<Scalar>::big();
const Scalar biginv = cephes_helper<Scalar>::biginv();
const Scalar inf = NumTraits<Scalar>::infinity();
Scalar ans, ax, c, yc, r, t, y, z;
Scalar pk, pkm1, pkm2, qk, qkm1, qkm2;
if (x == inf) return zero; // std::isinf crashes on CUDA
/* Compute x**a * exp(-x) / gamma(a) */
ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
if (ax < -maxlog) { // underflow
return zero;
}
ax = numext::exp(ax);
// continued fraction
y = one - a;
z = x + y + one;
c = zero;
pkm2 = one;
qkm2 = x;
pkm1 = x + one;
qkm1 = z * x;
ans = pkm1 / qkm1;
for (int i = 0; i < 2000; i++) {
c += one;
y += one;
z += two;
yc = y * c;
pk = pkm1 * z - pkm2 * yc;
qk = qkm1 * z - qkm2 * yc;
if (qk != zero) {
r = pk / qk;
t = numext::abs((ans - r) / r);
ans = r;
} else {
t = one;
}
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if (numext::abs(pk) > big) {
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
if (t <= machep) {
break;
}
}
return (ans * ax);
return igammac_cf_impl<Scalar, VALUE>::run(a, x);
}
};
@@ -704,15 +811,10 @@ struct igammac_impl {
* Implementation of igamma (incomplete gamma integral), based on Cephes but requires C++11/C99 *
************************************************************************************************/
template <typename Scalar>
struct igamma_retval {
typedef Scalar type;
};
#if !EIGEN_HAS_C99_MATH
template <typename Scalar>
struct igamma_impl {
template <typename Scalar, IgammaComputationMode mode>
struct igamma_generic_impl {
EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) {
EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false),
@@ -723,69 +825,17 @@ struct igamma_impl {
#else
template <typename Scalar>
struct igamma_impl {
template <typename Scalar, IgammaComputationMode mode>
struct igamma_generic_impl {
EIGEN_DEVICE_FUNC
static Scalar run(Scalar a, Scalar x) {
/* igam()
* Incomplete gamma integral
*
*
*
* SYNOPSIS:
*
* double a, x, y, igam();
*
* y = igam( a, x );
*
* DESCRIPTION:
*
* The function is defined by
*
* x
* -
* 1 | | -t a-1
* igam(a,x) = ----- | e t dt.
* - | |
* | (a) -
* 0
*
*
* In this implementation both arguments must be positive.
* The integral is evaluated by either a power series or
* continued fraction expansion, depending on the relative
* values of a and x.
*
* ACCURACY (double):
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 200000 3.6e-14 2.9e-15
* IEEE 0,100 300000 9.9e-14 1.5e-14
*
*
* ACCURACY (float):
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 20000 7.8e-6 5.9e-7
*
*/
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1985, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* left tail of incomplete gamma function:
*
* inf. k
* a -x - x
* x e > ----------
* - -
* k=0 | (a+k+1)
/* Depending on the mode, returns
* - VALUE: incomplete Gamma function igamma(a, x)
* - DERIVATIVE: derivative of incomplete Gamma function d/da igamma(a, x)
* - SAMPLE_DERIVATIVE: implicit derivative of a Gamma random variable
* x ~ Gamma(x | a, 1), dx/da = -1 / Gamma(x | a, 1) * d igamma(a, x) / dx
*
* Derivatives are implemented by forward-mode differentiation.
*/
const Scalar zero = 0;
const Scalar one = 1;
@@ -797,71 +847,167 @@ struct igamma_impl {
return nan;
}
if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans
if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans
return nan;
}
if ((x > one) && (x > a)) {
/* The checks above ensure that we meet the preconditions for
* igammac_impl::Impl(), so call it, rather than igammac_impl::Run().
* Calling Run() would also work, but in that case the compiler may not be
* able to prove that igammac_impl::Run and igamma_impl::Run are not
* mutually recursive. This leads to worse code, particularly on
* platforms like nvptx, where recursion is allowed only begrudgingly.
*/
return (one - igammac_impl<Scalar>::Impl(a, x));
}
return Impl(a, x);
}
private:
/* igammac_impl calls igamma_impl::Impl. */
friend struct igammac_impl<Scalar>;
/* Actually computes igam(a, x).
*
* Preconditions:
* x > 0
* a > 0
* !(x > 1 && x > a)
*/
EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) {
const Scalar zero = 0;
const Scalar one = 1;
const Scalar machep = cephes_helper<Scalar>::machep();
const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
Scalar ans, ax, c, r;
/* Compute x**a * exp(-x) / gamma(a) */
ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
if (ax < -maxlog) {
// underflow
return zero;
}
ax = numext::exp(ax);
/* power series */
r = a;
c = one;
ans = one;
for (int i = 0; i < 2000; i++) {
r += one;
c *= x/r;
ans += c;
if (c/ans <= machep) {
break;
Scalar ret = igammac_cf_impl<Scalar, mode>::run(a, x);
if (mode == VALUE) {
return one - ret;
} else {
return -ret;
}
}
return (ans * ax / a);
return igamma_series_impl<Scalar, mode>::run(a, x);
}
};
#endif // EIGEN_HAS_C99_MATH
template <typename Scalar>
struct igamma_retval {
typedef Scalar type;
};
template <typename Scalar>
struct igamma_impl : igamma_generic_impl<Scalar, VALUE> {
/* igam()
* Incomplete gamma integral.
*
* The CDF of Gamma(a, 1) random variable at the point x.
*
* Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample
* 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points.
* The ground truth is computed by mpmath. Mean absolute error:
* float: 1.26713e-05
* double: 2.33606e-12
*
* Cephes documentation below.
*
* SYNOPSIS:
*
* double a, x, y, igam();
*
* y = igam( a, x );
*
* DESCRIPTION:
*
* The function is defined by
*
* x
* -
* 1 | | -t a-1
* igam(a,x) = ----- | e t dt.
* - | |
* | (a) -
* 0
*
*
* In this implementation both arguments must be positive.
* The integral is evaluated by either a power series or
* continued fraction expansion, depending on the relative
* values of a and x.
*
* ACCURACY (double):
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 200000 3.6e-14 2.9e-15
* IEEE 0,100 300000 9.9e-14 1.5e-14
*
*
* ACCURACY (float):
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 20000 7.8e-6 5.9e-7
*
*/
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1985, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* left tail of incomplete gamma function:
*
* inf. k
* a -x - x
* x e > ----------
* - -
* k=0 | (a+k+1)
*
*/
};
template <typename Scalar>
struct igamma_der_a_retval : igamma_retval<Scalar> {};
template <typename Scalar>
struct igamma_der_a_impl : igamma_generic_impl<Scalar, DERIVATIVE> {
/* Derivative of the incomplete Gamma function with respect to a.
*
* Computes d/da igamma(a, x) by forward differentiation of the igamma code.
*
* Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample
* 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points.
* The ground truth is computed by mpmath. Mean absolute error:
* float: 6.17992e-07
* double: 4.60453e-12
*
* Reference:
* R. Moore. "Algorithm AS 187: Derivatives of the incomplete gamma
* integral". Journal of the Royal Statistical Society. 1982
*/
};
template <typename Scalar>
struct gamma_sample_der_alpha_retval : igamma_retval<Scalar> {};
template <typename Scalar>
struct gamma_sample_der_alpha_impl
: igamma_generic_impl<Scalar, SAMPLE_DERIVATIVE> {
/* Derivative of a Gamma random variable sample with respect to alpha.
*
* Consider a sample of a Gamma random variable with the concentration
* parameter alpha: sample ~ Gamma(alpha, 1). The reparameterization
* derivative that we want to compute is dsample / dalpha =
* d igammainv(alpha, u) / dalpha, where u = igamma(alpha, sample).
* However, this formula is numerically unstable and expensive, so instead
* we use implicit differentiation:
*
* igamma(alpha, sample) = u, where u ~ Uniform(0, 1).
* Apply d / dalpha to both sides:
* d igamma(alpha, sample) / dalpha
* + d igamma(alpha, sample) / dsample * dsample/dalpha = 0
* d igamma(alpha, sample) / dalpha
* + Gamma(sample | alpha, 1) dsample / dalpha = 0
* dsample/dalpha = - (d igamma(alpha, sample) / dalpha)
* / Gamma(sample | alpha, 1)
*
* Here Gamma(sample | alpha, 1) is the PDF of the Gamma distribution
* (note that the derivative of the CDF w.r.t. sample is the PDF).
* See the reference below for more details.
*
* The derivative of igamma(alpha, sample) is computed by forward
* differentiation of the igamma code. Division by the Gamma PDF is performed
* in the same code, increasing the accuracy and speed due to cancellation
* of some terms.
*
* Accuracy estimation. For each alpha in [10^-2, 10^-1...10^3] we sample
* 50 Gamma random variables sample ~ Gamma(sample | alpha, 1), a total of 300
* points. The ground truth is computed by mpmath. Mean absolute error:
* float: 2.1686e-06
* double: 1.4774e-12
*
* Reference:
* M. Figurnov, S. Mohamed, A. Mnih "Implicit Reparameterization Gradients".
* 2018
*/
};
/*****************************************************************************
* Implementation of Riemann zeta function of two arguments, based on Cephes *
*****************************************************************************/
@@ -1574,6 +1720,8 @@ struct betainc_impl<double> {
}
};
#endif // EIGEN_HAS_C99_MATH
/****************************************************************************
* Implementation of Bessel function, based on Cephes *
****************************************************************************/
@@ -1902,8 +2050,6 @@ struct i1e_impl<double> {
}
};
#endif // EIGEN_HAS_C99_MATH
} // end namespace internal
namespace numext {
@@ -1950,6 +2096,18 @@ EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma, Scalar)
return EIGEN_MATHFUNC_IMPL(igamma, Scalar)::run(a, x);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma_der_a, Scalar)
igamma_der_a(const Scalar& a, const Scalar& x) {
return EIGEN_MATHFUNC_IMPL(igamma_der_a, Scalar)::run(a, x);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(gamma_sample_der_alpha, Scalar)
gamma_sample_der_alpha(const Scalar& a, const Scalar& x) {
return EIGEN_MATHFUNC_IMPL(gamma_sample_der_alpha, Scalar)::run(a, x);
}
template <typename Scalar>
EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igammac, Scalar)
igammac(const Scalar& a, const Scalar& x) {

15
unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h
View File

@@ -42,6 +42,21 @@ Packet perfc(const Packet& a) { using numext::erfc; return erfc(a); }
template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Packet pigamma(const Packet& a, const Packet& x) { using numext::igamma; return igamma(a, x); }
/** \internal \returns the derivative of the incomplete gamma function
* igamma_der_a(\a a, \a x) */
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pigamma_der_a(const Packet& a, const Packet& x) {
using numext::igamma_der_a; return igamma_der_a(a, x);
}
/** \internal \returns compute the derivative of the sample
* of Gamma(alpha, 1) random variable with respect to the parameter a
* gamma_sample_der_alpha(\a alpha, \a sample) */
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pgamma_sample_der_alpha(const Packet& alpha, const Packet& sample) {
using numext::gamma_sample_der_alpha; return gamma_sample_der_alpha(alpha, sample);
}
/** \internal \returns the complementary incomplete gamma function igammac(\a a, \a x) */
template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Packet pigammac(const Packet& a, const Packet& x) { using numext::igammac; return igammac(a, x); }

35
unsupported/Eigen/src/SpecialFunctions/arch/CUDA/CudaSpecialFunctions.h
View File

@@ -120,6 +120,41 @@ double2 pigamma<double2>(const double2& a, const double2& x)
return make_double2(igamma(a.x, x.x), igamma(a.y, x.y));
}
template <>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pigamma_der_a<float4>(
const float4& a, const float4& x) {
using numext::igamma_der_a;
return make_float4(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y),
igamma_der_a(a.z, x.z), igamma_der_a(a.w, x.w));
}
template <>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
pigamma_der_a<double2>(const double2& a, const double2& x) {
using numext::igamma_der_a;
return make_double2(igamma_der_a(a.x, x.x), igamma_der_a(a.y, x.y));
}
template <>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pgamma_sample_der_alpha<float4>(
const float4& alpha, const float4& sample) {
using numext::gamma_sample_der_alpha;
return make_float4(
gamma_sample_der_alpha(alpha.x, sample.x),
gamma_sample_der_alpha(alpha.y, sample.y),
gamma_sample_der_alpha(alpha.z, sample.z),
gamma_sample_der_alpha(alpha.w, sample.w));
}
template <>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
pgamma_sample_der_alpha<double2>(const double2& alpha, const double2& sample) {
using numext::gamma_sample_der_alpha;
return make_double2(
gamma_sample_der_alpha(alpha.x, sample.x),
gamma_sample_der_alpha(alpha.y, sample.y));
}
template<> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
float4 pigammac<float4>(const float4& a, const float4& x)
{

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

@@ -181,7 +181,7 @@ namespace Eigen
* \ingroup Splines_Module
*
* \param[in] pts The data points to which a spline should be fit.
* \param[out] chord_lengths The resulting chord lenggth vector.
* \param[out] chord_lengths The resulting chord length vector.
*
* \sa Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data
**/