GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 19:59:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

synch with main branch

Browse Source
This commit is contained in:
Gael Guennebaud 2009-07-28 17:35:07 +02:00
commit 54804eb626
19 changed files with 623 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
CTestConfig.cmake
View File

@ -5,9 +5,9 @@
## ENABLE_TESTING()
## INCLUDE(CTest)
set(CTEST_PROJECT_NAME "Eigen")
set(CTEST_NIGHTLY_START_TIME "05:00:00 UTC")
set(CTEST_NIGHTLY_START_TIME "06:00:00 UTC")
set(CTEST_DROP_METHOD "http")
set(CTEST_DROP_SITE "my.cdash.org")
set(CTEST_DROP_LOCATION "/submit.php?project=Eigen")
set(CTEST_DROP_SITE "eigen.tuxfamily.org")
set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen")
set(CTEST_DROP_SITE_CDASH TRUE)

1
Eigen/Core
View File

@ -161,6 +161,7 @@ namespace Eigen {
#include "src/Core/CwiseNullaryOp.h"
#include "src/Core/CwiseUnaryView.h"
#include "src/Core/Dot.h"
#include "src/Core/StableNorm.h"
#include "src/Core/MapBase.h"
#include "src/Core/Map.h"
#include "src/Core/Block.h"

2
Eigen/QR
View File

@ -41,7 +41,7 @@ namespace Eigen {
// declare all classes for a given matrix type
#define EIGEN_QR_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \
PREFIX template class QR<MATRIXTYPE>; \
PREFIX template class HouseholderQR<MATRIXTYPE>; \
PREFIX template class Tridiagonalization<MATRIXTYPE>; \
PREFIX template class HessenbergDecomposition<MATRIXTYPE>; \
PREFIX template class SelfAdjointEigenSolver<MATRIXTYPE>

22
Eigen/src/Cholesky/LDLT.h
View File

@ -62,7 +62,7 @@ template<typename MatrixType> class LDLT
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<int, 1, MatrixType::RowsAtCompileTime> IntRowVectorType;
/**
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
@ -83,7 +83,7 @@ template<typename MatrixType> class LDLT
inline TriangularView<MatrixType, UnitLowerTriangular> matrixL(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix;
return m_matrix;
}
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
@ -97,24 +97,24 @@ template<typename MatrixType> class LDLT
}
/** \returns the coefficients of the diagonal matrix D */
inline Diagonal<MatrixType,0> vectorD(void) const
{
inline Diagonal<MatrixType,0> vectorD(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix.diagonal();
return m_matrix.diagonal();
}
/** \returns true if the matrix is positive (semidefinite) */
inline bool isPositive(void) const
{
inline bool isPositive(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == 1;
return m_sign == 1;
}
/** \returns true if the matrix is negative (semidefinite) */
inline bool isNegative(void) const
{
inline bool isNegative(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == -1;
return m_sign == -1;
}
template<typename RhsDerived, typename ResultType>

4
Eigen/src/Core/Block.h
View File

@ -33,8 +33,8 @@
* \param MatrixType the type of the object in which we are taking a block
* \param BlockRows the number of rows of the block we are taking at compile time (optional)
* \param BlockCols the number of columns of the block we are taking at compile time (optional)
* \param _PacketAccess allows to enforce aligned loads and stores if set to ForceAligned.
* The default is AsRequested. This parameter is internaly used by Eigen
* \param _PacketAccess allows to enforce aligned loads and stores if set to \b ForceAligned.
* The default is \b AsRequested. This parameter is internaly used by Eigen
* in expressions such as \code mat.block() += other; \endcode and most of
* the time this is the only way it is used.
* \param _DirectAccessStatus \internal used for partial specialization

148
Eigen/src/Core/Dot.h
View File

@ -292,154 +292,6 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<
return ei_sqrt(squaredNorm());
}
/** \returns the \em l2 norm of \c *this using a numerically more stable
* algorithm.
*
* \sa norm(), dot(), squaredNorm(), blueNorm()
*/
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::stableNorm() const
{
return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>());
}
/** \internal Computes ibeta^iexp by binary expansion of iexp,
* exact if ibeta is the machine base */
template<typename T> inline T bexp(int ibeta, int iexp)
{
T tbeta = T(ibeta);
T res = 1.0;
int n = iexp;
if (n<0)
{
n = - n;
tbeta = 1.0/tbeta;
}
for(;;)
{
if ((n % 2)==0)
res = res * tbeta;
n = n/2;
if (n==0) return res;
tbeta = tbeta*tbeta;
}
return res;
}
/** \returns the \em l2 norm of \c *this using the Blue's algorithm.
* A Portable Fortran Program to Find the Euclidean Norm of a Vector,
* ACM TOMS, Vol 4, Issue 1, 1978.
*
* \sa norm(), dot(), squaredNorm(), stableNorm()
*/
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::blueNorm() const
{
static int nmax;
static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
int n;
Scalar ax, abig, amed, asml;
if(nmax <= 0)
{
int nbig, ibeta, it, iemin, iemax, iexp;
Scalar abig, eps;
// This program calculates the machine-dependent constants
// bl, b2, slm, s2m, relerr overfl, nmax
// from the "basic" machine-dependent numbers
// nbig, ibeta, it, iemin, iemax, rbig.
// The following define the basic machine-dependent constants.
// For portability, the PORT subprograms "ilmaeh" and "rlmach"
// are used. For any specific computer, each of the assignment
// statements can be replaced
nbig = std::numeric_limits<int>::max(); // largest integer
ibeta = NumTraits<Scalar>::Base; // base for floating-point numbers
it = NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent
iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent
rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number
// Check the basic machine-dependent constants.
if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
|| (it<=4 && ibeta <= 3 ) || it<2)
{
ei_assert(false && "the algorithm cannot be guaranteed on this computer");
}
iexp = -((1-iemin)/2);
b1 = bexp<Scalar>(ibeta, iexp); // lower boundary of midrange
iexp = (iemax + 1 - it)/2;
b2 = bexp<Scalar>(ibeta,iexp); // upper boundary of midrange
iexp = (2-iemin)/2;
s1m = bexp<Scalar>(ibeta,iexp); // scaling factor for lower range
iexp = - ((iemax+it)/2);
s2m = bexp<Scalar>(ibeta,iexp); // scaling factor for upper range
overfl = rbig*s2m; // overfow boundary for abig
eps = bexp<Scalar>(ibeta, 1-it);
relerr = ei_sqrt(eps); // tolerance for neglecting asml
abig = 1.0/eps - 1.0;
if (Scalar(nbig)>abig) nmax = abig; // largest safe n
else nmax = nbig;
}
n = size();
if(n==0)
return 0;
asml = Scalar(0);
amed = Scalar(0);
abig = Scalar(0);
for(int j=0; j<n; ++j)
{
ax = ei_abs(coeff(j));
if(ax > b2) abig += ei_abs2(ax*s2m);
else if(ax < b1) asml += ei_abs2(ax*s1m);
else amed += ei_abs2(ax);
}
if(abig > Scalar(0))
{
abig = ei_sqrt(abig);
if(abig > overfl)
{
ei_assert(false && "overflow");
return rbig;
}
if(amed > Scalar(0))
{
abig = abig/s2m;
amed = ei_sqrt(amed);
}
else
{
return abig/s2m;
}
}
else if(asml > Scalar(0))
{
if (amed > Scalar(0))
{
abig = ei_sqrt(amed);
amed = ei_sqrt(asml) / s1m;
}
else
{
return ei_sqrt(asml)/s1m;
}
}
else
{
return ei_sqrt(amed);
}
asml = std::min(abig, amed);
abig = std::max(abig, amed);
if(asml <= abig*relerr)
return abig;
else
return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig));
}
/** \returns an expression of the quotient of *this by its own norm.
*
* \only_for_vectors

3
Eigen/src/Core/Functors.h
View File

@ -124,9 +124,6 @@ template<typename Scalar> struct ei_scalar_hypot_op EIGEN_EMPTY_STRUCT {
// typedef typename NumTraits<Scalar>::Real result_type;
EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const
{
// typedef typename NumTraits<T>::Real RealScalar;
// RealScalar _x = ei_abs(x);
// RealScalar _y = ei_abs(y);
Scalar p = std::max(_x, _y);
Scalar q = std::min(_x, _y);
Scalar qp = q/p;

10
Eigen/src/Core/MapBase.h
View File

@ -173,8 +173,14 @@ template<typename Derived> class MapBase
using Base::operator=;
using Base::operator*=;
using Base::operator+=;
using Base::operator-=;
template<typename Lhs,typename Rhs>
Derived& operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProduct>, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other)
{ return Base::operator+=(other); }
template<typename Lhs,typename Rhs>
Derived& operator-=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProduct>, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other)
{ return Base::operator-=(other); }
template<typename OtherDerived>
Derived& operator+=(const MatrixBase<OtherDerived>& other)

10
Eigen/src/Core/Matrix.h
View File

@ -571,10 +571,16 @@ class Matrix
template<typename OtherDerived>
EIGEN_STRONG_INLINE Matrix& _set(const MatrixBase<OtherDerived>& other)
{
_resize_to_match(other);
return Base::operator=(other);
_set_selector(other.derived(), typename ei_meta_if<(int(OtherDerived::Flags) & EvalBeforeAssigningBit), ei_meta_true, ei_meta_false>::ret());
return *this;
}
template<typename OtherDerived>
EIGEN_STRONG_INLINE void _set_selector(const OtherDerived& other, const ei_meta_true&) { _set_noalias(other.eval()); }
template<typename OtherDerived>
EIGEN_STRONG_INLINE void _set_selector(const OtherDerived& other, const ei_meta_false&) { _set_noalias(other); }
/** \internal Like _set() but additionally makes the assumption that no aliasing effect can happen (which
* is the case when creating a new matrix) so one can enforce lazy evaluation.
*

1
Eigen/src/Core/MatrixBase.h
View File

@ -437,6 +437,7 @@ template<typename Derived> class MatrixBase
RealScalar norm() const;
RealScalar stableNorm() const;
RealScalar blueNorm() const;
RealScalar hypotNorm() const;
const PlainMatrixType normalized() const;
void normalize();

8
Eigen/src/Core/NumTraits.h
View File

@ -70,9 +70,7 @@ template<> struct NumTraits<float>
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1,
Base = 2,
Mantissa = 23
MulCost = 1
};
};
@ -85,9 +83,7 @@ template<> struct NumTraits<double>
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1,
Base = 2,
Mantissa = 52
MulCost = 1
};
};

12
Eigen/src/Core/Product.h
View File

@ -287,6 +287,18 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
}
/** replaces \c *this by \c *this * \a other.
*
* \returns a reference to \c *this
*/
template<typename Derived>
template<typename OtherDerived>
inline Derived &
MatrixBase<Derived>::operator*=(const MatrixBase<OtherDerived> &other)
{
return derived() = derived() * other.derived();
}
/***************************************************************************
* Normal product .coeff() implementation (with meta-unrolling)
***************************************************************************/

