Ublas :

Date view	Thread view	Subject view	Author view

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

Next message: Gunter Winkler: "Re: [ublas] BOOST_UBLAS_USE_ITERATING vs lu_factorize: broken?"
Previous message: Rutger ter Borg: "Re: [ublas] bindings"
Next in thread: Gunter Winkler: "Re: [ublas] BOOST_UBLAS_USE_ITERATING vs lu_factorize: broken?"
Reply: Gunter Winkler: "Re: [ublas] BOOST_UBLAS_USE_ITERATING vs lu_factorize: broken?"

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

(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}

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();
  }

Next message: Gunter Winkler: "Re: [ublas] BOOST_UBLAS_USE_ITERATING vs lu_factorize: broken?"
Previous message: Rutger ter Borg: "Re: [ublas] bindings"
Next in thread: Gunter Winkler: "Re: [ublas] BOOST_UBLAS_USE_ITERATING vs lu_factorize: broken?"
Reply: Gunter Winkler: "Re: [ublas] BOOST_UBLAS_USE_ITERATING vs lu_factorize: broken?"

Date view	Thread view	Subject view	Author view