mirror of
https://gitlab.com/libeigen/eigen.git
synced 2026-04-10 11:34:33 +08:00
Clean up informal language, vague TODOs, and dead code in comments
libeigen/eigen!2191 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
@@ -225,7 +225,7 @@ class LDLT : public SolverBase<LDLT<MatrixType_, UpLo_> > {
|
||||
|
||||
/** \returns the internal LDLT decomposition matrix
|
||||
*
|
||||
* TODO: document the storage layout
|
||||
* TODO: document the storage layout.
|
||||
*/
|
||||
inline const MatrixType& matrixLDLT() const {
|
||||
eigen_assert(m_isInitialized && "LDLT is not initialized.");
|
||||
|
||||
@@ -178,9 +178,6 @@ class ArrayBase : public DenseBase<Derived> {
|
||||
return MatrixWrapper<const Derived>(derived());
|
||||
}
|
||||
|
||||
// template<typename Dest>
|
||||
// inline void evalTo(Dest& dst) const { dst = matrix(); }
|
||||
|
||||
protected:
|
||||
EIGEN_DEFAULT_COPY_CONSTRUCTOR(ArrayBase)
|
||||
EIGEN_DEFAULT_EMPTY_CONSTRUCTOR_AND_DESTRUCTOR(ArrayBase)
|
||||
|
||||
@@ -63,7 +63,7 @@ struct copy_using_evaluator_traits {
|
||||
static constexpr int RestrictedLinearSize = min_size_prefer_fixed(MaxSizeAtCompileTime, MaxPacketSize);
|
||||
static constexpr int OuterStride = outer_stride_at_compile_time<Dst>::ret;
|
||||
|
||||
// TODO distinguish between linear traversal and inner-traversals
|
||||
// TODO: distinguish between linear traversal and inner-traversal packet types.
|
||||
using LinearPacketType = typename find_best_packet<DstScalar, RestrictedLinearSize>::type;
|
||||
using InnerPacketType = typename find_best_packet<DstScalar, RestrictedInnerSize>::type;
|
||||
|
||||
@@ -1016,7 +1016,7 @@ struct Assignment<DstXprType, CwiseNullaryOp<scalar_zero_op<typename DstXprType:
|
||||
};
|
||||
|
||||
// Generic assignment through evalTo.
|
||||
// TODO: not sure we have to keep that one, but it helps porting current code to new evaluator mechanism.
|
||||
// TODO: evaluate whether this generic evalTo-based assignment path is still needed.
|
||||
// Note that the last template argument "Weak" is needed to make it possible to perform
|
||||
// both partial specialization+SFINAE without ambiguous specialization
|
||||
template <typename DstXprType, typename SrcXprType, typename Functor, typename Weak>
|
||||
|
||||
@@ -98,7 +98,7 @@ class CwiseBinaryOp : public CwiseBinaryOpImpl<BinaryOp, LhsType, RhsType,
|
||||
typedef std::remove_reference_t<RhsNested> RhsNested_;
|
||||
|
||||
#if EIGEN_COMP_MSVC
|
||||
// Required for Visual Studio or the Copy constructor will probably not get inlined!
|
||||
// Required for Visual Studio, which may fail to inline the copy constructor otherwise.
|
||||
EIGEN_STRONG_INLINE CwiseBinaryOp(const CwiseBinaryOp<BinaryOp, LhsType, RhsType>&) = default;
|
||||
#endif
|
||||
|
||||
|
||||
@@ -431,8 +431,7 @@ class DenseBase
|
||||
|
||||
// By default, the fastest version with undefined NaN propagation semantics is
|
||||
// used.
|
||||
// TODO(rmlarsen): Replace with default template argument when we move to
|
||||
// c++11 or beyond.
|
||||
// TODO(rmlarsen): Replace with default template argument (C++14 is now the minimum standard).
|
||||
EIGEN_DEVICE_FUNC inline typename internal::traits<Derived>::Scalar minCoeff() const {
|
||||
return minCoeff<PropagateFast>();
|
||||
}
|
||||
@@ -449,7 +448,7 @@ class DenseBase
|
||||
template <int NaNPropagation, typename IndexType>
|
||||
EIGEN_DEVICE_FUNC typename internal::traits<Derived>::Scalar maxCoeff(IndexType* index) const;
|
||||
|
||||
// TODO(rmlarsen): Replace these methods with a default template argument.
|
||||
// TODO(rmlarsen): Replace these methods with a default template argument (C++14 is now the minimum standard).
|
||||
template <typename IndexType>
|
||||
EIGEN_DEVICE_FUNC inline typename internal::traits<Derived>::Scalar minCoeff(IndexType* row, IndexType* col) const {
|
||||
return minCoeff<PropagateFast>(row, col);
|
||||
@@ -580,12 +579,12 @@ class DenseBase
|
||||
#else
|
||||
typedef std::conditional_t<(Flags & DirectAccessBit) == DirectAccessBit,
|
||||
internal::pointer_based_stl_iterator<Derived>,
|
||||
internal::generic_randaccess_stl_iterator<Derived> >
|
||||
internal::generic_randaccess_stl_iterator<Derived>>
|
||||
iterator_type;
|
||||
|
||||
typedef std::conditional_t<(Flags & DirectAccessBit) == DirectAccessBit,
|
||||
internal::pointer_based_stl_iterator<const Derived>,
|
||||
internal::generic_randaccess_stl_iterator<const Derived> >
|
||||
internal::generic_randaccess_stl_iterator<const Derived>>
|
||||
const_iterator_type;
|
||||
|
||||
// Stl-style iterators are supported only for vectors.
|
||||
|
||||
@@ -89,7 +89,7 @@ struct product_type {
|
||||
/* The following allows to select the kind of product at compile time
|
||||
* based on the three dimensions of the product.
|
||||
* This is a compile time mapping from {1,Small,Large}^3 -> {product types} */
|
||||
// FIXME I'm not sure the current mapping is the ideal one.
|
||||
// FIXME: the current compile-time product-type mapping may not be optimal.
|
||||
template <int M, int N>
|
||||
struct product_type_selector<M, N, 1> {
|
||||
enum { ret = OuterProduct };
|
||||
@@ -193,12 +193,11 @@ struct product_type_selector<Large, Large, Small> {
|
||||
* Implementation of Inner Vector Vector Product
|
||||
***********************************************************************/
|
||||
|
||||
// FIXME : maybe the "inner product" could return a Scalar
|
||||
// instead of a 1x1 matrix ??
|
||||
// Pro: more natural for the user
|
||||
// Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix
|
||||
// product ends up to a row-vector times col-vector product... To tackle this use
|
||||
// case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x);
|
||||
// FIXME: consider returning a Scalar instead of a 1x1 matrix for inner products.
|
||||
// Pro: more natural for the user.
|
||||
// Con: in a meta-unrolled algorithm a matrix-matrix product may reduce to a
|
||||
// row-vector times column-vector product. To handle this, we could specialize
|
||||
// Block<MatrixType,1,1> with operator=(Scalar x).
|
||||
|
||||
/***********************************************************************
|
||||
* Implementation of Outer Vector Vector Product
|
||||
|
||||
@@ -1329,9 +1329,7 @@ EIGEN_DEVICE_FUNC inline typename unpacket_traits<Packet>::type predux_max(const
|
||||
/** \internal \returns true if all coeffs of \a a means "true"
|
||||
* It is supposed to be called on values returned by pcmp_*.
|
||||
*/
|
||||
// not needed yet
|
||||
// template<typename Packet> EIGEN_DEVICE_FUNC inline bool predux_all(const Packet& a)
|
||||
// { return bool(a); }
|
||||
// TODO: implement predux_all when needed.
|
||||
|
||||
/** \internal \returns true if any coeffs of \a a means "true"
|
||||
* It is supposed to be called on values returned by pcmp_*.
|
||||
|
||||
@@ -11,7 +11,7 @@
|
||||
#ifndef EIGEN_MATHFUNCTIONS_H
|
||||
#define EIGEN_MATHFUNCTIONS_H
|
||||
|
||||
// TODO this should better be moved to NumTraits
|
||||
// TODO: consider moving these constants to NumTraits.
|
||||
// Source: WolframAlpha
|
||||
#define EIGEN_PI 3.141592653589793238462643383279502884197169399375105820974944592307816406L
|
||||
#define EIGEN_LOG2E 1.442695040888963407359924681001892137426645954152985934135449406931109219L
|
||||
@@ -390,7 +390,7 @@ struct cast_impl<OldType, NewType,
|
||||
}
|
||||
};
|
||||
|
||||
// here, for once, we're plainly returning NewType: we don't want cast to do weird things.
|
||||
// Returns NewType directly to avoid unintended intermediate conversions.
|
||||
|
||||
template <typename OldType, typename NewType>
|
||||
EIGEN_DEVICE_FUNC inline NewType cast(const OldType& x) {
|
||||
@@ -832,8 +832,8 @@ EIGEN_DEVICE_FUNC std::enable_if_t<(std::numeric_limits<T>::has_infinity && !Num
|
||||
|
||||
template <typename T>
|
||||
EIGEN_DEVICE_FUNC
|
||||
std::enable_if_t<!(std::numeric_limits<T>::has_quiet_NaN || std::numeric_limits<T>::has_signaling_NaN), bool>
|
||||
isnan_impl(const T&) {
|
||||
std::enable_if_t<!(std::numeric_limits<T>::has_quiet_NaN || std::numeric_limits<T>::has_signaling_NaN), bool>
|
||||
isnan_impl(const T&) {
|
||||
return false;
|
||||
}
|
||||
|
||||
@@ -1450,9 +1450,9 @@ EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
EIGEN_DEVICE_FUNC
|
||||
EIGEN_ALWAYS_INLINE std::enable_if_t<!(NumTraits<T>::IsSigned || NumTraits<T>::IsComplex), typename NumTraits<T>::Real>
|
||||
abs(const T& x) {
|
||||
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
|
||||
std::enable_if_t<!(NumTraits<T>::IsSigned || NumTraits<T>::IsComplex), typename NumTraits<T>::Real>
|
||||
abs(const T& x) {
|
||||
return x;
|
||||
}
|
||||
|
||||
|
||||
@@ -182,7 +182,7 @@ struct Assignment<DstXprType,
|
||||
|
||||
//----------------------------------------
|
||||
// Catch "Dense ?= xpr + Product<>" expression to save one temporary
|
||||
// FIXME we could probably enable these rules for any product, i.e., not only Dense and DefaultProduct
|
||||
// FIXME: consider enabling these rules for all product types, not only Dense and DefaultProduct.
|
||||
|
||||
template <typename OtherXpr, typename Lhs, typename Rhs>
|
||||
struct evaluator_assume_aliasing<
|
||||
@@ -1158,7 +1158,7 @@ struct generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, Pr
|
||||
* Products with transpositions matrices
|
||||
***************************************************************************/
|
||||
|
||||
// FIXME could we unify Transpositions and Permutation into a single "shape"??
|
||||
// FIXME: consider unifying Transpositions and Permutation into a single shape.
|
||||
|
||||
/** \internal
|
||||
* \class transposition_matrix_product
|
||||
|
||||
@@ -43,7 +43,7 @@ struct traits<Ref<PlainObjectType_, Options_, StrideType_> >
|
||||
OuterStrideMatch = IsVectorAtCompileTime || int(OuterStrideAtCompileTime) == int(Dynamic) ||
|
||||
int(OuterStrideAtCompileTime) == int(Derived::OuterStrideAtCompileTime),
|
||||
// NOTE, this indirection of evaluator<Derived>::Alignment is needed
|
||||
// to workaround a very strange bug in MSVC related to the instantiation
|
||||
// to work around an MSVC bug related to the instantiation
|
||||
// of has_*ary_operator in evaluator<CwiseNullaryOp>.
|
||||
// This line is surprisingly very sensitive. For instance, simply adding parenthesis
|
||||
// as "DerivedAlignment = (int(evaluator<Derived>::Alignment))," will make MSVC fail...
|
||||
|
||||
@@ -40,8 +40,7 @@ inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& sc
|
||||
scale = maxCoeff;
|
||||
}
|
||||
|
||||
// TODO if the maxCoeff is much much smaller than the current scale,
|
||||
// then we can neglect this sub vector
|
||||
// TODO: skip sub-vector when maxCoeff << current scale.
|
||||
if (scale > Scalar(0)) // if scale==0, then bl is 0
|
||||
ssq += (bl * invScale).squaredNorm();
|
||||
}
|
||||
|
||||
@@ -1407,12 +1407,12 @@ EIGEN_STRONG_INLINE Packet8l preverse(const Packet8l& a) {
|
||||
|
||||
template <>
|
||||
EIGEN_STRONG_INLINE Packet16f pabs(const Packet16f& a) {
|
||||
// _mm512_abs_ps intrinsic not found, so hack around it
|
||||
// _mm512_abs_ps intrinsic not found, so implement via bitwise AND with sign-bit mask.
|
||||
return _mm512_castsi512_ps(_mm512_and_si512(_mm512_castps_si512(a), _mm512_set1_epi32(0x7fffffff)));
|
||||
}
|
||||
template <>
|
||||
EIGEN_STRONG_INLINE Packet8d pabs(const Packet8d& a) {
|
||||
// _mm512_abs_ps intrinsic not found, so hack around it
|
||||
// _mm512_abs_pd intrinsic not found, so implement via bitwise AND with sign-bit mask.
|
||||
return _mm512_castsi512_pd(_mm512_and_si512(_mm512_castpd_si512(a), _mm512_set1_epi64(0x7fffffffffffffff)));
|
||||
}
|
||||
template <>
|
||||
|
||||
@@ -55,7 +55,7 @@ EIGEN_STRONG_INLINE int64_t predux(const Packet8l& a) {
|
||||
// MSVC's _mm512_reduce_mul_epi64 is borked, at least up to and including 1939.
|
||||
// alignas(64) int64_t data[] = { 1,1,-1,-1,1,-1,-1,-1 };
|
||||
// int64_t out = _mm512_reduce_mul_epi64(_mm512_load_epi64(data));
|
||||
// produces garbage: 4294967295. It seems to happen whenever the output is supposed to be negative.
|
||||
// produces garbage: 4294967295. This occurs when the result should be negative.
|
||||
// Fall back to a manual approach:
|
||||
template <>
|
||||
EIGEN_STRONG_INLINE int64_t predux_mul(const Packet8l& a) {
|
||||
|
||||
@@ -294,7 +294,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_impl_double(cons
|
||||
Packet x2 = pmul(x, x);
|
||||
Packet x3 = pmul(x2, x);
|
||||
|
||||
// Evaluate the polynomial approximant , probably to improve instruction-level parallelism.
|
||||
// Evaluate the polynomial in factored form for better instruction-level parallelism.
|
||||
// y = x - 0.5*x^2 + x^3 * polevl( x, P, 5 ) / p1evl( x, Q, 5 ) );
|
||||
Packet y, y1, y_;
|
||||
y = pmadd(cst_cephes_log_p0, x, cst_cephes_log_p1);
|
||||
@@ -1861,8 +1861,8 @@ struct accurate_log2<double> {
|
||||
// This function implements the non-trivial case of pow(x,y) where x is
|
||||
// positive and y is (possibly) non-integer.
|
||||
// Formally, pow(x,y) = exp2(y * log2(x)), where exp2(x) is shorthand for 2^x.
|
||||
// TODO(rmlarsen): We should probably add this as a packet up 'ppow', to make it
|
||||
// easier to specialize or turn off for specific types and/or backends.x
|
||||
// TODO(rmlarsen): We should probably add this as a packet op 'ppow', to make it
|
||||
// easier to specialize or turn off for specific types and/or backends.
|
||||
template <typename Packet>
|
||||
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_pow_impl(const Packet& x, const Packet& y) {
|
||||
typedef typename unpacket_traits<Packet>::type Scalar;
|
||||
|
||||
@@ -249,8 +249,8 @@ void evaluateProductBlockingSizesHeuristic(Index& k, Index& m, Index& n, Index n
|
||||
|
||||
// Here, nc is chosen such that a block of kc x nc of the rhs fit within half of L2.
|
||||
// The second half is implicitly reserved to access the result and lhs coefficients.
|
||||
// When k<max_kc, then nc can arbitrarily growth. In practice, it seems to be fruitful
|
||||
// to limit this growth: we bound nc to growth by a factor x1.5.
|
||||
// When k<max_kc, then nc can grow without bound. In practice, it seems to be fruitful
|
||||
// to limit this growth: we bound nc growth to a factor of 1.5x.
|
||||
// However, if the entire lhs block fit within L1, then we are not going to block on the rows at all,
|
||||
// and it becomes fruitful to keep the packed rhs blocks in L1 if there is enough remaining space.
|
||||
Index max_nc;
|
||||
@@ -587,8 +587,7 @@ class gebp_traits<std::complex<RealScalar>, RealScalar, ConjLhs_, false, Arch, P
|
||||
}
|
||||
|
||||
EIGEN_STRONG_INLINE void loadRhsQuad_impl(const RhsScalar* b, RhsPacket& dest, const true_type&) const {
|
||||
// FIXME we can do better!
|
||||
// what we want here is a ploadheight
|
||||
// FIXME: replace with a dedicated ploadheight operation for more efficient quad loading.
|
||||
RhsScalar tmp[4] = {b[0], b[0], b[1], b[1]};
|
||||
dest = ploadquad<RhsPacket>(tmp);
|
||||
}
|
||||
@@ -669,7 +668,7 @@ DoublePacket<typename unpacket_traits<Packet>::half> predux_half(
|
||||
const DoublePacket<Packet>& a,
|
||||
std::enable_if_t<unpacket_traits<Packet>::size >= 16 &&
|
||||
!NumTraits<typename unpacket_traits<Packet>::type>::IsComplex>* = 0) {
|
||||
// yes, that's pretty hackish :(
|
||||
// Workaround: reduce real packets to half size by reinterpreting as complex.
|
||||
DoublePacket<typename unpacket_traits<Packet>::half> res;
|
||||
typedef std::complex<typename unpacket_traits<Packet>::type> Cplx;
|
||||
typedef typename packet_traits<Cplx>::type CplxPacket;
|
||||
@@ -689,7 +688,7 @@ void loadQuadToDoublePacket(const Scalar* b, DoublePacket<RealPacket>& dest,
|
||||
template <typename Scalar, typename RealPacket>
|
||||
void loadQuadToDoublePacket(const Scalar* b, DoublePacket<RealPacket>& dest,
|
||||
std::enable_if_t<unpacket_traits<RealPacket>::size == 16>* = 0) {
|
||||
// yes, that's pretty hackish too :(
|
||||
// Workaround: load quad elements by reinterpreting real packets as complex.
|
||||
typedef typename NumTraits<Scalar>::Real RealScalar;
|
||||
RealScalar r[4] = {numext::real(b[0]), numext::real(b[0]), numext::real(b[1]), numext::real(b[1])};
|
||||
RealScalar i[4] = {numext::imag(b[0]), numext::imag(b[0]), numext::imag(b[1]), numext::imag(b[1])};
|
||||
|
||||
@@ -383,8 +383,8 @@ struct generic_product_impl<Lhs, Rhs, DenseShape, DenseShape, GemmProduct>
|
||||
// to determine the following heuristic.
|
||||
// EIGEN_GEMM_TO_COEFFBASED_THRESHOLD is typically defined to 20 in GeneralProduct.h,
|
||||
// unless it has been specialized by the user or for a given architecture.
|
||||
// Note that the condition rhs.rows()>0 was required because lazy product is (was?) not happy with empty inputs.
|
||||
// I'm not sure it is still required.
|
||||
// Note that the condition rhs.rows()>0 was required because lazy product did not handle empty inputs
|
||||
// correctly. It is unclear whether this guard is still necessary.
|
||||
if ((rhs.rows() + dst.rows() + dst.cols()) < EIGEN_GEMM_TO_COEFFBASED_THRESHOLD && rhs.rows() > 0)
|
||||
lazyproduct::eval_dynamic(dst, lhs, rhs, internal::assign_op<typename Dst::Scalar, Scalar>());
|
||||
else {
|
||||
|
||||
@@ -182,7 +182,7 @@ EIGEN_STRONG_INLINE void parallelize_gemm(const Functor& func, Index rows, Index
|
||||
|
||||
// compute the maximal number of threads from the total amount of work:
|
||||
double work = static_cast<double>(rows) * static_cast<double>(cols) * static_cast<double>(depth);
|
||||
double kMinTaskSize = 50000; // FIXME improve this heuristic.
|
||||
double kMinTaskSize = 50000; // FIXME: tune this minimum task-size heuristic based on architecture and scalar type.
|
||||
pb_max_threads = std::max<Index>(1, std::min<Index>(pb_max_threads, static_cast<Index>(work / kMinTaskSize)));
|
||||
|
||||
// compute the number of threads we are going to use
|
||||
|
||||
@@ -212,7 +212,7 @@ struct trmv_selector<Mode, ColMajor> {
|
||||
ResScalar actualAlpha = alpha * lhs_alpha * rhs_alpha;
|
||||
|
||||
// FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
|
||||
// on, the other hand it is good for the cache to pack the vector anyways...
|
||||
// On the other hand, it is good for the cache to pack the vector anyways...
|
||||
constexpr bool EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime == 1;
|
||||
constexpr bool ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex);
|
||||
constexpr bool MightCannotUseDest = (Dest::InnerStrideAtCompileTime != 1) || ComplexByReal;
|
||||
|
||||
@@ -357,7 +357,7 @@
|
||||
// notice that since these are C headers, the extern "C" is theoretically needed anyways.
|
||||
extern "C" {
|
||||
// In theory we should only include immintrin.h and not the other *mmintrin.h header files directly.
|
||||
// Doing so triggers some issues with ICC. However old gcc versions seems to not have this file, thus:
|
||||
// Doing so triggers some issues with ICC. However old gcc versions may not have this file, thus:
|
||||
#if EIGEN_COMP_ICC >= 1110 || EIGEN_COMP_EMSCRIPTEN
|
||||
#include <immintrin.h>
|
||||
#else
|
||||
@@ -388,7 +388,7 @@ extern "C" {
|
||||
#define EIGEN_VECTORIZE_VSX 1
|
||||
#define EIGEN_VECTORIZE_FMA
|
||||
#include <altivec.h>
|
||||
// We need to #undef all these ugly tokens defined in <altivec.h>
|
||||
// We need to #undef macros defined by <altivec.h> that conflict with standard C++ names.
|
||||
// => use __vector instead of vector
|
||||
#undef bool
|
||||
#undef vector
|
||||
@@ -400,7 +400,7 @@ extern "C" {
|
||||
#define EIGEN_VECTORIZE_ALTIVEC
|
||||
#define EIGEN_VECTORIZE_FMA
|
||||
#include <altivec.h>
|
||||
// We need to #undef all these ugly tokens defined in <altivec.h>
|
||||
// We need to #undef macros defined by <altivec.h> that conflict with standard C++ names.
|
||||
// => use __vector instead of vector
|
||||
#undef bool
|
||||
#undef vector
|
||||
|
||||
@@ -1139,7 +1139,7 @@ EIGEN_DEVICE_FUNC constexpr void ignore_unused_variable(const T&) {}
|
||||
|
||||
#if EIGEN_COMP_MSVC
|
||||
// NOTE MSVC often gives C4127 warnings with compiletime if statements. See bug 1362.
|
||||
// This workaround is ugly, but it does the job.
|
||||
// This workaround suppresses MSVC C4127 warnings for compile-time conditionals.
|
||||
#define EIGEN_CONST_CONDITIONAL(cond) (void)0, cond
|
||||
#else
|
||||
#define EIGEN_CONST_CONDITIONAL(cond) cond
|
||||
|
||||
@@ -435,9 +435,8 @@ struct greater_equal_zero_op {
|
||||
|
||||
/* reductions for lists */
|
||||
|
||||
// using auto -> return value spec makes ICC 13.0 and 13.1 crash here, so we have to hack it
|
||||
// together in front... (13.0 doesn't work with array_prod/array_reduce/... anyway, but 13.1
|
||||
// does...
|
||||
// Using auto -> return value spec makes ICC 13.0 and 13.1 crash here,
|
||||
// so the return type is specified explicitly using decltype.
|
||||
template <typename... Ts>
|
||||
EIGEN_DEVICE_FUNC constexpr decltype(reduce<product_op, Ts...>::run((*((Ts*)0))...)) arg_prod(Ts... ts) {
|
||||
return reduce<product_op, Ts...>::run(ts...);
|
||||
|
||||
@@ -277,8 +277,8 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
|
||||
using std::abs;
|
||||
if ((iter == 10 || iter == 20) && iu > 1) {
|
||||
// exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
|
||||
return ComplexSchur<MatrixType>::ComplexScalar(
|
||||
abs(numext::real(m_matT.coeff(iu, iu - 1))) + abs(numext::real(m_matT.coeff(iu - 1, iu - 2))));
|
||||
return ComplexSchur<MatrixType>::ComplexScalar(abs(numext::real(m_matT.coeff(iu, iu - 1))) +
|
||||
abs(numext::real(m_matT.coeff(iu - 1, iu - 2))));
|
||||
}
|
||||
|
||||
// compute the shift as one of the eigenvalues of t, the 2x2
|
||||
@@ -363,7 +363,7 @@ struct complex_schur_reduce_to_hessenberg<MatrixType, false> {
|
||||
_this.m_hess.compute(matrix);
|
||||
_this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
|
||||
if (computeU) {
|
||||
// This may cause an allocation which seems to be avoidable
|
||||
// TODO: this temporary allocation could potentially be avoided.
|
||||
MatrixType Q = _this.m_hess.matrixQ();
|
||||
_this.m_matU = Q.template cast<ComplexScalar>();
|
||||
}
|
||||
|
||||
@@ -317,7 +317,7 @@ namespace internal {
|
||||
* HessenbergDecomposition class until the it is assigned or evaluated for
|
||||
* some other reason (the reference should remain valid during the life time
|
||||
* of this object). This class is the return type of
|
||||
* HessenbergDecomposition::matrixH(); there is probably no other use for this
|
||||
* HessenbergDecomposition::matrixH(); there is no other intended use for this
|
||||
* class.
|
||||
*/
|
||||
template <typename MatrixType>
|
||||
|
||||
@@ -111,8 +111,7 @@ template <typename Derived>
|
||||
inline typename MatrixBase<Derived>::RealScalar MatrixBase<Derived>::operatorNorm() const {
|
||||
using std::sqrt;
|
||||
typename Derived::PlainObject m_eval(derived());
|
||||
// FIXME if it is really guaranteed that the eigenvalues are already sorted,
|
||||
// then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
|
||||
// FIXME: if eigenvalues are guaranteed to be sorted, comparing the first and last is sufficient.
|
||||
return sqrt((m_eval * m_eval.adjoint()).eval().template selfadjointView<Lower>().eigenvalues().maxCoeff());
|
||||
}
|
||||
|
||||
|
||||
@@ -343,7 +343,7 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa
|
||||
template <typename MatrixType>
|
||||
inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() {
|
||||
const Index size = m_matT.cols();
|
||||
// FIXME to be efficient the following would requires a triangular reduxion code
|
||||
// FIXME: a triangular reduction would be more efficient here.
|
||||
// Scalar norm = m_matT.upper().cwiseAbs().sum()
|
||||
// + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
|
||||
Scalar norm(0);
|
||||
|
||||
@@ -548,8 +548,7 @@ EIGEN_DEVICE_FUNC ComputationInfo computeFromTridiagonal_impl(DiagType& diag, Su
|
||||
info = NoConvergence;
|
||||
|
||||
// Sort eigenvalues and corresponding vectors.
|
||||
// TODO make the sort optional ?
|
||||
// TODO use a better sort algorithm !!
|
||||
// TODO: make the sort optional and use a more efficient sorting algorithm.
|
||||
if (info == Success) {
|
||||
for (Index i = 0; i < n - 1; ++i) {
|
||||
Index k;
|
||||
@@ -653,12 +652,12 @@ struct direct_selfadjoint_eigenvalues<SolverType, 3, false> {
|
||||
|
||||
// Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
|
||||
Scalar shift = mat.trace() / Scalar(3);
|
||||
// TODO Avoid this copy. Currently it is necessary to suppress bogus values when determining maxCoeff and for
|
||||
// computing the eigenvectors later
|
||||
// TODO: avoid this copy. Currently necessary to suppress bogus values when determining maxCoeff and for
|
||||
// computing the eigenvectors later.
|
||||
MatrixType scaledMat = mat.template selfadjointView<Lower>();
|
||||
scaledMat.diagonal().array() -= shift;
|
||||
Scalar scale = scaledMat.cwiseAbs().maxCoeff();
|
||||
if (scale > 0) scaledMat /= scale; // TODO for scale==0 we could save the remaining operations
|
||||
if (scale > 0) scaledMat /= scale; // TODO: skip remaining operations when scale==0.
|
||||
|
||||
// compute the eigenvalues
|
||||
computeRoots(scaledMat, eivals);
|
||||
|
||||
@@ -247,7 +247,7 @@ struct homogeneous_left_product_impl<Homogeneous<MatrixType, Vertical>, Lhs>
|
||||
|
||||
template <typename Dest>
|
||||
EIGEN_DEVICE_FUNC void evalTo(Dest& dst) const {
|
||||
// FIXME investigate how to allow lazy evaluation of this product when possible
|
||||
// FIXME: investigate how to allow lazy evaluation of this product when possible.
|
||||
dst = Block < const LhsMatrixTypeNested, LhsMatrixTypeNested::RowsAtCompileTime,
|
||||
LhsMatrixTypeNested::ColsAtCompileTime == Dynamic
|
||||
? Dynamic
|
||||
@@ -278,7 +278,7 @@ struct homogeneous_right_product_impl<Homogeneous<MatrixType, Horizontal>, Rhs>
|
||||
|
||||
template <typename Dest>
|
||||
EIGEN_DEVICE_FUNC void evalTo(Dest& dst) const {
|
||||
// FIXME investigate how to allow lazy evaluation of this product when possible
|
||||
// FIXME: investigate how to allow lazy evaluation of this product when possible.
|
||||
dst = m_lhs * Block < const RhsNested,
|
||||
RhsNested::RowsAtCompileTime == Dynamic ? Dynamic : RhsNested::RowsAtCompileTime - 1,
|
||||
RhsNested::ColsAtCompileTime > (m_rhs, 0, 0, m_rhs.rows() - 1, m_rhs.cols());
|
||||
@@ -392,8 +392,7 @@ struct generic_product_impl<Lhs, Homogeneous<RhsArg, Vertical>, DenseShape, Homo
|
||||
}
|
||||
};
|
||||
|
||||
// TODO: the following specialization is to address a regression from 3.2 to 3.3
|
||||
// In the future, this path should be optimized.
|
||||
// TODO: this specialization addresses a performance regression from 3.2 to 3.3; optimize this path.
|
||||
template <typename Lhs, typename RhsArg, int ProductTag>
|
||||
struct generic_product_impl<Lhs, Homogeneous<RhsArg, Vertical>, TriangularShape, HomogeneousShape, ProductTag> {
|
||||
template <typename Dest>
|
||||
|
||||
@@ -111,7 +111,7 @@ class Hyperplane {
|
||||
* If the dimension of the ambient space is greater than 2, then there isn't uniqueness,
|
||||
* so an arbitrary choice is made.
|
||||
*/
|
||||
// FIXME to be consistent with the rest this could be implemented as a static Through function ??
|
||||
// FIXME: for consistency, consider implementing as a static Through function.
|
||||
EIGEN_DEVICE_FUNC explicit Hyperplane(const ParametrizedLine<Scalar, AmbientDimAtCompileTime>& parametrized) {
|
||||
normal() = parametrized.direction().unitOrthogonal();
|
||||
offset() = -parametrized.origin().dot(normal());
|
||||
|
||||
@@ -790,7 +790,7 @@ EIGEN_DEVICE_FUNC Quaternion<Scalar, Options> Quaternion<Scalar, Options>::FromT
|
||||
template <class Derived>
|
||||
EIGEN_DEVICE_FUNC inline Quaternion<typename internal::traits<Derived>::Scalar> QuaternionBase<Derived>::inverse()
|
||||
const {
|
||||
// FIXME should this function be called multiplicativeInverse and conjugate() be called inverse() or opposite() ??
|
||||
// FIXME: consider renaming to multiplicativeInverse() and renaming conjugate() to inverse() or opposite().
|
||||
Scalar n2 = this->squaredNorm();
|
||||
if (n2 > Scalar(0))
|
||||
return Quaternion<Scalar>(conjugate().coeffs() / n2);
|
||||
|
||||
@@ -80,7 +80,7 @@ class UniformScaling {
|
||||
}
|
||||
|
||||
/** Concatenates a uniform scaling and a linear transformation matrix */
|
||||
// TODO returns an expression
|
||||
// TODO: return an expression instead of a dense matrix.
|
||||
template <typename Derived>
|
||||
inline typename Eigen::internal::plain_matrix_type<Derived>::type operator*(const MatrixBase<Derived>& other) const {
|
||||
return other * m_factor;
|
||||
|
||||
@@ -36,7 +36,7 @@ void make_block_householder_triangular_factor(TriangularFactorType& triFactor, c
|
||||
triFactor.row(i).tail(rt).noalias() = -hCoeffs(i) * vectors.col(i).tail(rs).adjoint() *
|
||||
vectors.bottomRightCorner(rs, rt).template triangularView<UnitLower>();
|
||||
|
||||
// FIXME use the following line with .noalias() once the triangular product can work inplace
|
||||
// FIXME: use the following line with .noalias() once triangular product supports in-place operation.
|
||||
// triFactor.row(i).tail(rt) = triFactor.row(i).tail(rt) * triFactor.bottomRightCorner(rt,rt).template
|
||||
// triangularView<Upper>();
|
||||
for (Index j = nbVecs - 1; j > i; --j) {
|
||||
@@ -71,7 +71,7 @@ void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vec
|
||||
(VectorsType::MaxColsAtCompileTime == 1 && MatrixType::MaxColsAtCompileTime != 1) ? RowMajor : ColMajor,
|
||||
VectorsType::MaxColsAtCompileTime, MatrixType::MaxColsAtCompileTime>
|
||||
tmp = V.adjoint() * mat;
|
||||
// FIXME add .noalias() once the triangular product can work inplace
|
||||
// FIXME: add .noalias() once triangular product supports in-place operation.
|
||||
if (forward)
|
||||
tmp = T.template triangularView<Upper>() * tmp;
|
||||
else
|
||||
|
||||
@@ -264,7 +264,7 @@ void IncompleteCholesky<Scalar, UpLo_, OrderingType>::factorize(const MatrixType
|
||||
else
|
||||
m_scale(j) = 1;
|
||||
|
||||
// TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
|
||||
// TODO: disable scaling when roughly uniform to speed up solve().
|
||||
|
||||
// Scale and compute the shift for the matrix
|
||||
RealScalar mindiag = NumTraits<RealScalar>::highest();
|
||||
|
||||
@@ -244,8 +244,8 @@ void IncompleteLUT<Scalar, StorageIndex>::analyzePattern(const MatrixType_& amat
|
||||
// To this end, let's symmetrize the pattern and perform AMD on it.
|
||||
SparseMatrix<Scalar, ColMajor, StorageIndex> mat1 = amat;
|
||||
SparseMatrix<Scalar, ColMajor, StorageIndex> mat2 = amat.transpose();
|
||||
// FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice.
|
||||
// on the other hand for a really non-symmetric pattern, mat2*mat1 should be preferred...
|
||||
// FIXME: for a nearly symmetric pattern, mat2+mat1 is appropriate;
|
||||
// for a highly non-symmetric pattern, mat2*mat1 should be preferred.
|
||||
SparseMatrix<Scalar, ColMajor, StorageIndex> AtA = mat2 + mat1;
|
||||
AMDOrdering<StorageIndex> ordering;
|
||||
ordering(AtA, m_P);
|
||||
|
||||
@@ -571,7 +571,7 @@ MatrixType FullPivLU<MatrixType, PermutationIndex>::reconstructedMatrix() const
|
||||
const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
|
||||
// LU
|
||||
MatrixType res(m_lu.rows(), m_lu.cols());
|
||||
// FIXME the .toDenseMatrix() should not be needed...
|
||||
// FIXME: the .toDenseMatrix() calls should not be needed.
|
||||
res = m_lu.leftCols(smalldim).template triangularView<UnitLower>().toDenseMatrix() *
|
||||
m_lu.topRows(smalldim).template triangularView<Upper>().toDenseMatrix();
|
||||
|
||||
@@ -632,10 +632,10 @@ struct kernel_retval<FullPivLU<MatrixType_, PermutationIndex_> >
|
||||
if (abs(dec().matrixLU().coeff(i, i)) > premultiplied_threshold) pivots.coeffRef(p++) = i;
|
||||
eigen_internal_assert(p == rank());
|
||||
|
||||
// we construct a temporaty trapezoid matrix m, by taking the U matrix and
|
||||
// permuting the rows and cols to bring the nonnegligible pivots to the top of
|
||||
// the main diagonal. We need that to be able to apply our triangular solvers.
|
||||
// FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
|
||||
// Construct a temporary trapezoid matrix m by taking the U matrix and permuting
|
||||
// the rows and cols to bring the nonnegligible pivots to the top of the main diagonal.
|
||||
// This is needed to apply our triangular solvers.
|
||||
// FIXME: simplify once triangularView supports rectangular matrices.
|
||||
Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, traits<MatrixType>::Options, MaxSmallDimAtCompileTime,
|
||||
MatrixType::MaxColsAtCompileTime>
|
||||
m(dec().matrixLU().block(0, 0, rank(), cols));
|
||||
|
||||
@@ -22,7 +22,7 @@ namespace internal {
|
||||
* \ingroup OrderingMethods_Module
|
||||
* \param[in] A the input non-symmetric matrix
|
||||
* \param[out] symmat the symmetric pattern A^T+A from the input matrix \a A.
|
||||
* FIXME: The values should not be considered here
|
||||
* FIXME: only the sparsity pattern should be used here; values should be ignored.
|
||||
*/
|
||||
template <typename MatrixType>
|
||||
void ordering_helper_at_plus_a(const MatrixType& A, MatrixType& symmat) {
|
||||
|
||||
@@ -185,7 +185,7 @@ class PardisoImpl : public SparseSolverBase<Derived> {
|
||||
bool symmetric = std::abs(m_type) < 10;
|
||||
m_iparm[0] = 1; // No solver default
|
||||
m_iparm[1] = 2; // use Metis for the ordering
|
||||
m_iparm[2] = 0; // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??)
|
||||
m_iparm[2] = 0; // Reserved. Set to zero. (Was number of processors / OMP_NUM_THREADS.)
|
||||
m_iparm[3] = 0; // No iterative-direct algorithm
|
||||
m_iparm[4] = 0; // No user fill-in reducing permutation
|
||||
m_iparm[5] = 0; // Write solution into x, b is left unchanged
|
||||
|
||||
@@ -575,8 +575,7 @@ template <typename RhsType, typename DstType>
|
||||
void FullPivHouseholderQR<MatrixType_, PermutationIndex_>::_solve_impl(const RhsType& rhs, DstType& dst) const {
|
||||
const Index l_rank = rank();
|
||||
|
||||
// FIXME introduce nonzeroPivots() and use it here. and more generally,
|
||||
// make the same improvements in this dec as in FullPivLU.
|
||||
// FIXME: introduce nonzeroPivots() and apply the same improvements as in FullPivLU.
|
||||
if (l_rank == 0) {
|
||||
dst.setZero();
|
||||
return;
|
||||
|
||||
@@ -384,7 +384,7 @@ void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename
|
||||
}
|
||||
}
|
||||
|
||||
// TODO: add a corresponding public API for updating a QR factorization
|
||||
// TODO: expose a public API for rank-1 QR update.
|
||||
/** \internal
|
||||
* Basically a modified copy of @c Eigen::internal::householder_qr_inplace_unblocked that
|
||||
* performs a rank-1 update of the QR matrix in compact storage. This function assumes, that
|
||||
|
||||
@@ -379,7 +379,7 @@ BDCSVD<MatrixType, Options>& BDCSVD<MatrixType, Options>::compute_impl(const Mat
|
||||
//**** step 2 - Divide & Conquer
|
||||
m_naiveU.setZero();
|
||||
m_naiveV.setZero();
|
||||
// FIXME this line involves a temporary matrix
|
||||
// FIXME: this line involves a temporary matrix.
|
||||
m_computed.topRows(diagSize()) = bid.bidiagonal().toDenseMatrix().transpose();
|
||||
m_computed.template bottomRows<1>().setZero();
|
||||
divide(0, diagSize() - 1, 0, 0, 0);
|
||||
@@ -429,7 +429,7 @@ void BDCSVD<MatrixType, Options>::copyUV(const HouseholderU& householderU, const
|
||||
m_matrixU = MatrixX::Identity(rows(), Ucols);
|
||||
m_matrixU.topLeftCorner(diagSize(), diagSize()) =
|
||||
naiveV.template cast<Scalar>().topLeftCorner(diagSize(), diagSize());
|
||||
// FIXME the following conditionals involve temporary buffers
|
||||
// FIXME: the following conditionals involve temporary buffers.
|
||||
if (m_useQrDecomp)
|
||||
m_matrixU.topLeftCorner(householderU.cols(), diagSize()).applyOnTheLeft(householderU);
|
||||
else
|
||||
@@ -440,7 +440,7 @@ void BDCSVD<MatrixType, Options>::copyUV(const HouseholderU& householderU, const
|
||||
m_matrixV = MatrixX::Identity(cols(), Vcols);
|
||||
m_matrixV.topLeftCorner(diagSize(), diagSize()) =
|
||||
naiveU.template cast<Scalar>().topLeftCorner(diagSize(), diagSize());
|
||||
// FIXME the following conditionals involve temporary buffers
|
||||
// FIXME: the following conditionals involve temporary buffers.
|
||||
if (m_useQrDecomp)
|
||||
m_matrixV.topLeftCorner(householderV.cols(), diagSize()).applyOnTheLeft(householderV);
|
||||
else
|
||||
@@ -538,7 +538,7 @@ void BDCSVD<MatrixType, Options>::divide(Index firstCol, Index lastCol, Index fi
|
||||
// We use the other algorithm which is more efficient for small
|
||||
// matrices.
|
||||
if (n < m_algoswap) {
|
||||
// FIXME this block involves temporaries
|
||||
// FIXME: this block involves temporaries.
|
||||
if (m_compV) {
|
||||
JacobiSVD<MatrixXr, ComputeFullU | ComputeFullV> baseSvd;
|
||||
computeBaseCase(baseSvd, n, firstCol, firstRowW, firstColW, shift);
|
||||
@@ -681,8 +681,8 @@ void BDCSVD<MatrixType, Options>::divide(Index firstCol, Index lastCol, Index fi
|
||||
// order except for possibly the (0,0) entry. The computed SVD is stored U, singVals and V, except
|
||||
// that if m_compV is false, then V is not computed. Singular values are sorted in decreasing order.
|
||||
//
|
||||
// TODO Opportunities for optimization: better root finding algo, better stopping criterion, better
|
||||
// handling of round-off errors, be consistent in ordering
|
||||
// TODO: opportunities for optimization: better root-finding algorithm, better stopping criterion,
|
||||
// better handling of round-off errors, and consistent ordering.
|
||||
// For instance, to solve the secular equation using FMM, see
|
||||
// http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
|
||||
template <typename MatrixType, int Options>
|
||||
|
||||
@@ -17,8 +17,8 @@
|
||||
namespace Eigen {
|
||||
|
||||
namespace internal {
|
||||
// UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API.
|
||||
// At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class.
|
||||
// UpperBidiagonalization may be replaced by a Bidiagonalization class; not part of stable API.
|
||||
// Kept for now as it is one of the few tests exercising the BandMatrix class.
|
||||
|
||||
template <typename MatrixType_>
|
||||
class UpperBidiagonalization {
|
||||
|
||||
@@ -79,8 +79,8 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r
|
||||
const Index t200 = rows / 11; // 11 == (log2(200)*1.39)
|
||||
const Index t = (rows * 100) / 139;
|
||||
|
||||
// FIXME reserve nnz non zeros
|
||||
// FIXME implement faster sorting algorithms for very small nnz
|
||||
// FIXME: reserve space for the expected number of non-zeros.
|
||||
// FIXME: implement faster sorting for very small nnz counts.
|
||||
// if the result is sparse enough => use a quick sort
|
||||
// otherwise => loop through the entire vector
|
||||
// In order to avoid to perform an expensive log2 when the
|
||||
@@ -131,7 +131,7 @@ struct conservative_sparse_sparse_product_selector<Lhs, Rhs, ResultType, ColMajo
|
||||
|
||||
// If the result is tall and thin (in the extreme case a column vector)
|
||||
// then it is faster to sort the coefficients inplace instead of transposing twice.
|
||||
// FIXME, the following heuristic is probably not very good.
|
||||
// FIXME: this heuristic has known limitations and should be improved.
|
||||
if (lhs.rows() > rhs.cols()) {
|
||||
using ColMajorMatrix = typename sparse_eval<ColMajorMatrixAux, ResultType::RowsAtCompileTime,
|
||||
ResultType::ColsAtCompileTime, ColMajorMatrixAux::Flags>::type;
|
||||
|
||||
@@ -39,7 +39,11 @@ struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, Diagonal
|
||||
: public sparse_diagonal_product_evaluator<Rhs, typename Lhs::DiagonalVectorType,
|
||||
Rhs::Flags & RowMajorBit ? SDP_AsScalarProduct : SDP_AsCwiseProduct> {
|
||||
typedef Product<Lhs, Rhs, DefaultProduct> XprType;
|
||||
enum { CoeffReadCost = HugeCost, Flags = Rhs::Flags & RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags
|
||||
enum {
|
||||
CoeffReadCost = HugeCost,
|
||||
Flags = Rhs::Flags & RowMajorBit,
|
||||
Alignment = 0
|
||||
}; // FIXME: compute proper CoeffReadCost and propagate Flags.
|
||||
|
||||
typedef sparse_diagonal_product_evaluator<Rhs, typename Lhs::DiagonalVectorType,
|
||||
Rhs::Flags & RowMajorBit ? SDP_AsScalarProduct : SDP_AsCwiseProduct>
|
||||
@@ -52,7 +56,11 @@ struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseSh
|
||||
: public sparse_diagonal_product_evaluator<Lhs, Transpose<const typename Rhs::DiagonalVectorType>,
|
||||
Lhs::Flags & RowMajorBit ? SDP_AsCwiseProduct : SDP_AsScalarProduct> {
|
||||
typedef Product<Lhs, Rhs, DefaultProduct> XprType;
|
||||
enum { CoeffReadCost = HugeCost, Flags = Lhs::Flags & RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags
|
||||
enum {
|
||||
CoeffReadCost = HugeCost,
|
||||
Flags = Lhs::Flags & RowMajorBit,
|
||||
Alignment = 0
|
||||
}; // FIXME: compute proper CoeffReadCost and propagate Flags.
|
||||
|
||||
typedef sparse_diagonal_product_evaluator<Lhs, Transpose<const typename Rhs::DiagonalVectorType>,
|
||||
Lhs::Flags & RowMajorBit ? SDP_AsCwiseProduct : SDP_AsScalarProduct>
|
||||
|
||||
@@ -1548,7 +1548,7 @@ SparseMatrix<Scalar, Options_, StorageIndex_>::operator=(const SparseMatrixBase<
|
||||
Eigen::Map<IndexVector>(dest.m_outerIndex, dest.outerSize()).setZero();
|
||||
|
||||
// pass 1
|
||||
// FIXME the above copy could be merged with that pass
|
||||
// FIXME: merge the above copy into this pass to avoid iterating twice.
|
||||
for (Index j = 0; j < otherCopy.outerSize(); ++j)
|
||||
for (typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it) ++dest.m_outerIndex[it.index()];
|
||||
|
||||
|
||||
@@ -115,7 +115,7 @@ class SparseMatrixBase : public EigenBase<Derived> {
|
||||
typedef Transpose<Derived> TransposeReturnType;
|
||||
typedef Transpose<const Derived> ConstTransposeReturnType;
|
||||
|
||||
// FIXME storage order do not match evaluator storage order
|
||||
// FIXME: storage order may not match evaluator storage order.
|
||||
typedef SparseMatrix<Scalar, Flags & RowMajorBit ? RowMajor : ColMajor, StorageIndex> PlainObject;
|
||||
|
||||
/** This is the "real scalar" type; if the \a Scalar type is already real numbers
|
||||
@@ -203,7 +203,7 @@ class SparseMatrixBase : public EigenBase<Derived> {
|
||||
return derived();
|
||||
}
|
||||
|
||||
SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */
|
||||
SparseMatrixBase() : m_isRValue(false) { /* TODO: validate traits flags. */
|
||||
}
|
||||
|
||||
template <typename OtherDerived>
|
||||
|
||||
@@ -117,7 +117,8 @@ class Ref<SparseMatrix<MatScalar, MatOptions, MatIndex>, Options, StrideType>
|
||||
#else
|
||||
template <typename SparseMatrixType, int Options>
|
||||
class Ref<SparseMatrixType, Options>
|
||||
: public SparseMapBase<Derived, WriteAccessors> // yes, that's weird to use Derived here, but that works!
|
||||
: public SparseMapBase<Derived, WriteAccessors> // Note: 'Derived' is used here intentionally; it resolves
|
||||
// correctly via CRTP.
|
||||
#endif
|
||||
{
|
||||
typedef SparseMatrix<MatScalar, MatOptions, MatIndex> PlainObjectType;
|
||||
|
||||
@@ -56,14 +56,11 @@ static void sparse_sparse_product_with_pruning_impl(const Lhs& lhs, const Rhs& r
|
||||
res.reserve(estimated_nnz_prod);
|
||||
double ratioColRes = double(estimated_nnz_prod) / (double(lhs.rows()) * double(rhs.cols()));
|
||||
for (Index j = 0; j < cols; ++j) {
|
||||
// FIXME:
|
||||
// double ratioColRes = (double(rhs.innerVector(j).nonZeros()) +
|
||||
// double(lhs.nonZeros())/double(lhs.cols()))/double(lhs.rows());
|
||||
// let's do a more accurate determination of the nnz ratio for the current column j of res
|
||||
// FIXME: compute a more accurate per-column nnz ratio for res.
|
||||
tempVector.init(ratioColRes);
|
||||
tempVector.setZero();
|
||||
for (typename evaluator<Rhs>::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt) {
|
||||
// FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
|
||||
// FIXME: rewrite as tmp += rhsIt.value() * lhs.col(rhsIt.index()).
|
||||
tempVector.restart();
|
||||
RhsScalar x = rhsIt.value();
|
||||
for (typename evaluator<Lhs>::InnerIterator lhsIt(lhsEval, rhsIt.index()); lhsIt; ++lhsIt) {
|
||||
|
||||
@@ -105,7 +105,7 @@ struct sparse_eval<T, Rows, 1, Flags> {
|
||||
typedef SparseVector<Scalar_, ColMajor, StorageIndex_> type;
|
||||
};
|
||||
|
||||
// TODO this seems almost identical to plain_matrix_type<T, Sparse>
|
||||
// TODO: consider unifying with plain_matrix_type<T, Sparse>.
|
||||
template <typename T, int Rows, int Cols, int Flags>
|
||||
struct sparse_eval {
|
||||
typedef typename traits<T>::Scalar Scalar_;
|
||||
|
||||
@@ -195,7 +195,7 @@ struct sparse_solve_triangular_sparse_selector<Lhs, Rhs, Mode, UpLo, ColMajor> {
|
||||
res.reserve(other.nonZeros());
|
||||
|
||||
for (Index col = 0; col < other.cols(); ++col) {
|
||||
// FIXME estimate number of non zeros
|
||||
// FIXME: estimate the number of non-zeros per column for better allocation.
|
||||
tempVector.init(.99 /*float(other.col(col).nonZeros())/float(other.rows())*/);
|
||||
tempVector.setZero();
|
||||
tempVector.restart();
|
||||
@@ -230,16 +230,11 @@ struct sparse_solve_triangular_sparse_selector<Lhs, Rhs, Mode, UpLo, ColMajor> {
|
||||
}
|
||||
}
|
||||
|
||||
// Index count = 0;
|
||||
// FIXME compute a reference value to filter zeros
|
||||
// FIXME: compute a reference value to filter zeros.
|
||||
for (typename AmbiVector<Scalar, StorageIndex>::Iterator it(tempVector /*,1e-12*/); it; ++it) {
|
||||
// ++ count;
|
||||
// std::cerr << "fill " << it.index() << ", " << col << "\n";
|
||||
// std::cout << it.value() << " ";
|
||||
// FIXME use insertBack
|
||||
// FIXME: use insertBack for better performance.
|
||||
res.insert(it.index(), col) = it.value();
|
||||
}
|
||||
// std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n";
|
||||
}
|
||||
res.finalize();
|
||||
other = res.markAsRValue();
|
||||
|
||||
@@ -51,7 +51,7 @@ class SparseLUTransposeView : public SparseSolverBase<SparseLUTransposeView<Conj
|
||||
eigen_assert(m_sparseLU->info() == Success && "The matrix should be factorized first");
|
||||
EIGEN_STATIC_ASSERT((Dest::Flags & RowMajorBit) == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
|
||||
|
||||
// this ugly const_cast_derived() helps to detect aliasing when applying the permutations
|
||||
// const_cast_derived() is needed to enable aliasing detection when applying the permutations.
|
||||
for (Index j = 0; j < B.cols(); ++j) {
|
||||
X.col(j) = m_sparseLU->colsPermutation() * B.const_cast_derived().col(j);
|
||||
}
|
||||
@@ -344,7 +344,7 @@ class SparseLU : public SparseSolverBase<SparseLU<MatrixType_, OrderingType_>>,
|
||||
// on return, X is overwritten by the computed solution
|
||||
X.resize(B.rows(), B.cols());
|
||||
|
||||
// this ugly const_cast_derived() helps to detect aliasing when applying the permutations
|
||||
// const_cast_derived() is needed to enable aliasing detection when applying the permutations.
|
||||
for (Index j = 0; j < B.cols(); ++j) X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
|
||||
|
||||
// Forward substitution with L
|
||||
@@ -603,7 +603,7 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) {
|
||||
* > A->ncol: number of bytes allocated when memory allocation failure occurred, plus A->ncol.
|
||||
* If lwork = -1, it is the estimated amount of space needed, plus A->ncol.
|
||||
*
|
||||
* It seems that A was the name of the matrix in the past.
|
||||
* Note: 'A' in the above description refers to the factored matrix (historical naming from SuperLU).
|
||||
*
|
||||
* \sa analyzePattern(), compute(), SparseLU(), info(), lastErrorMessage()
|
||||
*/
|
||||
@@ -616,7 +616,6 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) {
|
||||
m_isInitialized = true;
|
||||
|
||||
// Apply the column permutation computed in analyzepattern()
|
||||
// m_mat = matrix * m_perm_c.inverse();
|
||||
m_mat = matrix;
|
||||
if (m_perm_c.size()) {
|
||||
m_mat.uncompress(); // NOTE: The effect of this command is only to create the InnerNonzeros pointers.
|
||||
@@ -779,7 +778,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) {
|
||||
}
|
||||
|
||||
// Update the determinant of the row permutation matrix
|
||||
// FIXME: the following test is not correct, we should probably take iperm_c into account and pivrow is not
|
||||
// FIXME: the following test is not correct; it should account for iperm_c, and pivrow is not
|
||||
// directly the row pivot.
|
||||
if (pivrow != jj) m_detPermR = -m_detPermR;
|
||||
|
||||
|
||||
@@ -27,11 +27,7 @@ namespace internal {
|
||||
* NOTE : This class corresponds to the SCformat structure in SuperLU
|
||||
*
|
||||
*/
|
||||
/* TODO
|
||||
* InnerIterator as for sparsematrix
|
||||
* SuperInnerIterator to iterate through all supernodes
|
||||
* Function for triangular solve
|
||||
*/
|
||||
// TODO: add InnerIterator, SuperInnerIterator, and triangular solve support.
|
||||
template <typename Scalar_, typename StorageIndex_>
|
||||
class MappedSuperNodalMatrix {
|
||||
public:
|
||||
|
||||
@@ -52,7 +52,7 @@ void SparseLUImpl<Scalar, StorageIndex>::heap_relax_snode(const Index n, IndexVe
|
||||
IndexVector post;
|
||||
internal::treePostorder(StorageIndex(n), et, post); // Post order etree
|
||||
IndexVector inv_post(n + 1);
|
||||
for (StorageIndex i = 0; i < n + 1; ++i) inv_post(post(i)) = i; // inv_post = post.inverse()???
|
||||
for (StorageIndex i = 0; i < n + 1; ++i) inv_post(post(i)) = i; // Compute the inverse postorder permutation.
|
||||
|
||||
// Renumber etree in postorder
|
||||
IndexVector iwork(n);
|
||||
|
||||
@@ -136,12 +136,9 @@ void SparseLUImpl<Scalar, StorageIndex>::dfs_kernel(const StorageIndex jj, Index
|
||||
// segment is seen for the first time. (Note that
|
||||
// "repfnz(krep)" may change later.)
|
||||
// Baktrack dfs to its parent
|
||||
if (traits.update_segrep(krep, jj))
|
||||
// if (marker1(krep) < jcol )
|
||||
{
|
||||
if (traits.update_segrep(krep, jj)) {
|
||||
segrep(nseg) = krep;
|
||||
++nseg;
|
||||
// marker1(krep) = jj;
|
||||
}
|
||||
|
||||
kpar = parent(krep); // Pop recursion, mimic recursion
|
||||
|
||||
@@ -338,8 +338,7 @@ void SparseQR<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) {
|
||||
m_Q.resize(m, diagSize);
|
||||
|
||||
// Allocate space for nonzero elements: rough estimation
|
||||
m_R.reserve(2 * mat.nonZeros()); // FIXME Get a more accurate estimation through symbolic factorization with the
|
||||
// etree
|
||||
m_R.reserve(2 * mat.nonZeros()); // FIXME: get a tighter bound via symbolic factorization using the etree.
|
||||
m_Q.reserve(2 * mat.nonZeros());
|
||||
m_hcoeffs.resize(diagSize);
|
||||
m_analysisIsok = true;
|
||||
@@ -502,7 +501,7 @@ void SparseQR<MatrixType, OrderingType>::factorize(const MatrixType& mat) {
|
||||
|
||||
if (nonzeroCol < diagSize) {
|
||||
// Compute the Householder reflection that eliminate the current column
|
||||
// FIXME this step should call the Householder module.
|
||||
// FIXME: refactor to use the Householder module's reflector computation.
|
||||
Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
|
||||
|
||||
// First, the squared norm of Q((col+1):m, col)
|
||||
@@ -627,7 +626,7 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
|
||||
|
||||
const SparseQRType& m_qr;
|
||||
const Derived& m_other;
|
||||
bool m_transpose; // TODO this actually means adjoint
|
||||
bool m_transpose; // TODO: rename to m_adjoint; this flag controls adjoint application.
|
||||
};
|
||||
|
||||
template <typename SparseQRType>
|
||||
@@ -646,14 +645,14 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp
|
||||
}
|
||||
inline Index rows() const { return m_qr.rows(); }
|
||||
inline Index cols() const { return m_qr.rows(); }
|
||||
// To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment
|
||||
// To use for operations with the transpose of Q. FIXME: currently identical to adjoint(); specialize for complex.
|
||||
SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const {
|
||||
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
|
||||
}
|
||||
const SparseQRType& m_qr;
|
||||
};
|
||||
|
||||
// TODO this actually represents the adjoint of Q
|
||||
// TODO: rename to SparseQRMatrixQAdjointReturnType; this represents the adjoint of Q.
|
||||
template <typename SparseQRType>
|
||||
struct SparseQRMatrixQTransposeReturnType {
|
||||
explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
|
||||
|
||||
@@ -583,7 +583,7 @@ void SuperLU<MatrixType>::factorize(const MatrixType &a) {
|
||||
|
||||
m_extractedDataAreDirty = true;
|
||||
|
||||
// FIXME how to better check for errors ???
|
||||
// FIXME: implement more detailed error checking based on SuperLU info codes.
|
||||
m_info = info == 0 ? Success : NumericalIssue;
|
||||
m_factorizationIsOk = true;
|
||||
}
|
||||
@@ -872,7 +872,7 @@ void SuperILU<MatrixType>::factorize(const MatrixType &a) {
|
||||
&info, Scalar());
|
||||
StatFree(&m_sluStat);
|
||||
|
||||
// FIXME how to better check for errors ???
|
||||
// FIXME: implement more detailed error checking based on SuperLU info codes.
|
||||
m_info = info == 0 ? Success : NumericalIssue;
|
||||
m_factorizationIsOk = true;
|
||||
}
|
||||
|
||||
@@ -150,9 +150,9 @@ class BenchTimer {
|
||||
#define BENCH(TIMER, TRIES, REP, CODE) \
|
||||
{ \
|
||||
TIMER.reset(); \
|
||||
for (int uglyvarname1 = 0; uglyvarname1 < TRIES; ++uglyvarname1) { \
|
||||
for (int bench_tries_ = 0; bench_tries_ < TRIES; ++bench_tries_) { \
|
||||
TIMER.start(); \
|
||||
for (int uglyvarname2 = 0; uglyvarname2 < REP; ++uglyvarname2) { \
|
||||
for (int bench_reps_ = 0; bench_reps_ < REP; ++bench_reps_) { \
|
||||
CODE; \
|
||||
} \
|
||||
TIMER.stop(); \
|
||||
|
||||
@@ -257,8 +257,8 @@ void alignedboxRotatable(const BoxType& box,
|
||||
// box((-3, -2, -2), (-1, 0, 0))
|
||||
|
||||
IsometryTransform tf2 = IsometryTransform::Identity();
|
||||
// for some weird reason the following statement has to be put separate from
|
||||
// the following rotate call, otherwise precision problems arise...
|
||||
// The following statement must be separate from the rotate call below,
|
||||
// otherwise precision problems arise.
|
||||
Rotation rot = rotate(NonInteger(EIGEN_PI));
|
||||
tf2.rotate(rot);
|
||||
|
||||
|
||||
@@ -356,7 +356,7 @@ void check_indexed_view() {
|
||||
VERIFY_IS_CWISE_EQUAL(R_ref(eigen_matrix_rows, eigen_matrix_cols), R_ref(c_array_rows, c_array_cols));
|
||||
}
|
||||
|
||||
// check mat(i,j) with weird types for i and j
|
||||
// check mat(i,j) with unusual types for i and j
|
||||
{
|
||||
VERIFY_IS_APPROX(A(B.RowsAtCompileTime - 1, 1), A(3, 1));
|
||||
VERIFY_IS_APPROX(A(B.RowsAtCompileTime, 1), A(4, 1));
|
||||
@@ -420,7 +420,7 @@ void check_indexed_view() {
|
||||
// check symbolic indices
|
||||
a(last) = 1.0;
|
||||
A(last, last) = 1;
|
||||
// check weird non-const, non-lvalue scenarios
|
||||
// check unusual non-const, non-lvalue scenarios
|
||||
{
|
||||
// in these scenarios, the objects are not declared 'const', and the compiler will attempt to use the non-const
|
||||
// overloads without intervention
|
||||
|
||||
@@ -271,7 +271,7 @@ void nullary_internal_logic() {
|
||||
VERIFY((!internal::has_binary_operator<internal::linspaced_op<float> >::value));
|
||||
VERIFY((internal::functor_has_linear_access<internal::linspaced_op<float> >::ret));
|
||||
|
||||
// Regression unit test for a weird MSVC bug.
|
||||
// Regression unit test for an MSVC bug.
|
||||
// Search "nullary_wrapper_workaround_msvc" in CoreEvaluators.h for the details.
|
||||
// See also traits<Ref>::match.
|
||||
{
|
||||
|
||||
@@ -124,7 +124,7 @@ void svd_least_square(const MatrixType& m) {
|
||||
if (internal::is_same<RealScalar, double>::value || svd.rank() == m.diagonal().size()) {
|
||||
using std::sqrt;
|
||||
// This test is not stable with single precision.
|
||||
// This is probably because squaring m signicantly affects the precision.
|
||||
// This is likely because squaring m significantly affects the precision.
|
||||
if (internal::is_same<RealScalar, float>::value) ++g_test_level;
|
||||
|
||||
VERIFY_IS_APPROX(m.adjoint() * (m * x), m.adjoint() * rhs);
|
||||
|
||||
@@ -120,8 +120,7 @@ enum EulerAxis {
|
||||
template <int _AlphaAxis, int _BetaAxis, int _GammaAxis>
|
||||
class EulerSystem {
|
||||
public:
|
||||
// It's defined this way and not as enum, because I think
|
||||
// that enum is not guerantee to support negative numbers
|
||||
// Defined as static constexpr rather than enum to ensure support for negative values.
|
||||
|
||||
/** The first rotation axis */
|
||||
static constexpr int AlphaAxis = _AlphaAxis;
|
||||
|
||||
@@ -16,8 +16,8 @@ namespace Eigen {
|
||||
|
||||
namespace internal {
|
||||
|
||||
// FFTW uses non-const arguments
|
||||
// so we must use ugly const_cast calls for all the args it uses
|
||||
// FFTW uses non-const arguments,
|
||||
// so const_cast is needed for all the args it uses.
|
||||
//
|
||||
// This should be safe as long as
|
||||
// 1. we use FFTW_ESTIMATE for all our planning
|
||||
|
||||
@@ -21,8 +21,7 @@ namespace internal {
|
||||
// post: sqrtT.block(i,i,2,2) is square root of T.block(i,i,2,2)
|
||||
template <typename MatrixType, typename ResultType>
|
||||
void matrix_sqrt_quasi_triangular_2x2_diagonal_block(const MatrixType& T, Index i, ResultType& sqrtT) {
|
||||
// TODO: This case (2-by-2 blocks with complex conjugate eigenvalues) is probably hidden somewhere
|
||||
// in EigenSolver. If we expose it, we could call it directly from here.
|
||||
// TODO: this 2x2 complex-conjugate eigenvalue case could reuse logic from EigenSolver if exposed.
|
||||
typedef typename traits<MatrixType>::Scalar Scalar;
|
||||
Matrix<Scalar, 2, 2> block = T.template block<2, 2>(i, i);
|
||||
EigenSolver<Matrix<Scalar, 2, 2> > es(block);
|
||||
|
||||
@@ -67,8 +67,7 @@ void lmpar(Matrix<Scalar, Dynamic, Dynamic> &r, const VectorXi &ipvt, const Matr
|
||||
l = ipvt[j];
|
||||
wa1[j] = diag[l] * (wa2[l] / dxnorm);
|
||||
}
|
||||
// it's actually a triangularView.solveInplace(), though in a weird
|
||||
// way:
|
||||
// Triangular solve (forward substitution):
|
||||
for (j = 0; j < n; ++j) {
|
||||
Scalar sum = 0.;
|
||||
for (i = 0; i < j; ++i) sum += r(i, j) * wa1[i];
|
||||
|
||||
Reference in New Issue
Block a user