Fix bug in alternate lower bound calculation due to missing parentheses.

Make a few expressions more concise.
This commit is contained in:
Rasmus Munk Larsen 2016-04-05 16:40:48 -07:00
parent d7eeee0c1d
commit 4d07064a3d

View File

@ -134,8 +134,7 @@ typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
if (n == 0) { if (n == 0) {
return 0; return 0;
} }
Vector v = Vector::Ones(n) / n; Vector v = dec.solve(Vector::Ones(n) / n);
v = dec.solve(v);
// lower_bound is a lower bound on // lower_bound is a lower bound on
// ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1
@ -155,12 +154,10 @@ typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
int old_v_max_abs_index = v_max_abs_index; int old_v_max_abs_index = v_max_abs_index;
for (int k = 0; k < 4; ++k) { for (int k = 0; k < 4; ++k) {
sign_vector = internal::SignOrUnity<Vector, RealVector, is_complex>::run(v); sign_vector = internal::SignOrUnity<Vector, RealVector, is_complex>::run(v);
if (k > 0 && !is_complex) { if (k > 0 && !is_complex && sign_vector == old_sign_vector) {
if (sign_vector == old_sign_vector) {
// Break if the solution stagnated. // Break if the solution stagnated.
break; break;
} }
}
// v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )| // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )|
v = dec.adjoint().solve(sign_vector); v = dec.adjoint().solve(sign_vector);
v.real().cwiseAbs().maxCoeff(&v_max_abs_index); v.real().cwiseAbs().maxCoeff(&v_max_abs_index);
@ -193,8 +190,9 @@ typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
// Hager's algorithm to vastly underestimate ||matrix||_1. // Hager's algorithm to vastly underestimate ||matrix||_1.
Scalar alternating_sign = 1; Scalar alternating_sign = 1;
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
v[i] = alternating_sign * static_cast<RealScalar>(1) + v[i] = alternating_sign *
(static_cast<RealScalar>(i) / (static_cast<RealScalar>(n - 1))); (static_cast<RealScalar>(1) +
(static_cast<RealScalar>(i) / (static_cast<RealScalar>(n - 1))));
alternating_sign = -alternating_sign; alternating_sign = -alternating_sign;
} }
v = dec.solve(v); v = dec.solve(v);