[Trilinos-Users] Tpetra and boost::interval [Update]
Sensei
senseiwa at gmail.com
Thu Aug 6 06:58:04 EDT 2015
Thanks to Mark & Eric for their hints, they've been very useful!
Right now, I am able to compile this project. The most basic operation
works, just assigning numbers and being able to retrieve them, and of
course, print and randomize.
My implementation deals with Teuchos traits for serialization and scalar
operations, and Kokkos operations, following
Trilinos/packages/kokkos/core/src/Kokkos_Complex.hpp.
I'm getting errors now trying to use the dot product, and creating a CRS
matrix. I haven't clear if the magnitude type should be a double or an
interval yet, however, I'm moving forward.
The source code is at the bottom.
The errors for the dot product seem to regard a failed conversion
between my scalar type (interval) and bool, which is weird to me: why an
inner product needs a boolean? It should be IMHO a morphism between the
vector space over intervals and the interval field. But I don't know the
inner workings of Tpetra and Kokkos.
I know I'm missing something here. As usual, I'm just asking for directions!
Thanks for any pointers!
Here's the complete error log following.
/usr/local/include/trilinos/Teuchos_CommHelpers.hpp:1017:26: Value of
type 'boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > >' is not
contextually convertible to 'bool'
main.cpp:8:10: In file included from main.cpp:8:
/usr/local/include/trilinos/Tpetra_DefaultPlatform.hpp:46:10: In file
included from /usr/local/include/trilinos/Tpetra_DefaultPlatform.hpp:46:
/usr/local/include/trilinos/Tpetra_ConfigDefs.hpp:137:10: In file
included from /usr/local/include/trilinos/Tpetra_ConfigDefs.hpp:137:
/usr/local/include/trilinos/Teuchos_CommHelpers.hpp:1087:13: In
instantiation of member function 'Teuchos::MixMaxUtilities::AND<true,
int, boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > > >::andOp'
requested here
/usr/local/include/trilinos/Teuchos_CommHelpers.hpp:1131:18: In
instantiation of member function 'Teuchos::ANDValueReductionOp<int,
boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > > >::reduce'
requested here
/usr/local/include/trilinos/Teuchos_CommHelpers.hpp:1429:14: In
instantiation of function template specialization
'Teuchos::createOp<int, boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > > >' requested here
/usr/local/include/trilinos/Teuchos_CommHelpers.hpp:1860:3: In
instantiation of function template specialization
'Teuchos::reduceAll<int, boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > > >' requested here
/usr/local/include/trilinos/Tpetra_Vector_def.hpp:171:16: In
instantiation of function template specialization
'Teuchos::reduceAll<int, boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > > >' requested here
main.cpp:507:8: In instantiation of member function
'Tpetra::Vector<boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > >, int, int,
Kokkos::SerialNode>::norm2' requested here
/usr/local/include/trilinos/Teuchos_CommHelpers.hpp:1017:41: Invalid
operands to binary expression ('boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > >' and 'const
boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::save_state<boost::numeric::interval_lib::rounded_transc_std<double,
boost::numeric::interval_lib::rounded_arith_std<double,
boost::numeric::interval_lib::rounding_control<double> > > >,
boost::numeric::interval_lib::checking_base<double> > >')
main.cpp:8:10: In file included from main.cpp:8:
/usr/local/include/trilinos/Tpetra_DefaultPlatform.hpp:46:10: In file
included from /usr/local/include/trilinos/Tpetra_DefaultPlatform.hpp:46:
/usr/local/include/trilinos/Tpetra_ConfigDefs.hpp:137:10: In file
included from /usr/local/include/trilinos/Tpetra_ConfigDefs.hpp:137:
#include <iostream>
#include <limits>
#include <cstring>
#include <mpi.h>
#include <boost/numeric/interval.hpp>
#include <Tpetra_DefaultPlatform.hpp>
#include <Tpetra_Vector.hpp>
#include <Tpetra_CrsMatrix.hpp>
#include <Tpetra_Version.hpp>
#include <Teuchos_GlobalMPISession.hpp>
#include <Teuchos_oblackholestream.hpp>
#include <Teuchos_SerializationTraits.hpp>
typedef boost::numeric::interval<double,
boost::numeric::interval_lib::policies<
boost::numeric::interval_lib::save_state<
boost::numeric::interval_lib::rounded_transc_std<double>
>,
boost::numeric::interval_lib::checking_base<double>
>
> scalar;
// Stream
namespace boost
{
namespace numeric
{
template <typename T, typename P>
std::ostream& operator<< (std::ostream& s, const interval<T,P>& v)
{
s << "[" << v.lower() << "," << v.upper() << "]";
return s;
}
}
}
// Traits
namespace Teuchos
{
template <>
struct ScalarTraits<scalar>
{
//! Mandatory typedef for result of magnitude
typedef scalar magnitudeType;
//! Typedef for half precision
typedef double halfPrecision;
//! Typedef for double precision
typedef double doublePrecision;
//! Determines if scalar type is std::complex
static const bool isComplex = false;
//! Determines if scalar type is an ordinal type
static const bool isOrdinal = false;
//! Determines if scalar type supports relational operators such
as <, >, <=, >=.
static const bool isComparable = true;
/** \brief Determines if scalar type have machine-specific
parameters
* (i.e. eps(), sfmin(), base(), prec(), t(), rnd(), emin(),
rmin(), emax(),
* rmax() are supported).
*/
static const bool hasMachineParameters = false;
//! Returns relative machine precision.
static inline double eps()
{
return std::numeric_limits<double>::epsilon();
}
//! Returns safe minimum (sfmin), such that 1/sfmin does not
overflow.
static inline double sfmin()
{
return std::numeric_limits<double>::min();
}
//! Returns the base of the machine.
static inline double base()
{
return std::numeric_limits<double>::radix;
}
//! Returns \c eps*base.
static inline double prec()
{
return eps() * base();
}
//! Returns the number of (base) digits in the mantissa.
static inline double t()
{
return std::numeric_limits<double>::digits;
}
//! Returns 1.0 when rounding occurs in addition, 0.0 otherwise
static inline double rnd()
{
return (std::numeric_limits<double>::round_style ==
std::round_to_nearest ? double(1.0) : double(0.0));
}
//! Returns the minimum exponent before (gradual) underflow.
static inline double emin()
{
return std::numeric_limits<double>::min_exponent;
}
//! Returns the underflow threshold - \c base^(emin-1)
static inline double rmin()
{
return std::numeric_limits<double>::min();
}
//! Returns the largest exponent before overflow.
static inline double emax()
{
return std::numeric_limits<double>::max_exponent;
}
//! Overflow theshold - \c (base^emax)*(1-eps)
static inline double rmax()
{
return std::numeric_limits<double>::max();
}
//! Returns the magnitudeType of the scalar type \c a.
static inline magnitudeType magnitude(scalar a)
{
return scalar(0.0, std::fabs(a.upper() - a.lower()));
}
//! Returns representation of zero for this scalar type.
static inline scalar zero()
{
return scalar(-123.4, 432.1);
}
//! Returns representation of one for this scalar type.
static inline scalar one()
{
return scalar(-1.0, +1.0);
}
//! Returns the real part of the scalar type \c a.
static inline scalar real(scalar a)
{
return a;
}
//! Returns the imaginary part of the scalar type \c a.
static inline scalar imag(scalar a)
{
return zero();
}
//! Returns the conjugate of the scalar type \c a.
static inline scalar conjugate(scalar a)
{
return a;
}
//! Returns a number that represents NaN.
static inline scalar nan()
{
return scalar(dbl_nan, dbl_nan);
}
//! Returns <tt>true</tt> if <tt>x</tt> is NaN or Inf.
static inline bool isnaninf(const scalar& x)
{
return generic_real_isnaninf(x.lower()) ||
generic_real_isnaninf(x.upper());
}
//! Seed the random number generator returned by <tt>random()</tt>.
static inline void seedrandom(unsigned int s)
{
// Copied from teuchos
std::srand(s);
#if __APPLE__
random();
#endif
}
//! Returns a random number (between -one() and +one()) of this
scalar type.
static inline scalar random()
{
// Copied from teuchos, now we just have two random numbers
double d1, d2;
d1 = -1.0 + 2.0 * ((double) std::rand() / RAND_MAX);
d2 = -1.0 + 2.0 * ((double) std::rand() / RAND_MAX);
if (d2 < d1) std::swap(d1, d2);
return scalar(d1, d2);
}
//! Returns the name of this scalar type.
static inline std::string name()
{
return "scalar aka 'boost::numeric::interval<double>'";
}
//! Returns a number of magnitudeType that is the square root of
this scalar type \c x.
static inline scalar squareroot(scalar x)
{
return boost::numeric::sqrt(x);
}
//! Returns the result of raising one scalar \c x to the power \c y.
static inline scalar pow(scalar x, scalar y)
{
return boost::numeric::exp(y * boost::numeric::log(y));
}
};
template <typename Ordinal>
class SerializationTraits<Ordinal, scalar> {
public:
//! @name Serialization type selection
//@{
/// \brief Whether the type T supports direct serialization.
///
/// See the class documentation for definitions of "direct" and
/// "indirect" serialization.
static const bool supportsDirectSerialization = false;
//@}
//! @name Direct serialization functions (not defined if
supportsDirectSerialization==false)
//@{
/** \brief Return the number of bytes for <tt>count</tt> objects. */
static Ordinal fromCountToDirectBytes(const Ordinal count)
{
return 2 * sizeof(double) * count;
}
/** \brief Convert the pointer type to <tt>char*</tt>. */
static char* convertToCharPtr(scalar* ptr)
{
void *p = ptr;
return static_cast<char*>(p);
}
/** \brief Convert the pointer type to <tt>const char*</tt>. */
static const char* convertToCharPtr( const scalar* ptr )
{
const void *p = ptr;
return static_cast<const char*>(p);
}
/** \brief Return the number of objects for <tt>bytes</tt> of
storage. */
static Ordinal fromDirectBytesToCount(const Ordinal bytes)
{
return bytes / (2 * sizeof(double));
}
/** \brief Convert the pointer type from <tt>char*</tt>. */
static scalar* convertFromCharPtr( char* ptr )
{
void *p = ptr;
return static_cast<scalar*>(p);
}
/** \brief Convert the pointer type from <tt>char*</tt>. */
static const scalar* convertFromCharPtr( const char* ptr )
{
const void *p = ptr;
return static_cast<const scalar*>(p);
}
//@}
//! @name Indirect serialization functions (always defined and
supported)
//@{
/** \brief Return the number of bytes for <tt>count</tt> objects. */
static Ordinal fromCountToIndirectBytes(const Ordinal count,
const scalar buffer[])
{
return 2 * sizeof(double) * count;
}
/** \brief Serialize to an indirect <tt>char[]</tt> buffer.
*
* \param count
* [in] The number of objects to serialize.
* \param buffer
* [in] The objects to serialize.
* \param bytes
* [in] Number of bytes in <tt>charBuffer[]</tt>
* \param charBuffer
* [out] Array (length <tt>bytes</tt>) containing the
serialized objects.
*
* <b>Preconditions:</b><ul>
* <li><tt>bytes==fromCountToIndirectBytes(count)</tt>
* </ul>
*/
static void serialize (const Ordinal count,
const scalar buffer[],
const Ordinal bytes,
char charBuffer[])
{
std::size_t b = 2 * sizeof(double) * count, i = 0;
if (b > bytes) throw "Serialization needs more bytes!";
while (i < b)
{
double tmp[2];
tmp[0] = buffer[i].lower();
tmp[1] = buffer[i].upper();
std::memcpy(charBuffer + i, tmp, 2 * sizeof(double));
i += 2 * sizeof(double);
}
}
/** \brief Return the number of objects for <tt>bytes</tt> of
storage. */
static Ordinal fromIndirectBytesToCount(const Ordinal bytes,
const char charBuffer[])
{
return bytes / (2 * sizeof(double));
}
/** \brief Deserialize from an indirect <tt>char[]</tt> buffer.
*
* \param bytes
* [in] Number of bytes in <tt>charBuffer[]</tt>
* \param charBuffer
* [in] Array (length <tt>bytes</tt>) containing the
serialized objects.
* \param count
* [in] The number of objects to deserialize.
* \param buffer
* [out] The deserialized objects.
* <b>Preconditions:</b><ul>
* <li><tt>count==fromIndirectBytesToCount(bytes)</tt>
* </ul>
*/
static void deserialize (const Ordinal bytes,
const char charBuffer[],
const Ordinal count,
scalar buffer[])
{
std::size_t c = bytes / (2 * sizeof(double)), i = 0;
if (c > count) throw "Deserialization needs more bytes!";
while (i < c)
{
double tmp[2];
tmp[0] = buffer[i].lower();
tmp[1] = buffer[i].upper();
std::memcpy(tmp, charBuffer + i, 2 * sizeof(double));
buffer[i] = scalar(tmp[0], tmp[1]);
i++;
}
}
};
}
// Operations
namespace Kokkos
{
//! Binary + operator for scalar.
scalar operator + (const scalar& x, const scalar& y)
{
return x + y;
}
//! Unary + operator for scalar.
scalar operator + (const scalar& x)
{
return x;
}
//! Binary - operator for scalar.
scalar operator - (const scalar& x, const scalar& y)
{
return x - y;
}
//! Unary - operator for scalar.
scalar operator - (const scalar& x)
{
return -x;
}
//! Binary * operator for scalar.
scalar operator * (const scalar& x, const scalar& y)
{
return x * y;
}
//! Imaginary part of a scalar.
scalar imag (const scalar& x)
{
return scalar(0.0, 0.0);
}
//! Real part of a scalar.
scalar real (const scalar& x)
{
return x;
}
//! Absolute value (magnitude) of a scalar.
double abs (const scalar& x)
{
return std::fabs(x.upper() - x.lower());
}
//! Conjugate of a scalar number.
scalar conj (const scalar& x)
{
return x;
}
//! Binary operator / for scalar and real numbers
scalar operator / (const scalar& x, const double& y)
{
return x / y;
}
//! Binary operator / for scalar.
scalar operator / (const scalar& x, const scalar& y)
{
return x / y;
}
//! Equality operator for two scalars.
bool operator == (const scalar& x, const scalar& y)
{
return x == y;
}
//! Inequality operator for two scalars.
bool operator != (const scalar& x, const scalar& y)
{
return x != y;
}
std::ostream& operator << (std::ostream& os, const scalar& x)
{
os << "[" << x.lower() << "," << x.upper() << "]";
return os;
}
// NO C++11 FOR NOW
void atomic_add (scalar* const dest, const scalar src)
{
*dest += src;
}
void atomic_assign (scalar* const dest, const scalar src)
{
*dest = src;
}
}
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;
auto v = Teuchos::rcp(new Tpetra::Vector<scalar>(map));
auto w = Teuchos::rcp(new Tpetra::Vector<scalar>(map));
auto p = Teuchos::rcp(new Tpetra::Vector<scalar>(map));
//auto A = Teuchos::rcp(new Tpetra::CrsMatrix<scalar>(map, 0));
v->putScalar(scalar(-3.1415, 2.7182));
w->putScalar(scalar(-2.0, 2.0 ));
v->print(std::cout);
w->print(std::cout);
//A->setAllToScalar(scalar(-1234, 1234));
v->norm2();
if ((false)) // PRINT V
{
std::cout << "gets in here" << std::endl;
auto data = v->getData();
for (auto i = data.lowerOffset(); i <= data.upperOffset(); i++)
{
std::cout << "V " << i << "@" << rank << "[" <<
map->getGlobalElement(i) << "] total " << data.size() << " -> " <<
data[i] << std::endl;
}
}
comm->barrier();
if ((true)) // SUM
{
std::cout << "SUM" << std::endl;
auto datav = v->getData();
auto dataw = w->getData();
auto datap = p->getData();
for (auto i = datav.lowerOffset(); i <= datav.upperOffset(); i++)
{
p->replaceLocalValue(i, datav[i] + dataw[i]);
}
}
comm->barrier();
//auto q = v->dot(*w);
if ((false)) // PRINT
{
auto data = p->getData();
for (auto i = data.lowerOffset(); i <= data.upperOffset(); i++)
{
std::cout << "P " << i << "@" << rank << "[" <<
map->getGlobalElement(i) << "] total " << data.size() << " -> " <<
data[i] << std::endl;
}
}
return 0;
}
More information about the Trilinos-Users
mailing list