[Trilinos-Users] [EXTERNAL] Re: Tpetra and boost::interval [Update]
Phipps, Eric T
etphipp at sandia.gov
Tue Aug 4 10:42:25 EDT 2015
I think you are defining operator<< in the wrong namespace. I think you
want instead:
namespace boost {
namespace numeric {
template <typename T, typename P, typename S>
std::ostream& operator<< (ostream& s, const interval<T,P,S>& v) {
s << "[" << v.lower() << "," << v.upper() << "]";
return s;
}
}
}
That way the relevant overload will be found by argument-dependent lookup
(the compiler looks in the namespace(s) of the argument(s)).
I am somewhat surprised interval doesn¹t already provide such an overload.
-Eric
On 8/4/15, 4:01 AM, "Sensei" <senseiwa at gmail.com> wrote:
>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::round
>ed_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::round
>ed_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::round
>ed_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::round
>ed_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::round
>ed_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::round
>ed_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::round
>ed_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::round
>ed_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