[Trilinos-Users] Tpetra and boost::interval

[Sensei]

Dear all,

I'm playing with Tpetra to see if I could move away from Epetra. So far 
so good, I can use vectors of double, but then I thought I'd like to see 
how it works with custom classes.

The first thing that came to my mind, is using the Boost Interval 
library, to avoid making errors in a class of mine. The source code is 
at the end of this email.

I am getting compile errors regarding boost and teuchos (and tpetra), as 
follows:



/usr/local/include/trilinos/Teuchos_ScalarTraitsDecl.hpp:59:44: No 
member named 'this_type_is_missing_a_specialization' in 
'boost::numeric::interval<double, 
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, 
boost::numeric::interval_lib::checking_strict<double> > >'

/usr/local/include/trilinos/Tpetra_MultiVector_def.hpp:2520:25: Invalid 
operands to binary expression ('Teuchos::FancyOStream' (aka 
'basic_FancyOStream<char>') and 'const boost::numeric::interval<double, 
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, 
boost::numeric::interval_lib::checking_strict<double> > >')

/usr/local/include/trilinos/Teuchos_SerializationTraits.hpp:67:43: No 
member named 'this_type_is_missing_a_specialization' in 
'boost::numeric::interval<double, 
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, 
boost::numeric::interval_lib::checking_strict<double> > >'

/usr/local/include/trilinos/Tpetra_Vector_def.hpp:259:23: Invalid 
operands to binary expression ('basic_ostream<char, 
std::__1::char_traits<char> >' and 'const 
boost::numeric::interval<double, 
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, 
boost::numeric::interval_lib::checking_strict<double> > >')



Of course, if I set INTERVAL to 0, then compiling with a simple double 
scalar type works flawlessly... but I'd like to see if I could have a 
custom scalar type.

Is it possible to make vectors/matrices (and then solve linear problems) 
with a custom scalar?

I might be missing something trivial here, so please correct my naive 
errors, if any!


Thank You!



#define INTERVAL 1

#include <iostream>

#include <mpi.h>

#include <boost/numeric/interval.hpp>
#include <Tpetra_DefaultPlatform.hpp>
#include <Tpetra_Vector.hpp>
#include <Tpetra_Version.hpp>
#include <Teuchos_GlobalMPISession.hpp>
#include <Teuchos_oblackholestream.hpp>

typedef boost::numeric::interval<double> scalar;

//std::ostream& operator<< (std::ostream &s, const scalar &v)
//{
//    s << "[" << v.lower() << "," << v.upper() << "]";
//    return s;
//}

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

     Teuchos::oblackholestream blackHole;

     Teuchos::GlobalMPISession session(&argc, &argv, &blackHole);

     auto comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();

     const int rank = comm->getRank ();

     std::ostream& log = (rank == 0) ? std::cout : blackHole;

     auto map = Teuchos::rcp(new Tpetra::Map<int, int>(20, 0, comm));

     std::cout << " rank " << rank << " : LOCAL " << 
map->getMinLocalIndex() << " to " << map->getMaxLocalIndex()
               << " GLOBAL " << map->getMinGlobalIndex() << " to " << 
map->getMaxGlobalIndex() << std::endl;

#if INTERVAL
     auto v = Teuchos::rcp(new Tpetra::Vector<scalar>(map));
#else
     auto v = Teuchos::rcp(new Tpetra::Vector<double>(map));
#endif

     v->randomize();
     v->print(std::cout);

     for (int i = map->getMinLocalIndex(); i <= map->getMaxLocalIndex(); 
i++)
#if INTERVAL
         v->replaceLocalValue(i, scalar(-1.0 * (rank + 1), i));
#else
         v->replaceLocalValue(i, 100 * (rank + 1) + i);
#endif

     auto data = v->getData();

     for (auto i = data.lowerOffset(); i <= data.upperOffset(); i++)
     {
#if INTERVAL
         std::cout << i << "@" << rank << "[" << 
map->getGlobalElement(i) << "] total " << data.size() << " -> " << 
data[i].lower() << " to " << data[i].upper() << std::endl;
#else
         std::cout << i << "@" << rank << "[" << 
map->getGlobalElement(i) << "] total " << data.size() << " -> " << 
data[i] << std::endl;
#endif
     }

     return 0;
}