194
Eigen/src/Core/StableNorm.h Normal file
View File

@ -0,0 +1,194 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_STABLENORM_H
#define EIGEN_STABLENORM_H
template<typename ExpressionType, typename Scalar>
inline void ei_stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
{
Scalar max = bl.cwise().abs().maxCoeff();
if (max>scale)
{
ssq = ssq * ei_abs2(scale/max);
scale = max;
invScale = Scalar(1)/scale;
}
// TODO if the max is much much smaller than the current scale,
// then we can neglect this sub vector
ssq += (bl*invScale).squaredNorm();
}
/** \returns the \em l2 norm of \c *this avoiding underflow and overflow.
* This version use a blockwise two passes algorithm:
* 1 - find the absolute largest coefficient \c s
* 2 - compute \f$ s \Vert \frac{*this}{s} \Vert \f$ in a standard way
*
* For architecture/scalar types supporting vectorization, this version
* is faster than blueNorm(). Otherwise the blueNorm() is much faster.
*
* \sa norm(), blueNorm(), hypotNorm()
*/
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::stableNorm() const
{
const int blockSize = 4096;
RealScalar scale = 0;
RealScalar invScale;
RealScalar ssq = 0; // sum of square
enum {
Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? ForceAligned : AsRequested
};
int n = size();
int bi=0;
if ((int(Flags)&DirectAccessBit) && !(int(Flags)&AlignedBit))
{
bi = ei_alignmentOffset(&const_cast_derived().coeffRef(0), n);
if (bi>0)
ei_stable_norm_kernel(start(bi), ssq, scale, invScale);
}
for (; bi<n; bi+=blockSize)
ei_stable_norm_kernel(VectorBlock<Derived,Dynamic,Alignment>(derived(),bi,std::min(blockSize, n - bi)), ssq, scale, invScale);
return scale * ei_sqrt(ssq);
}
/** \returns the \em l2 norm of \c *this using the Blue's algorithm.
* A Portable Fortran Program to Find the Euclidean Norm of a Vector,
* ACM TOMS, Vol 4, Issue 1, 1978.
*
* For architecture/scalar types without vectorization, this version
* is much faster than stableNorm(). Otherwise the stableNorm() is faster.
*
* \sa norm(), stableNorm(), hypotNorm()
*/
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::blueNorm() const
{
static int nmax = -1;
static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
if(nmax <= 0)
{
int nbig, ibeta, it, iemin, iemax, iexp;
RealScalar abig, eps;
// This program calculates the machine-dependent constants
// bl, b2, slm, s2m, relerr overfl, nmax
// from the "basic" machine-dependent numbers
// nbig, ibeta, it, iemin, iemax, rbig.
// The following define the basic machine-dependent constants.
// For portability, the PORT subprograms "ilmaeh" and "rlmach"
// are used. For any specific computer, each of the assignment
// statements can be replaced
nbig = std::numeric_limits<int>::max(); // largest integer
ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers
it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa
iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent
iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
rbig = std::numeric_limits<RealScalar>::max(); // largest floating-point number
// Check the basic machine-dependent constants.
if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
|| (it<=4 && ibeta <= 3 ) || it<2)
{
ei_assert(false && "the algorithm cannot be guaranteed on this computer");
}
iexp = -((1-iemin)/2);
b1 = std::pow(double(ibeta),iexp); // lower boundary of midrange
iexp = (iemax + 1 - it)/2;
b2 = std::pow(double(ibeta),iexp); // upper boundary of midrange
iexp = (2-iemin)/2;
s1m = std::pow(double(ibeta),iexp); // scaling factor for lower range
iexp = - ((iemax+it)/2);
s2m = std::pow(double(ibeta),iexp); // scaling factor for upper range
overfl = rbig*s2m; // overfow boundary for abig
eps = std::pow(double(ibeta), 1-it);
relerr = ei_sqrt(eps); // tolerance for neglecting asml
abig = 1.0/eps - 1.0;
if (RealScalar(nbig)>abig) nmax = abig; // largest safe n
else nmax = nbig;
}
int n = size();
RealScalar ab2 = b2 / RealScalar(n);
RealScalar asml = RealScalar(0);
RealScalar amed = RealScalar(0);
RealScalar abig = RealScalar(0);
for(int j=0; j<n; ++j)
{
RealScalar ax = ei_abs(coeff(j));
if(ax > ab2) abig += ei_abs2(ax*s2m);
else if(ax < b1) asml += ei_abs2(ax*s1m);
else amed += ei_abs2(ax);
}
if(abig > RealScalar(0))
{
abig = ei_sqrt(abig);
if(abig > overfl)
{
ei_assert(false && "overflow");
return rbig;
}
if(amed > RealScalar(0))
{
abig = abig/s2m;
amed = ei_sqrt(amed);
}
else
return abig/s2m;
}
else if(asml > RealScalar(0))
{
if (amed > RealScalar(0))
{
abig = ei_sqrt(amed);
amed = ei_sqrt(asml) / s1m;
}
else
return ei_sqrt(asml)/s1m;
}
else
return ei_sqrt(amed);
asml = std::min(abig, amed);
abig = std::max(abig, amed);
if(asml <= abig*relerr)
return abig;
else
return abig * ei_sqrt(RealScalar(1) + ei_abs2(asml/abig));
}
/** \returns the \em l2 norm of \c *this avoiding undeflow and overflow.
* This version use a concatenation of hypot() calls, and it is very slow.
*
* \sa norm(), stableNorm()
*/
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::hypotNorm() const
{
return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>());
}
#endif // EIGEN_STABLENORM_H

