# Ublas :

Subject: [ublas] BOOST_UBLAS_USE_ITERATING vs lu_factorize: broken?
From: Paul C. Leopardi (paul.leopardi_at_[hidden])
Date: 2009-08-01 01:27:33

Has anyone checked lu_factorize with -DBOOST_UBLAS_USE_ITERATING?
It's giving me stranger results. I'm using Boost 1.37.0 with g++ 4.4.0.
Best, Paul

The test function from GluCat 0.4.1:

int test13()
{
using namespace glucat;
using namespace std;
cout << "Programming example 13 : Multiplication and division" << endl;
typedef matrix_multi<long double> cm;
typedef framed_multi<long double> cf;
const cf a("{-3}+{-2}+{-1}");
const cf b("{-2}+1.e-8{1}+{2}+1.e8{3}");

const cf::index_set_t sub = a.frame() | b.frame();
const cm A( a, sub );
const cm B( b, sub );
cout << "a: " << endl << a << endl;
cout << "A = cm(a): " << endl << A << endl;
cout << "b: " << endl << b << endl;
cout << "B = cm(b): " << endl << B << endl;

// Multiplication
const cf c = a*b;
cout << "a*b: " << endl << c << endl;
const cm C = A * B;
cout << "A*B: " << endl << C << endl;

// Division
cout << "B/B: " << endl << (B / B) << endl;
cout << "(B/B)*B: " << endl << (B / B) * B << endl;
[snip...]

Output without -DBOOST_UBLAS_USE_ITERATING:
Programming example 13 : Multiplication and division
a:
{-3}+{-2}+{-1}
A = cm(a):
{-3}+{-2}+{-1}
b:
{-2}+1e-08{1}+{2}+1e+08{3}
B = cm(b):
{-2}+1e-08{1}+{2}+1e+08{3}
a*b:
-1+{-3,-2}-{-2,-1}+1e-08{-3,1}+1e-08{-2,1}+1e-08
{-1,1}+{-3,2}+{-2,2}+{-1,2}+1e+08{-3,3}+1e+08{-2,3}+1e+08{-1,3}
A*B:
-1+{-3,-2}-{-2,-1}+9.99717e-09{-3,1}+1e-08{-2,1}+9.99717e-09
{-1,1}+{-3,2}+{-2,2}+{-1,2}+1e+08{-3,3}+1e+08{-2,3}+1e+08{-1,3}

AT==[8,8]((0,1e-08,2,0,1e+08,0,0,0),(1e-08,0,0,2,0,1e+08,0,0),
(0,0,0,-1e-08,0,0,1e+08,0),(0,0,-1e-08,0,0,0,0,1e+08),
(1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0))
BT==[8,8]((0,1e-08,2,0,1e+08,0,0,0),(1e-08,0,0,2,0,1e+08,0,0),
(0,0,0,-1e-08,0,0,1e+08,0),(0,0,-1e-08,0,0,0,0,1e+08),
(1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0))
LU==[8,8]((1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0),
(0,1e-16,2e-08,0,1e+08,0,0,0),(1e-16,0,0,2e-08,0,1e+08,0,0),
(0,0,0,-1e-16,0,0,1e+08,0),(0,0,-1e-16,0,0,0,0,1e+08))
XT==[8,8]((1,0,0,0,0,0,0,0),(0,1,0,0,0,0,0,0),(0,0,1,0,0,0,0,0),
(0,0,0,1,0,0,0,0),(0,0,0,0,1,0,0,0),(0,0,0,0,0,1,0,0),(0,0,0,0,0,0,1,0),
(0,0,0,0,0,0,0,1))
XT==[8,8]((1,0,0,0,0,0,0,0),(0,1,0,0,0,0,0,0),(0,0,1,0,0,0,0,0),
(0,0,0,1,0,0,0,0),(0,0,0,0,1,0,0,0),(0,0,0,0,0,1,0,0),(0,0,0,0,0,0,1,0),
(0,0,0,0,0,0,0,1))

B/B:
1

AT==[8,8]((0,1e-08,2,0,1e+08,0,0,0),(1e-08,0,0,2,0,1e+08,0,0),
(0,0,0,-1e-08,0,0,1e+08,0),(0,0,-1e-08,0,0,0,0,1e+08),
(1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0))
BT==[8,8]((0,1e-08,2,0,1e+08,0,0,0),(1e-08,0,0,2,0,1e+08,0,0),
(0,0,0,-1e-08,0,0,1e+08,0),(0,0,-1e-08,0,0,0,0,1e+08),
(1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0))
LU==[8,8]((1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0),
(0,1e-16,2e-08,0,1e+08,0,0,0),(1e-16,0,0,2e-08,0,1e+08,0,0),
(0,0,0,-1e-16,0,0,1e+08,0),(0,0,-1e-16,0,0,0,0,1e+08))
XT==[8,8]((1,0,0,0,0,0,0,0),(0,1,0,0,0,0,0,0),(0,0,1,0,0,0,0,0),
(0,0,0,1,0,0,0,0),(0,0,0,0,1,0,0,0),(0,0,0,0,0,1,0,0),(0,0,0,0,0,0,1,0),
(0,0,0,0,0,0,0,1))
XT==[8,8]((1,0,0,0,0,0,0,0),(0,1,0,0,0,0,0,0),(0,0,1,0,0,0,0,0),
(0,0,0,1,0,0,0,0),(0,0,0,0,1,0,0,0),(0,0,0,0,0,1,0,0),(0,0,0,0,0,0,1,0),
(0,0,0,0,0,0,0,1))

(B/B)*B:
{-2}+1e-08{1}+{2}+1e+08{3}
[snip...]

Output with -DBOOST_UBLAS_USE_ITERATING:
Programming example 13 : Multiplication and division
a:
{-3}+{-2}+{-1}
A = cm(a):
{-3}+{-2}+{-1}
b:
{-2}+1e-08{1}+{2}+1e+08{3}
B = cm(b):
{-2}+1e-08{1}+{2}+1e+08{3}
a*b:
-1+{-3,-2}-{-2,-1}+1e-08{-3,1}+1e-08{-2,1}+1e-08
{-1,1}+{-3,2}+{-2,2}+{-1,2}+1e+08{-3,3}+1e+08{-2,3}+1e+08{-1,3}
A*B:
-1+{-3,-2}-{-2,-1}+9.99717e-09{-3,1}+1e-08{-2,1}+9.99717e-09
{-1,1}+{-3,2}+{-2,2}+{-1,2}+1e+08{-3,3}+1e+08{-2,3}+1e+08{-1,3}

AT==[8,8]((0,1e-08,2,0,1e+08,0,0,0),(1e-08,0,0,2,0,1e+08,0,0),
(0,0,0,-1e-08,0,0,1e+08,0),(0,0,-1e-08,0,0,0,0,1e+08),
(1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0))
B/B:
nan

AT==[8,8]((0,1e-08,2,0,1e+08,0,0,0),(1e-08,0,0,2,0,1e+08,0,0),
(0,0,0,-1e-08,0,0,1e+08,0),(0,0,-1e-08,0,0,0,0,1e+08),
(1e+08,0,0,0,0,-1e-08,-2,0),(0,1e+08,0,0,-1e-08,0,0,-2),
(0,0,1e+08,0,0,0,0,1e-08),(0,0,0,1e+08,0,0,1e-08,0))
(B/B)*B:
0
[snip...]

The division function from GluCat 0.4.1, including extra std::cout statements
for debugging:

/// Quotient of multivector and scalar
template< typename Scalar_T, const index_t LO, const index_t HI >
inline
matrix_multi<Scalar_T,LO,HI>&
matrix_multi<Scalar_T,LO,HI>::
operator/= (const Scalar_T& scr)
{ return *this *= Scalar_T(1)/scr; }

/// Geometric quotient
template< typename Scalar_T, const index_t LO, const index_t HI >
const matrix_multi<Scalar_T,LO,HI>
operator/ (const matrix_multi<Scalar_T,LO,HI>& lhs,
const matrix_multi<Scalar_T,LO,HI>& rhs)
{
#ifndef _GLUCAT_USE_DENSE_MATRICES
if (lhs.isnan() || rhs.isnan())
return numeric_traits<Scalar_T>::NaN();
#endif

if (rhs == Scalar_T(0))
return numeric_traits<Scalar_T>::NaN();

typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
typedef typename multivector_t::index_set_t index_set_t;
typedef typename multivector_t::framed_multi_t framed_multi_t;

// Operate only within a common frame
const index_set_t our_frame = lhs.m_frame | rhs.m_frame;
const multivector_t& lhs_ref = (lhs.m_frame == our_frame)
? lhs
: multivector_t(framed_multi_t(lhs),
our_frame, true);
const multivector_t& rhs_ref = (rhs.m_frame == our_frame)
? rhs
: multivector_t(framed_multi_t(rhs),
our_frame, true);

// Solve result == lhs_ref/rhs_ref <=> result*rhs_ref == lhs_ref
// We now solve X == B/A
// (where X == result, B == lhs_ref.m_matrix and A == rhs_ref.m_matrix)
// X == B/A <=> X*A == B <=> AT*XT == BT
// So, we solve AT*XT == BT

typedef typename multivector_t::matrix_t matrix_t;
typedef typename matrix_t::size_type matrix_index_t;

const matrix_t& AT = ublas::trans(rhs_ref.m_matrix);
matrix_t LU = AT;
std::cout << std::endl;
std::cout << "AT==" << AT << std::endl;

typedef ublas::permutation_matrix<matrix_index_t> permutation_t;

permutation_t pvector(AT.size1());
if (! ublas::lu_factorize(LU, pvector))
{
const matrix_t& BT = ublas::trans(lhs_ref.m_matrix);
std::cout << "BT==" << BT << std::endl;
matrix_t XT = BT;
std::cout << "LU==" << LU << std::endl;
ublas::lu_substitute(LU, pvector, XT);
#ifndef _GLUCAT_USE_DENSE_MATRICES
if (matrix::isnan(XT))
return numeric_traits<Scalar_T>::NaN();
#endif
std::cout << "XT==" << XT << std::endl;

// Iterative refinement.
// Reference: Nicholas J. Higham,
// "Accuracy and Stability of Numerical Algorithms",
// SIAM, 1996, ISBN 0-89871-355-2, Chapter 11
if (Tune_P::div_max_steps > 0)
{
// matrix_t R = ublas::prod(AT, XT) - BT;
matrix_t R = -BT;
ublas::axpy_prod(AT, XT, R, false);
#ifndef _GLUCAT_USE_DENSE_MATRICES
if (matrix::isnan(R))
return numeric_traits<Scalar_T>::NaN();
#endif

Scalar_T nr = ublas::norm_inf(R);
if ( nr != Scalar_T(0)
&& !numeric_traits<Scalar_T>::isNaN_or_isInf(nr) )
{
matrix_t XTnew = XT;
Scalar_T nrold = nr + Scalar_T(1);
for (int
step = 0;
step != Tune_P::div_max_steps &&
nr < nrold &&
nr != Scalar_T(0) &&
nr == nr;
++step)
{
nrold = nr;
if (step != 0)
XT = XTnew;
matrix_t& D = R;
ublas::lu_substitute(LU, pvector, D);
XTnew -= D;
// noalias(R) = ublas::prod(AT, XTnew) - BT;
R = -BT;
ublas::axpy_prod(AT, XTnew, R, false);
nr = ublas::norm_inf(R);
}
}
}
std::cout << "XT==" << XT << std::endl;
std::cout << std::endl;
return multivector_t(ublas::trans(XT), our_frame);
}
else
// AT is singular. Return NaN
return numeric_traits<Scalar_T>::NaN();
}