10#ifndef EIGEN_CONDITIONESTIMATOR_H
11#define EIGEN_CONDITIONESTIMATOR_H
17template <
typename Vector,
typename RealVector,
bool IsComplex>
18struct rcond_compute_sign {
20 const RealVector v_abs = v.cwiseAbs();
21 return (v_abs.array() ==
static_cast<typename Vector::RealScalar
>(0))
22 .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs));
27template <
typename Vector>
30 return (v.array() <
static_cast<typename Vector::RealScalar
>(0))
31 .select(-Vector::Ones(v.size()), Vector::Ones(v.size()));
55template <
typename Decomposition>
56typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(
const Decomposition& dec)
58 typedef typename Decomposition::MatrixType MatrixType;
59 typedef typename Decomposition::Scalar Scalar;
60 typedef typename Decomposition::RealScalar RealScalar;
61 typedef typename internal::plain_col_type<MatrixType>::type
Vector;
62 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector;
63 const bool is_complex = (NumTraits<Scalar>::IsComplex != 0);
65 eigen_assert(dec.rows() == dec.cols());
66 const Index n = dec.rows();
71#ifdef __INTEL_COMPILER
73 #pragma warning ( disable : 2259 )
75 Vector v = dec.solve(Vector::Ones(n) / Scalar(n));
76#ifdef __INTEL_COMPILER
84 RealScalar lower_bound = v.template lpNorm<1>();
91 RealScalar old_lower_bound = lower_bound;
94 Index v_max_abs_index = -1;
95 Index old_v_max_abs_index = v_max_abs_index;
96 for (
int k = 0; k < 4; ++k)
98 sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v);
99 if (k > 0 && !is_complex && sign_vector == old_sign_vector) {
104 v = dec.adjoint().solve(sign_vector);
105 v.real().cwiseAbs().maxCoeff(&v_max_abs_index);
106 if (v_max_abs_index == old_v_max_abs_index) {
111 v = dec.solve(Vector::Unit(n, v_max_abs_index));
112 lower_bound = v.template lpNorm<1>();
113 if (lower_bound <= old_lower_bound) {
118 old_sign_vector = sign_vector;
120 old_v_max_abs_index = v_max_abs_index;
121 old_lower_bound = lower_bound;
133 Scalar alternating_sign(RealScalar(1));
134 for (
Index i = 0; i < n; ++i) {
136 v[i] = alternating_sign *
static_cast<RealScalar
>(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1))));
137 alternating_sign = -alternating_sign;
140 const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n));
141 return numext::maxi(lower_bound, alternate_lower_bound);
157template <
typename Decomposition>
158typename Decomposition::RealScalar
159rcond_estimate_helper(
typename Decomposition::RealScalar matrix_norm,
const Decomposition& dec)
161 typedef typename Decomposition::RealScalar RealScalar;
162 eigen_assert(dec.rows() == dec.cols());
163 if (dec.rows() == 0)
return NumTraits<RealScalar>::infinity();
164 if (matrix_norm == RealScalar(0))
return RealScalar(0);
165 if (dec.rows() == 1)
return RealScalar(1);
166 const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);
167 return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0)
168 : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
Matrix< Type, Size, 1 > Vector
[c++11]
Definition: Matrix.h:551
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74