Re: [ublas] Does uBLAS work with dd_real and qd_real ?

Subject: Re: [ublas] Does uBLAS work with dd_real and qd_real ?
From: Paul Leopardi (paul.leopardi_at_[hidden])
Date: 2010-03-25 10:33:16

On Friday 26 March 2010 00:27:09 Nasos Iliopoulos wrote:
> Paul,
> the problem seems to be in the implementation of the dd_real in the qd
> library. I noticed that the default constructor does not initialize the
> number class, hence constructing it in an ill-formed state. If you
> recompile the qd library changing dd_real() { }
> to:
>   dd_real() { x[0]=0.0; x[1]=0.0;}

Hi Nasos,
Yes, I noticed this too. To my mind this is strange behaviour on the part of 
uBLAS. Specifically, if the type of m is matrix<scalar_t> then m.clear() does 
not set the values of m to scalar_t(0) but rather to scalar_t(), i.e. it calls 
the default constructor rather than the constructor from 0 (if any).

What does the C++ standard say about the result of calling the default 
constructor for float or double [ i.e. float(), double() ]? It says that since 
float and double are POD (Plain Old Data) classes, the result is zero-
initialization. But for non-POD classes T, calling T() just results in default 
initialization.

http://www.fnal.gov/docs/working-groups/fpcltf/Pkg/ISOcxx/doc/POD.html
http://www.boost.org/doc/libs/1_33_1/libs/utility/value_init.htm

Should dd_real be changed so that dd_real() yields 0, thus forcing default 
initialization to zero, or should clear() be changed so that it calls 
value_type(0)? In either case, since I am writing a third library which is 
intended to be used with both, my next course of action would be to submit a 
patch to have the behaviour changed. 

My vote would be to change uBLAS so that it tests a  zero_initialization trait 
at compile time before calling value_type() [implicit initialization to zero], 
otherwise tests a zero_constructible trait before calling value_type(0) 
[explicit initialization to zero], and only otherwise calls value_type() 
[default initialization]. In other words, clear() should try to clear to zero, 
implicitly if possible and explictly where necessary and possible. My GluCat 
library could then define the missing traits for dd_real and qd_real to force 
uBLAS to do the right thing. 

The alternative would be to submit a patch to qd to ask that dd_real and 
qd_real be changed to do default-initialization to zero.
Best, Paul

PS:

I created the file clear_test.cpp with the following contents,

#include <qd/qd_real.h>

#include <boost/numeric/ublas/fwd.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

template < typename scalar_t >
void
test_clear(int dim)
{
  namespace ublas = boost::numeric::ublas;
  typedef ublas::matrix< scalar_t, ublas::row_major > lhs_t;

  std::cout << "default constructor == " << scalar_t() << std::endl;

  lhs_t lhs(dim, dim);
  lhs.clear();

  std::cout << "lhs:" << lhs << std::endl;
}

int main(int argc, char ** argv)
{ 
  std::cout << "float:" << std::endl;
  test_clear< float >(1);
  std::cout << std::endl;

  std::cout << "double:" << std::endl;
  test_clear< double >(1);
  std::cout << std::endl;

  unsigned int old_fpu_control;
  fpu_fix_start(&old_fpu_control);

  std::cout << "dd_real:" << std::endl;
  test_clear< dd_real >(1);
  std::cout << std::endl;

  std::cout << "qd_real:" << std::endl;
  test_clear< qd_real >(1);
  std::cout << std::endl;

  return 0;
}

It produces the following output:

 ./clear_test.exe
float:
default constructor == 0
lhs:[1,1]((0))

double:
default constructor == 0
lhs:[1,1]((0))

dd_real:
default constructor == ERROR (dd_real::to_digits): can't compute exponent.
.e+00
lhs:ERROR (dd_real::to_digits): can't compute exponent.
[1,1]((h. <strange characters>

qd_real:
default constructor == 2.083705e-309
terminate called after throwing an instance of 'std::logic_error'
  what():  (qd_real::to_digits): can't compute exponent.
lhs:Aborted

After changing all instances of clear() in matrix.h to use value_type(0) 
rather than value_type(), and changing the "default constructor" output line 
in clear_test.cpp to:

  std::cout << "zero constructor == " << scalar_t(0) << std::endl;

the output becomes:

 ./clear_test.exe
float:
zero constructor == 0
lhs:[1,1]((0))

double:
zero constructor == 0
lhs:[1,1]((0))

dd_real:
zero constructor == 0.000000e+00
lhs:[1,1]((0.000000e+00))

qd_real:
zero constructor == 0.000000e+00
lhs:[1,1]((0.000000e+00))