[Trilinos-Users] Tpetra and boost::interval [Update]
Sensei
senseiwa at gmail.com
Tue Aug 4 06:01:12 EDT 2015
Dear all,
Following the advice from Eric, I've read Sacado_MP_VectorTraits.hpp and
implemented all traits I could see. Two traits, the ScalarTraits and
SerializableTraits.
I've added the operator<< into the Teuchos namespace, due to the last
error about operator<< not being defined, but I get the following errors:
/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> > >')
main.cpp:10:10: In file included from main.cpp:10:
/usr/local/include/trilinos/Tpetra_Vector.hpp:2:10: In file included
from /usr/local/include/trilinos/Tpetra_Vector.hpp:2:
/usr/local/include/trilinos/Tpetra_Vector_decl.hpp:69:12: In
instantiation of member function
'Tpetra::Vector<boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>,
boost::numeric::interval_lib::checking_strict<double> > >, int, int,
Kokkos::SerialNode>::describe' requested here
main.cpp:410:31: In instantiation of member function
'Tpetra::Vector<boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>,
boost::numeric::interval_lib::checking_strict<double> > >, int, int,
Kokkos::SerialNode>::Vector' requested here
ostream:195:20: Candidate function not viable: no known conversion from
'const boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>,
boost::numeric::interval_lib::checking_strict<double> > >' to 'const
void *' for 1st argument; take the address of the argument with &
ostream:179:20: Candidate function not viable: no known conversion from
'const boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>,
boost::numeric::interval_lib::checking_strict<double> > >' to
'std::__1::basic_ostream<char> &(*)(std::__1::basic_ostream<char> &)'
for 1st argument
ostream:180:20: Candidate function not viable: no known conversion from
'const boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>,
boost::numeric::interval_lib::checking_strict<double> > >' to
'basic_ios<char_type, traits_type> &(*)(basic_ios<char_type,
traits_type> &)' for 1st argument
...
...
...
ostream:827:1: Candidate template ignored: could not match 'const _CharT
*' against 'boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>,
boost::numeric::interval_lib::checking_strict<double> > >'
ostream:1065:1: Candidate template ignored: could not match
'basic_string' against 'interval'
ostream:1082:1: Candidate template ignored: could not match 'shared_ptr'
against 'interval'
ostream:1051:5: Candidate template ignored: disabled by 'enable_if'
[with _Stream = std::__1::basic_ostream<char> &, _Tp =
boost::numeric::interval<double,
boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>,
boost::numeric::interval_lib::checking_strict<double> > >]
ostream:1089:1: Candidate template ignored: could not match 'bitset'
against 'interval'
iomanip:361:1: Candidate template ignored: could not match '__iom_t8'
against 'interval'
iomanip:476:5: Candidate template ignored: could not match '__iom_t10'
against 'interval'
complex:1514:1: Candidate template ignored: could not match 'complex'
against 'interval'
The code is below.
Any hints on how to proceed is very, very welcome. My view of Tpetra is
to allow general types to be used, not only PODs, so I'm willing to try.
Thanks!
#define INTERVAL 1
#include <iostream>
#include <limits>
#include <cstring>
#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>
#include <Teuchos_SerializationTraits.hpp>
typedef boost::numeric::interval<double> scalar;
bool isempty(const scalar& x)
{
return std::fabs(x.upper() - x.lower()) <
std::numeric_limits<double>::epsilon();
}
scalar exp(const scalar& x)
{
if (isempty(x)) return scalar::empty();
return scalar(std::exp(x.lower()), std::exp(x.upper()));
}
scalar log(const scalar& x)
{
if (isempty(x) || (x.upper() <=
std::numeric_limits<double>::epsilon())) return scalar::empty();
double l = !(x.lower() > std::numeric_limits<double>::epsilon()) ?
(-1.0 * std::numeric_limits<double>::infinity()) : std::log(x.lower());
return scalar(l, std::log(x.upper()));
}
scalar pow(const scalar& x, const scalar& y)
{
return exp(y * log(x));
}
namespace Teuchos
{
std::ostream& operator<< (std::ostream &s, const scalar &v)
{
s << "[" << v.lower() << "," << v.upper() << "]";
return s;
}
// Teuchos::FancyOStream& operator<< (Teuchos::FancyOStream &s, const
scalar &v)
// {
// return s;
// }
template <>
struct ScalarTraits<scalar>
{
//! Mandatory typedef for result of magnitude
typedef double 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 double magnitude(scalar a)
{
return std::fabs(a.upper() - a.lower());
}
//! Returns representation of zero for this scalar type.
static inline scalar zero()
{
return scalar(0.0, 0.0);
}
//! 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()
{
// SHOULD THIS BE THE ARCHETYPE FOR REAL/IMAG/CONJUGATE ?
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);
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)
{
// Defined by boost
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 exp(y * log(x));
}
};
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++;
}
}
};
}
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;
}
More information about the Trilinos-Users
mailing list