mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Rewrite balancer to avoid overflows.

Browse Source 
The previous balancer overflowed for large row/column norms.
Modified to prevent that.

Fixes #2273.
This commit is contained in:
Antonio Sanchez 2021-06-18 14:24:11 -07:00 committed by Rasmus Munk Larsen
parent 35a367d557
commit e9ab4278b7
1 changed files with 28 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -20,12 +20,6 @@ namespace internal {
#ifndef EIGEN_PARSED_BY_DOXYGEN
template <typename T>
T radix(){ return 2; }
template <typename T>
T radix2(){ return radix<T>()*radix<T>(); }
template<int Size>
struct decrement_if_fixed_size
{
@ -141,7 +135,10 @@ inline
bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, RealScalar& colB, RealScalar& rowB )
{
if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; }
if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm
|| !(numext::isfinite)(colNorm) || !(numext::isfinite)(rowNorm)){
return true;
}
else
{
//To find the balancing coefficients, if the radix is 2,
@ -149,33 +146,41 @@ bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
// \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$
// then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$
// and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
rowB = rowNorm / radix<RealScalar>();
const RealScalar radix = RealScalar(2);
const RealScalar radix2 = RealScalar(4);
rowB = rowNorm / radix;
colB = RealScalar(1);
const RealScalar s = colNorm + rowNorm;
while (colNorm < rowB)
// Find sigma s.t. rowNorm / 2 <= 2^(2*sigma) * colNorm
RealScalar scout = colNorm;
while (scout < rowB)
{
colB *= radix<RealScalar>();
colNorm *= radix2<RealScalar>();
colB *= radix;
scout *= radix2;
}
// We now have an upper-bound for sigma, try to lower it.
// Find sigma s.t. 2^(2*sigma) * colNorm / 2 < rowNorm
scout = colNorm * (colB / radix) * colB; // Avoid overflow.
while (scout >= rowNorm)
{
colB /= radix;
scout /= radix2;
}
rowB = rowNorm * radix<RealScalar>();
while (colNorm >= rowB)
{
colB /= radix<RealScalar>();
colNorm /= radix2<RealScalar>();
}
//This line is used to avoid insubstantial balancing
if ((rowNorm + colNorm) < RealScalar(0.95) * s * colB)
// This line is used to avoid insubstantial balancing.
if ((rowNorm + radix * scout) < RealScalar(0.95) * s * colB)
{
isBalanced = false;
rowB = RealScalar(1) / colB;
return false;
}
else{
return true; }
else
{
return true;
}
}
}