[Trilinos-Users] [EXTERNAL] Re: Tpetra and boost::interval [Update]
Phipps, Eric T
etphipp at sandia.gov
Tue Aug 4 10:44:33 EDT 2015
Actually I think I messed that up as interval only has two template
parameters not three:
namespace boost {
namespace numeric {
template <typename T, typename P>
std::ostream& operator<< (ostream& s, const interval<T,P>& v) {
s << "[" << v.lower() << "," << v.upper() << "]";
return s;
}
}
}
-Eric
On 8/4/15, 8:42 AM, "Phipps, Eric T" <etphipp at sandia.gov> wrote:
>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::roun
>>d
>>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::roun
>>d
>>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::roun
>>d
>>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::roun
>>d
>>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::roun
>>d
>>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::roun
>>d
>>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::roun
>>d
>>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::roun
>>d
>>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;
>>}
>>
>
>_______________________________________________
>Trilinos-Users mailing list
>Trilinos-Users at trilinos.org
>https://trilinos.org/mailman/listinfo/trilinos-users
More information about the Trilinos-Users
mailing list