30
Eigen/src/QR/SelfAdjointEigenSolver.h
View File

@ -383,18 +383,24 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st
int kn1 = (k+1)*n;
#endif
// let's do the product manually to avoid the need of temporaries...
for (int i=0; i<n; ++i)
{
#ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
Scalar matrixQ_i_k = matrixQ[i*n+k];
matrixQ[i*n+k] = c * matrixQ_i_k - s * matrixQ[i*n+k+1];
matrixQ[i*n+k+1] = s * matrixQ_i_k + c * matrixQ[i*n+k+1];
#else
Scalar matrixQ_i_k = matrixQ[i+kn];
matrixQ[i+kn] = c * matrixQ_i_k - s * matrixQ[i+kn1];
matrixQ[i+kn1] = s * matrixQ_i_k + c * matrixQ[i+kn1];
#endif
}
Matrix<Scalar,Dynamic,1> aux = Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn,n);
Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn,n)
= Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn,n) * c - s * Map<Matrix<Scalar,Dynamic,1> >(matrixQ+kn1,n);
Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn1,n)
= Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn1,n) * c - s * aux;//Map<Matrix<Scalar,Dynamic,1> >(matrixQ+kn,n);
// for (int i=0; i<n; ++i)
// {
// #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
// Scalar matrixQ_i_k = matrixQ[i*n+k];
// matrixQ[i*n+k] = c * matrixQ_i_k - s * matrixQ[i*n+k+1];
// matrixQ[i*n+k+1] = s * matrixQ_i_k + c * matrixQ[i*n+k+1];
// #else
// Scalar matrixQ_i_k = matrixQ[i+kn];
// matrixQ[i+kn] = c * matrixQ_i_k - s * matrixQ[i+kn1];
// matrixQ[i+kn1] = s * matrixQ_i_k + c * matrixQ[i+kn1];
// #endif
// }
}
}
}

