[ublas] Floating point exceptions and unexpected call to sqrt with sparse integer vectors

Subject: [ublas] Floating point exceptions and unexpected call to sqrt with sparse integer vectors
From: Mark Johnson (mj1_at_[hidden])
Date: 2009-08-10 15:48:12

Thanks Jesse for your answer to my previous question.

A lot of my work involves statistics of sparse vectors of positive 
integers, so it seems ublas sparse vectors and matrices are going to be 
very useful for me.

However, I find that sparse vectors of integers produce strange floating 
point exceptions (they work great with floats and doubles).  When I use 
gdb to trace these exceptions I find unexplained calls to sqrt(), even 
though my program doesn't involve sqrt() at all.  I'd appreciate any advice!

Here's a simple program that demonstrates the problem; I use 
feenableexcept() to trap floating point exceptions.

    #include <boost/numeric/ublas/io.hpp>
    #include <boost/numeric/ublas/vector.hpp>
    #include <boost/numeric/ublas/vector_expression.hpp>
    #include <boost/numeric/ublas/vector_proxy.hpp>
    #include <boost/numeric/ublas/vector_sparse.hpp>

    #include <iostream>

    // _GNU_SOURCE must be defined and fenv.h included to
    // trap floating-point exceptions (fpes).
    //
    #define _GNU_SOURCE 1
    #include <fenv.h>

    namespace ublas = boost::numeric::ublas;

    int main(int argc, char** argv) {

       // trap floating point exceptions
       feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);

       ublas::mapped_vector<int> v = ublas::unit_vector<int>(2, 1);
       std::cout << "v = " << v << std::endl;
    }

When I compile and run this program I get a floating point exception, 
which is curious as there's no floating point math in the program at all.

    [mj_at_lugha tmp]$ g++ -g -O0 ublas-int.cc -o ublas-int
    [mj_at_lugha tmp]$ ./ublas-int
    Floating point exception

Gdb shows that the exception is generated in a call to sqrt(-2147483648)

    [mj_at_lugha tmp]$ gdb ublas-int
    ...(gdb) run
    Starting program: /home/mj/tmp/ublas-int

    Program received signal SIGFPE, Arithmetic exception.
    0x0000000000405468 in std::sqrt<int> (__x=-2147483648) at
    /usr/lib/gcc/x86_64-redhat-linux/4.4.0/../../../../include/c++/4.4.0/cmath:439
    439        { return __builtin_sqrt(__x); }
    Missing separate debuginfos, use: debuginfo-install
    glibc-2.10.1-2.x86_64 libgcc-4.4.0-4.x86_64 libstdc++-4.4.0-4.x86_64
    (gdb) bt
    #0  0x0000000000405468 in std::sqrt<int> (__x=-2147483648) at
    /usr/lib/gcc/x86_64-redhat-linux/4.4.0/../../../../include/c++/4.4.0/cmath:439
    #1  0x0000000000401ac0 in (anonymous
    namespace)::boost_numeric_ublas_sqrt<int> (t=@0x7fffffffdd98) at
    /usr/include/boost/numeric/ublas/traits.hpp:32
    #2  0x0000000000404759 in
    boost::numeric::ublas::scalar_traits<int>::type_sqrt
    (t=@0x7fffffffdd98) at /usr/include/boost/numeric/ublas/traits.hpp:107
    #3  0x0000000000403fbe in
    boost::numeric::ublas::detail::expression_type_check<boost::numeric::ublas::mapped_vector<int,
    boost::numeric::ublas::map_std<unsigned long, int,
    std::allocator<std::pair<unsigned long const, int> > > >,
    boost::numeric::ublas::vector<int,
    boost::numeric::ublas::unbounded_array<int, std::allocator<int> > >
     > (e1=@0x7fffffffdee0, e2=@0x7fffffffde00) at
    /usr/include/boost/numeric/ublas/detail/vector_assign.hpp:42
    #4  0x0000000000403466 in
    boost::numeric::ublas::vector_assign<boost::numeric::ublas::scalar_assign,
    boost::numeric::ublas::mapped_vector<int,
    boost::numeric::ublas::map_std<unsigned long, int,
    std::allocator<std::pair<unsigned long const, int> > > >,
    boost::numeric::ublas::unit_vector<int, std::allocator<int> > >
    (v=@0x7fffffffdee0, e=@0x7fffffffdf20) at
    /usr/include/boost/numeric/ublas/detail/vector_assign.hpp:369
    #5  0x0000000000402e5c in
    boost::numeric::ublas::vector_assign<boost::numeric::ublas::scalar_assign,
    boost::numeric::ublas::mapped_vector<int,
    boost::numeric::ublas::map_std<unsigned long, int,
    std::allocator<std::pair<unsigned long const, int> > > >,
    boost::numeric::ublas::unit_vector<int, std::allocator<int> > >
    (v=@0x7fffffffdee0, e=@0x7fffffffdf20) at
    /usr/include/boost/numeric/ublas/detail/vector_assign.hpp:448
    #6  0x0000000000402a19 in boost::numeric::ublas::mapped_vector<int,
    boost::numeric::ublas::map_std<unsigned long, int,
    std::allocator<std::pair<unsigned long const, int> > >
     >::mapped_vector<boost::numeric::ublas::unit_vector<int,
    std::allocator<int> > > (this=0x7fffffffdee0, ae=@0x7fffffffdf20,
    non_zeros=0) at /usr/include/boost/numeric/ublas/vector_sparse.hpp:312
    #7  0x0000000000401a16 in main (argc=1, argv=0x7fffffffe028) at
    ublas-int.cc:22
    (gdb)

Any ideas as to what's happening?  I must admit I don't understand why 
the program is calling sqrt()!

Thanks in advance,

Mark