13
Eigen/src/SVD/SVD.h
View File

@ -125,7 +125,7 @@ template<typename MatrixType> class SVD
{
return (b >= Scalar(0.0) ? ei_abs(a) : -ei_abs(a));
}
protected:
/** \internal */
MatrixUType m_matU;
@ -254,11 +254,14 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
if (g != Scalar(0.0))
{
g = Scalar(1.0)/g;
for (j=l; j<n; j++)
if (m-l)
{
s = A.col(i).end(m-l).dot(A.col(j).end(m-l));
f = (s/A(i,i))*g;
A.col(j).end(m-i) += f * A.col(i).end(m-i);
for (j=l; j<n; j++)
{
s = A.col(i).end(m-l).dot(A.col(j).end(m-l));
f = (s/A(i,i))*g;
A.col(j).end(m-i) += f * A.col(i).end(m-i);
}
}
A.col(i).end(m-i) *= g;
}

344
bench/bench_norm.cpp Normal file
View File

@ -0,0 +1,344 @@
#include <typeinfo>
#include <Eigen/Array>
#include "BenchTimer.h"
using namespace Eigen;
using namespace std;
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v)
{
return v.norm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v)
{
return v.hypotNorm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar blueNorm(const T& v)
{
return v.blueNorm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar lapackNorm(T& v)
{
typedef typename T::Scalar Scalar;
int n = v.size();
Scalar scale = 0;
Scalar ssq = 1;
for (int i=0;i<n;++i)
{
Scalar ax = ei_abs(v.coeff(i));
if (scale >= ax)
{
ssq += ei_abs2(ax/scale);
}
else
{
ssq = Scalar(1) + ssq * ei_abs2(scale/ax);
scale = ax;
}
}
return scale * ei_sqrt(ssq);
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar twopassNorm(T& v)
{
typedef typename T::Scalar Scalar;
Scalar s = v.cwise().abs().maxCoeff();
return s*(v/s).norm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar bl2passNorm(T& v)
{
return v.stableNorm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v)
{
int n =v.size() / 2;
for (int i=0;i<n;++i)
v(i) = v(2*i)*v(2*i) + v(2*i+1)*v(2*i+1);
n = n/2;
while (n>0)
{
for (int i=0;i<n;++i)
v(i) = v(2*i) + v(2*i+1);
n = n/2;
}
return ei_sqrt(v(0));
}
#ifdef EIGEN_VECTORIZE
Packet4f ei_plt(const Packet4f& a, Packet4f& b) { return _mm_cmplt_ps(a,b); }
Packet2d ei_plt(const Packet2d& a, Packet2d& b) { return _mm_cmplt_pd(a,b); }
Packet4f ei_pandnot(const Packet4f& a, Packet4f& b) { return _mm_andnot_ps(a,b); }
Packet2d ei_pandnot(const Packet2d& a, Packet2d& b) { return _mm_andnot_pd(a,b); }
#endif
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
{
#ifndef EIGEN_VECTORIZE
return v.blueNorm();
#else
typedef typename T::Scalar Scalar;
static int nmax = 0;
static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
int n;
if(nmax <= 0)
{
int nbig, ibeta, it, iemin, iemax, iexp;
Scalar abig, eps;
nbig = std::numeric_limits<int>::max(); // largest integer
ibeta = std::numeric_limits<Scalar>::radix; //NumTraits<Scalar>::Base; // base for floating-point numbers
it = std::numeric_limits<Scalar>::digits; //NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent
iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent
rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number
// Check the basic machine-dependent constants.
if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
|| (it<=4 && ibeta <= 3 ) || it<2)
{
ei_assert(false && "the algorithm cannot be guaranteed on this computer");
}
iexp = -((1-iemin)/2);
b1 = std::pow(ibeta, iexp); // lower boundary of midrange
iexp = (iemax + 1 - it)/2;
b2 = std::pow(ibeta,iexp); // upper boundary of midrange
iexp = (2-iemin)/2;
s1m = std::pow(ibeta,iexp); // scaling factor for lower range
iexp = - ((iemax+it)/2);
s2m = std::pow(ibeta,iexp); // scaling factor for upper range
overfl = rbig*s2m; // overfow boundary for abig
eps = std::pow(ibeta, 1-it);
relerr = ei_sqrt(eps); // tolerance for neglecting asml
abig = 1.0/eps - 1.0;
if (Scalar(nbig)>abig) nmax = abig; // largest safe n
else nmax = nbig;
}
typedef typename ei_packet_traits<Scalar>::type Packet;
const int ps = ei_packet_traits<Scalar>::size;
Packet pasml = ei_pset1(Scalar(0));
Packet pamed = ei_pset1(Scalar(0));
Packet pabig = ei_pset1(Scalar(0));
Packet ps2m = ei_pset1(s2m);
Packet ps1m = ei_pset1(s1m);
Packet pb2 = ei_pset1(b2);
Packet pb1 = ei_pset1(b1);
for(int j=0; j<v.size(); j+=ps)
{
Packet ax = ei_pabs(v.template packet<Aligned>(j));
Packet ax_s2m = ei_pmul(ax,ps2m);
Packet ax_s1m = ei_pmul(ax,ps1m);
Packet maskBig = ei_plt(pb2,ax);
Packet maskSml = ei_plt(ax,pb1);
// Packet maskMed = ei_pand(maskSml,maskBig);
// Packet scale = ei_pset1(Scalar(0));
// scale = ei_por(scale, ei_pand(maskBig,ps2m));
// scale = ei_por(scale, ei_pand(maskSml,ps1m));
// scale = ei_por(scale, ei_pandnot(ei_pset1(Scalar(1)),maskMed));
// ax = ei_pmul(ax,scale);
// ax = ei_pmul(ax,ax);
// pabig = ei_padd(pabig, ei_pand(maskBig, ax));
// pasml = ei_padd(pasml, ei_pand(maskSml, ax));
// pamed = ei_padd(pamed, ei_pandnot(ax,maskMed));
pabig = ei_padd(pabig, ei_pand(maskBig, ei_pmul(ax_s2m,ax_s2m)));
pasml = ei_padd(pasml, ei_pand(maskSml, ei_pmul(ax_s1m,ax_s1m)));
pamed = ei_padd(pamed, ei_pandnot(ei_pmul(ax,ax),ei_pand(maskSml,maskBig)));
}
Scalar abig = ei_predux(pabig);
Scalar asml = ei_predux(pasml);
Scalar amed = ei_predux(pamed);
if(abig > Scalar(0))
{
abig = ei_sqrt(abig);
if(abig > overfl)
{
ei_assert(false && "overflow");
return rbig;
}
if(amed > Scalar(0))
{
abig = abig/s2m;
amed = ei_sqrt(amed);
}
else
{
return abig/s2m;
}
}
else if(asml > Scalar(0))
{
if (amed > Scalar(0))
{
abig = ei_sqrt(amed);
amed = ei_sqrt(asml) / s1m;
}
else
{
return ei_sqrt(asml)/s1m;
}
}
else
{
return ei_sqrt(amed);
}
asml = std::min(abig, amed);
abig = std::max(abig, amed);
if(asml <= abig*relerr)
return abig;
else
return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig));
#endif
}
#define BENCH_PERF(NRM) { \
Eigen::BenchTimer tf, td, tcf; tf.reset(); td.reset(); tcf.reset();\
for (int k=0; k<tries; ++k) { \
tf.start(); \
for (int i=0; i<iters; ++i) NRM(vf); \
tf.stop(); \
} \
for (int k=0; k<tries; ++k) { \
td.start(); \
for (int i=0; i<iters; ++i) NRM(vd); \
td.stop(); \
} \
for (int k=0; k<std::max(1,tries/3); ++k) { \
tcf.start(); \
for (int i=0; i<iters; ++i) NRM(vcf); \
tcf.stop(); \
} \
std::cout << #NRM << "\t" << tf.value() << " " << td.value() << " " << tcf.value() << "\n"; \
}
void check_accuracy(double basef, double based, int s)
{
double yf = basef * ei_abs(ei_random<double>());
double yd = based * ei_abs(ei_random<double>());
VectorXf vf = VectorXf::Ones(s) * yf;
VectorXd vd = VectorXd::Ones(s) * yd;
std::cout << "reference\t" << ei_sqrt(double(s))*yf << "\t" << ei_sqrt(double(s))*yd << "\n";
std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\n";
std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\n";
std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\n";
std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\n";
std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\n";
std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\n";
std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\n";
}
void check_accuracy_var(int ef0, int ef1, int ed0, int ed1, int s)
{
VectorXf vf(s);
VectorXd vd(s);
for (int i=0; i<s; ++i)
{
vf[i] = ei_abs(ei_random<double>()) * std::pow(double(10), ei_random<int>(ef0,ef1));
vd[i] = ei_abs(ei_random<double>()) * std::pow(double(10), ei_random<int>(ed0,ed1));
}
//std::cout << "reference\t" << ei_sqrt(double(s))*yf << "\t" << ei_sqrt(double(s))*yd << "\n";
std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\t" << sqsumNorm(vf.cast<long double>()) << "\t" << sqsumNorm(vd.cast<long double>()) << "\n";
std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\t" << hypotNorm(vf.cast<long double>()) << "\t" << hypotNorm(vd.cast<long double>()) << "\n";
std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast<long double>()) << "\t" << lapackNorm(vd.cast<long double>()) << "\n";
std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\t" << twopassNorm(vf.cast<long double>()) << "\t" << twopassNorm(vd.cast<long double>()) << "\n";
// std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\t" << bl2passNorm(vf.cast<long double>()) << "\t" << bl2passNorm(vd.cast<long double>()) << "\n";
}
int main(int argc, char** argv)
{
int tries = 10;
int iters = 100000;
double y = 1.1345743233455785456788e12 * ei_random<double>();
VectorXf v = VectorXf::Ones(1024) * y;
// return 0;
int s = 10000;
double basef_ok = 1.1345743233455785456788e15;
double based_ok = 1.1345743233455785456788e95;
double basef_under = 1.1345743233455785456788e-27;
double based_under = 1.1345743233455785456788e-303;
double basef_over = 1.1345743233455785456788e+27;
double based_over = 1.1345743233455785456788e+302;
std::cout.precision(20);
std::cerr << "\nNo under/overflow:\n";
check_accuracy(basef_ok, based_ok, s);
std::cerr << "\nUnderflow:\n";
check_accuracy(basef_under, based_under, s);
std::cerr << "\nOverflow:\n";
check_accuracy(basef_over, based_over, s);
std::cerr << "\nVarying (over):\n";
for (int k=0; k<1; ++k)
{
check_accuracy_var(20,27,190,302,s);
std::cout << "\n";
}
std::cerr << "\nVarying (under):\n";
for (int k=0; k<1; ++k)
{
check_accuracy_var(-27,20,-302,-190,s);
std::cout << "\n";
}
std::cout.precision(4);
std::cerr << "Performance (out of cache):\n";
{
int iters = 1;
VectorXf vf = VectorXf::Random(1024*1024*32) * y;
VectorXd vd = VectorXd::Random(1024*1024*32) * y;
VectorXcf vcf = VectorXcf::Random(1024*1024*32) * y;
BENCH_PERF(sqsumNorm);
BENCH_PERF(blueNorm);
// BENCH_PERF(pblueNorm);
// BENCH_PERF(lapackNorm);
// BENCH_PERF(hypotNorm);
// BENCH_PERF(twopassNorm);
BENCH_PERF(bl2passNorm);
}
std::cerr << "\nPerformance (in cache):\n";
{
int iters = 100000;
VectorXf vf = VectorXf::Random(512) * y;
VectorXd vd = VectorXd::Random(512) * y;
VectorXcf vcf = VectorXcf::Random(512) * y;
BENCH_PERF(sqsumNorm);
BENCH_PERF(blueNorm);
// BENCH_PERF(pblueNorm);
// BENCH_PERF(lapackNorm);
// BENCH_PERF(hypotNorm);
// BENCH_PERF(twopassNorm);
BENCH_PERF(bl2passNorm);
}
}

4
test/adjoint.cpp
View File

@ -76,8 +76,8 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
{
VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1));
VERIFY_IS_APPROX(v1.norm(), v1.stableNorm());
// NOTE disabled because it currently compiles for float and double only
// VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm());
VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm());
VERIFY_IS_APPROX(v1.hypotNorm(), v1.stableNorm());
}
// check compatibility of dot and adjoint

4
test/mixingtypes.cpp
View File

@ -82,7 +82,7 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType)
void test_mixingtypes()
{
// check that our operator new is indeed called:
CALL_SUBTEST(mixingtypes<3>());
CALL_SUBTEST(mixingtypes<4>());
CALL_SUBTEST(mixingtypes<3>(3));
CALL_SUBTEST(mixingtypes<4>(4));
CALL_SUBTEST(mixingtypes<Dynamic>(20));
}

6
test/product_large.cpp
View File

@ -42,4 +42,10 @@ void test_product_large()
m = (v+v).asDiagonal() * m;
VERIFY_IS_APPROX(m, MatrixXf::Constant(N,3,2));
}
{
// test deferred resizing in Matrix::operator=
MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a;
VERIFY_IS_APPROX((a = a * b), (c * b).eval());
}
}