From 8edc8e8175d367228a84ea43abc11fd233a253a4 Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Tue, 22 Apr 2025 20:47:39 +0100 Subject: [PATCH 1/7] Add daceAbsolute function to take DA absolute value --- core/dacemath.c | 22 +++++++++++++++++++--- core/include/dace/dacebase.h | 1 + 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/core/dacemath.c b/core/dacemath.c index 195f29d..92e872e 100755 --- a/core/dacemath.c +++ b/core/dacemath.c @@ -328,7 +328,7 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb) continue; const double c = ia->cc*ckon; - if(!(fabs(c) <= DACECom_t.eps)) + if(LIKELY(!(fabs(c) <= DACECom_t.eps))) { ib->cc = c; ib->ii = ia->ii; @@ -345,7 +345,7 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb) continue; const double c = ia->cc*ckon; - if(!(fabs(c) <= DACECom_t.eps)) + if(LIKELY(!(fabs(c) <= DACECom_t.eps))) { if(UNLIKELY(ib >= ibmax)) { @@ -653,6 +653,22 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc) * DACE intrinsic function routines *********************************************************************************/ +/*! Absolute value of a DA object. + Returns either the DA or the negative of the DA based on the constant part. + \param[in] ina Pointer to the DA object to operate on + \param[out] inc Pointer to the DA object to store the result in + \note This routine is aliasing safe, i.e. inc can be the same as ina. + \sa daceAbsoluteValue + \sa daceNorm + */ +void daceAbsolute(const DACEDA *ina, DACEDA *inc) +{ + daceCopy(ina, inc); + const double a0 = daceGetConstant(inc); + if(a0 < 0.0) + daceMultiplyDouble(inc, -1.0, inc); +} + /*! Truncate the constant part of a DA object to an integer. \param[in] ina Pointer to the DA object to operate on \param[out] inc Pointer to the DA object to store the result in @@ -684,7 +700,7 @@ void daceRound(const DACEDA *ina, DACEDA *inc) void daceModulo(const DACEDA *ina, const double p, DACEDA *inc) { daceCopy(ina, inc); - daceSetCoefficient0(inc, 0, fmod(daceGetConstant(inc),p)); + daceSetCoefficient0(inc, 0, fmod(daceGetConstant(inc), p)); } /*! Raise a DA object to the p-th power. diff --git a/core/include/dace/dacebase.h b/core/include/dace/dacebase.h index ae6c6c3..3939ba1 100644 --- a/core/include/dace/dacebase.h +++ b/core/include/dace/dacebase.h @@ -170,6 +170,7 @@ DACE_API void daceIntegrate(const unsigned int iint, const DACEDA REF(ina), DACE /******************************************************************************** * DACE intrinsic function routines *********************************************************************************/ +DACE_API void daceAbsolute(const DACEDA REF(ina), DACEDA REF(inc)); DACE_API void daceTruncate(const DACEDA REF(ina), DACEDA REF(inc)); DACE_API void daceRound(const DACEDA REF(ina), DACEDA REF(inc)); DACE_API void daceModulo(const DACEDA REF(ina), const double p, DACEDA REF(inc)); From 4d10f57cfd90631e506629fe11e54bef45c1e282 Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Tue, 22 Apr 2025 21:23:30 +0100 Subject: [PATCH 2/7] expose daceAbsolute in C++ interface. Can't use DA abs() to avoid confusion with builtin double abs() --- interfaces/cxx/DA.cpp | 27 +++++++++++++++++++++++++++ interfaces/cxx/include/dace/DA.h | 2 ++ 2 files changed, 29 insertions(+) diff --git a/interfaces/cxx/DA.cpp b/interfaces/cxx/DA.cpp index e5b4a4b..b860077 100644 --- a/interfaces/cxx/DA.cpp +++ b/interfaces/cxx/DA.cpp @@ -966,6 +966,21 @@ DA DA::trim(const unsigned int min, const unsigned int max) const{ return temp; } +DA DA::absolute() const{ +/*! Absolute value of a DA object. + Returns either the DA or the negative of the DA based on the constant part. + \return A new DA object with the absolute value DA. + \throw DACE::DACEException + \sa DA::abs + \sa DA::norm + */ + DA temp; + daceAbsolute(m_index, temp.m_index); + if(daceGetError()) DACEException(); + + return temp; +} + DA DA::trunc() const{ /*! Truncate the constant part of a DA object to an integer. The result is copied in a new DA object. @@ -2070,6 +2085,18 @@ DA trim(const DA &da, const unsigned int min, const unsigned int max){ */ return da.trim(min,max);} +DA absolute(const DA &da){ +/*! Absolute value of a DA object. + Returns either the DA or the negative of the DA based on the constant part. + \param[in] da a given DA object. + \return A new DA object with the absolute value DA. + \throw DACE::DACEException + \sa DA::absolute + \sa DA::abs + \sa DA::norm + */ + return da.absolute();} + DA trunc(const DA &da){ /*! Truncate the constant part of a DA object to an integer. The result is copied in a new DA object. diff --git a/interfaces/cxx/include/dace/DA.h b/interfaces/cxx/include/dace/DA.h index 20d96c0..d511c54 100644 --- a/interfaces/cxx/include/dace/DA.h +++ b/interfaces/cxx/include/dace/DA.h @@ -162,6 +162,7 @@ class DACE_API DA DA integ(const std::vector ind) const; //!< Integral with respect to given variables DA trim(const unsigned int min, const unsigned int max = DA::getMaxOrder()) const; //!< Trim the coefficients only include orders between min and max, inclusively + DA absolute() const; //!< Absolute value of DA based on constant part DA trunc() const; //!< Truncate the constant part to an integer DA round() const; //!< Round the constant part to an integer DA mod(const double p) const; //!< Modulo of the constant part @@ -272,6 +273,7 @@ DACE_API DA deriv(const DA &da, const std::vector ind); DACE_API DA integ(const DA &da, const unsigned int i); DACE_API DA integ(const DA &da, const std::vector ind); DACE_API DA trim(const DA &da, const unsigned int min, const unsigned int max = DA::getMaxOrder()); +DACE_API DA absolute(const DA &da); DACE_API DA trunc(const DA &da); DACE_API DA round(const DA &da); DACE_API DA mod(const DA &da, const double p); From d603fc3ca0cc122d4be937b63fea3c712864ce3f Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Tue, 22 Apr 2025 21:28:28 +0100 Subject: [PATCH 3/7] add double absolute() overload extension --- interfaces/cxx/MathExtension.cpp | 7 +++++++ interfaces/cxx/include/dace/MathExtension.h | 1 + 2 files changed, 8 insertions(+) diff --git a/interfaces/cxx/MathExtension.cpp b/interfaces/cxx/MathExtension.cpp index dd530a5..a3e9dc9 100644 --- a/interfaces/cxx/MathExtension.cpp +++ b/interfaces/cxx/MathExtension.cpp @@ -35,6 +35,13 @@ namespace DACE{ +double absolute(const double x){ +/*! Absolute value. + \param[in] x Function argument. + */ + return std::abs(x); +} + double cons(const double x){ /*! Constant part. For double type this is just x. \param[in] x Function argument. diff --git a/interfaces/cxx/include/dace/MathExtension.h b/interfaces/cxx/include/dace/MathExtension.h index 69eff5b..f2bd00a 100644 --- a/interfaces/cxx/include/dace/MathExtension.h +++ b/interfaces/cxx/include/dace/MathExtension.h @@ -31,6 +31,7 @@ namespace DACE{ +DACE_API double absolute(const double x); //!< Absolute value DACE_API double cons(const double x); //!< Constant part (i.e. the value x) DACE_API double logb(const double x, const double b = 10.0); //!< Logarithm relative to base b DACE_API double isrt(const double x); //!< Inverse square root From b5e0c2f1842913c55fab660180f7c5b5a79212a5 Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Fri, 25 Apr 2025 17:17:08 +0100 Subject: [PATCH 4/7] The great abs() removal and AlgebraicVector cleanup (part 1) --- core/dacecompat.c | 8 +- core/dacenorm.c | 9 - core/include/dace/dacebase.h | 1 - interfaces/cxx/AlgebraicVector.cpp | 2 +- interfaces/cxx/DA.cpp | 63 +- interfaces/cxx/MathExtension.cpp | 7 + interfaces/cxx/include/dace/AlgebraicVector.h | 30 +- .../cxx/include/dace/AlgebraicVector_t.h | 878 +++++++++++------- interfaces/cxx/include/dace/DA.h | 23 +- interfaces/cxx/include/dace/MathExtension.h | 1 + 10 files changed, 615 insertions(+), 407 deletions(-) diff --git a/core/dacecompat.c b/core/dacecompat.c index a45453f..ef69589 100644 --- a/core/dacecompat.c +++ b/core/dacecompat.c @@ -63,14 +63,16 @@ void dacepek(const DACEDA *ina, const unsigned int jj[], double *cjj) } /*! Compute absolute value of a DA object. + Same as daceNorm(ina, 0). \param[in] ina Pointer to DA object to take absolute value of \param[out] anorm Pointer where to store the absolute value - \deprecated Has been replaced by daceAbsoluteValue() - \sa daceAbsoluteValue() + \deprecated Has been replaced by daceNorm() and daceAbsolute() + \sa daceNorm() + \sa daceAbsolute() */ void daceabs(const DACEDA *ina, double *anorm) { - *anorm = daceAbsoluteValue(ina); + *anorm = daceNorm(ina, 0); } /*! Compute absolute value of a DA object. diff --git a/core/dacenorm.c b/core/dacenorm.c index 2254a74..804941f 100644 --- a/core/dacenorm.c +++ b/core/dacenorm.c @@ -42,15 +42,6 @@ * DACE norm and norm estimation routines *********************************************************************************/ -/*! Compute the absolute value (maximum coefficient norm) of a DA object. - \param[in] ina Pointer to the DA object to take absolute value of - \return The absolute value of ina -*/ -double daceAbsoluteValue(const DACEDA *ina) -{ - return daceNorm(ina, 0); -} - /*! Compute a norm of a DA object. \param[in] ina Pointer to the DA object to take norm of \param[in] ityp Type of norm to compute. diff --git a/core/include/dace/dacebase.h b/core/include/dace/dacebase.h index 3939ba1..68f177a 100644 --- a/core/include/dace/dacebase.h +++ b/core/include/dace/dacebase.h @@ -214,7 +214,6 @@ DACE_API void dacePsiFunction(const DACEDA REF(ina), const unsigned int n, DACED /******************************************************************************** * DACE norm and norm estimation routines *********************************************************************************/ -DACE_API double daceAbsoluteValue(const DACEDA REF(ina)); DACE_API double daceNorm(const DACEDA REF(ina), const unsigned int ityp); DACE_API void daceOrderedNorm(const DACEDA REF(ina), const unsigned int ivar, const unsigned int ityp, double onorm[]); DACE_API void daceEstimate(const DACEDA REF(ina), const unsigned int ivar, const unsigned int ityp, double c[], double err[], const unsigned int nc); diff --git a/interfaces/cxx/AlgebraicVector.cpp b/interfaces/cxx/AlgebraicVector.cpp index 388a4c0..3fbf82f 100644 --- a/interfaces/cxx/AlgebraicVector.cpp +++ b/interfaces/cxx/AlgebraicVector.cpp @@ -37,7 +37,7 @@ #include "dace/AlgebraicVector.h" #include "dace/AlgebraicVector_t.h" -namespace DACE{ +namespace DACE { /*********************************************************************************** * Coefficient access routines diff --git a/interfaces/cxx/DA.cpp b/interfaces/cxx/DA.cpp index b860077..2298050 100644 --- a/interfaces/cxx/DA.cpp +++ b/interfaces/cxx/DA.cpp @@ -1534,25 +1534,12 @@ unsigned int DA::size() const{ \return The number of non-zero coefficients stored in the DA object. \throw DACE::DACEException */ - unsigned int res; - res=daceGetLength(m_index); + const unsigned int res = daceGetLength(m_index); if(daceGetError()) DACEException(); return res; } -double DA::abs() const{ -/*! Compute the max norm of a DA object. - \return A double corresponding to the result of the operation. - \throw DACE::DACEException - */ - double c; - c=daceAbsoluteValue(m_index); - if(daceGetError()) DACEException(); - - return c; -} - double DA::norm(const unsigned int type) const{ /*! Compute different types of norms for a DA object. \param[in] type type of norm to be computed. Possible norms are:\n @@ -1562,8 +1549,7 @@ double DA::norm(const unsigned int type) const{ \return A double corresponding to the result of the operation. \throw DACE::DACEException */ - double c; - c=daceNorm(m_index, type); + const double c = daceNorm(m_index, type); if(daceGetError()) DACEException(); return c; @@ -1582,7 +1568,7 @@ std::vector DA::orderNorm(const unsigned int var, const unsigned int typ \throw DACE::DACEException */ std::vector v(daceGetMaxOrder()+1); - daceOrderedNorm(m_index, var, type, v.data()); // Note: v.data() is C++11 + daceOrderedNorm(m_index, var, type, v.data()); if(daceGetError()) DACEException(); return v; @@ -1604,7 +1590,7 @@ std::vector DA::estimNorm(const unsigned int var, const unsigned int typ \note If estimation is not possible, zero is returned for all requested orders. */ std::vector v(nc+1); - daceEstimate(m_index, var, type, v.data(), NULL ,nc); // Note: v.data() is C++11 + daceEstimate(m_index, var, type, v.data(), NULL, nc); if(daceGetError()) DACEException(); return v; @@ -2582,15 +2568,6 @@ unsigned int size(const DA &da){ */ return da.size();} -double abs(const DA &da){ -/*! Compute the max norm of a DA object. - \param[in] da a given DA object. - \return A double corresponding to the result of the operation. - \throw DACE::DACEException - \sa DA::abs - */ - return da.abs();} - double norm(const DA &da, unsigned int type){ /*! Compute different types of norms for a DA object. \param[in] da a given DA object. @@ -2759,6 +2736,38 @@ void write(const DA &da, std::ostream &os){ */ return da.write(os);} +namespace abs_cons { + double abs(const DA &da){ +/*! Absolute value of constant part. + \param[in] da a given DA object. + \throw DACE::DACEException + \sa DA::cons + */ + return std::abs(da.cons()); + } +} + +namespace abs_max { + double abs(const DA &da){ +/*! Largest coefficient in absolute value. + \param[in] da a given DA object. + \throw DACE::DACEException + \sa DA::norm + */ + return da.norm(0); + } +} + +namespace abs_sum { + double abs(const DA &da){ +/*! Sum of absolute values of all coefficients. + \param[in] da a given DA object. + \throw DACE::DACEException + \sa DA::norm + */ + return da.norm(1); + } +} // static class variables const unsigned int storedDA::headerSize = daceBlobSize(NULL); diff --git a/interfaces/cxx/MathExtension.cpp b/interfaces/cxx/MathExtension.cpp index a3e9dc9..e381542 100644 --- a/interfaces/cxx/MathExtension.cpp +++ b/interfaces/cxx/MathExtension.cpp @@ -64,6 +64,13 @@ double isrt(const double x){ return 1.0/std::sqrt(x); } +double icbrt(const double x){ +/*! Inverse cube root 1/cbrt(x). + \param[in] x Function argument. + */ + return 1.0/std::cbrt(x); +} + double sqr(const double x){ /*! Square of x. \param[in] x Function argument. diff --git a/interfaces/cxx/include/dace/AlgebraicVector.h b/interfaces/cxx/include/dace/AlgebraicVector.h index b7b113c..546f670 100644 --- a/interfaces/cxx/include/dace/AlgebraicVector.h +++ b/interfaces/cxx/include/dace/AlgebraicVector.h @@ -88,11 +88,25 @@ template class AlgebraicVector : public std::vector /*********************************************************************************** * Math routines ************************************************************************************/ - // Included also in the std library cmath + AlgebraicVector absolute() const; //!< Element-wise absolute value function. + AlgebraicVector trunc() const; //!< Element-wise truncation. + AlgebraicVector round() const; //!< Element-wise rounding. + AlgebraicVector mod() const; //!< Element-wise modulo. AlgebraicVector pow(const int p) const; //!< Element-wise exponentiation to (integer) power. + AlgebraicVector pow(const double p) const; //!< Element-wise exponentiation to (double) power. + AlgebraicVector root(const int p = 2) const; //!< Element-wise p-th root. + AlgebraicVector minv() const; //!< Element-wise multiplicative inverse. + AlgebraicVector sqr() const; //!< Element-wise square. AlgebraicVector sqrt() const; //!< Element-wise square root. + AlgebraicVector isrt() const; //!< Element-wise inverse square root. + AlgebraicVector cbrt() const; //!< Element-wise cube root. + AlgebraicVector icbrt() const; //!< Element-wise inverse cube root. + AlgebraicVector hypot(const AlgebraicVector &obj) const; //!< Element-wise hypotenuse. AlgebraicVector exp() const; //!< Element-wise exponential. AlgebraicVector log() const; //!< Element-wise natural logarithm. + AlgebraicVector logb(const double b = 10.0) const; //!< Element-wise logarithm wrt a given base. + AlgebraicVector log10() const; //!< Element-wise logarithm to base 10. + AlgebraicVector log2() const; //!< Element-wise logarithm to base 2. AlgebraicVector sin() const; //!< Element-wise sine. AlgebraicVector cos() const; //!< Element-wise cosine. AlgebraicVector tan() const; //!< Element-wise tangent. @@ -103,19 +117,10 @@ template class AlgebraicVector : public std::vector AlgebraicVector sinh() const; //!< Element-wise hyperbolic sine. AlgebraicVector cosh() const; //!< Element-wise hyperbolic cosine. AlgebraicVector tanh() const; //!< Element-wise hyperbolic tangent. - - // Available in cmath for double only with C++11 AlgebraicVector asinh() const; //!< Element-wise hyperbolic arcsine. AlgebraicVector acosh() const; //!< Element-wise hyperbolic arccosine. AlgebraicVector atanh() const; //!< Element-wise hyperbolic arctangent. - // Our own math extension routines not included in cmath - AlgebraicVector logb(const double b = 10.0) const; //!< Element-wise logarithm wrt a given base. - AlgebraicVector isrt() const; //!< Element-wise inverse square root. - AlgebraicVector sqr() const; //!< Element-wise square. - AlgebraicVector minv() const; //!< Element-wise multiplicative inverse. - AlgebraicVector root(const int p = 2) const; //!< Element-wise p-th root. - /*********************************************************************************** * Vector routines ************************************************************************************/ @@ -123,7 +128,7 @@ template class AlgebraicVector : public std::vector //!< Dot product (scalar product, inner product) of two vectors. template AlgebraicVector::returnType> cross(const AlgebraicVector &obj) const; //!< Cross product of two vectors of length 3. - T vnorm() const; //!< Euclidean vector norm (length). + T length() const; //!< Length of the vector in Euclidean norm. AlgebraicVector normalize() const; //!< Normalized vector of unit length along this vector. /*********************************************************************************** @@ -139,10 +144,9 @@ template class AlgebraicVector : public std::vector AlgebraicVector trim(const unsigned int min, const unsigned int max = DA::getMaxOrder()) const; //!< Trim the coefficients of each components to particular orders. DA only. AlgebraicVector invert() const; //!< Inverse function of the AlgebraicVector. DA only. - // XXX: define and add the norm estimation routines from DA including abs(), norm(p), convergence radius estimation + // XXX: define and add the norm estimation routines from DA including norm(p), convergence radius estimation // XXX: add jacobian routine (returns a DA matrix containing the jacobian) /* - double abs() const; //!< Maximum absolute value of all coefficients double norm(const unsigned int type = 0) const; //!< Different types of norms over all coefficients std::vector orderNorm(const unsigned int var = 0, const unsigned int type = 0) const; //!< Different types of norms over coefficients of each order separately diff --git a/interfaces/cxx/include/dace/AlgebraicVector_t.h b/interfaces/cxx/include/dace/AlgebraicVector_t.h index a6c26a9..138228c 100644 --- a/interfaces/cxx/include/dace/AlgebraicVector_t.h +++ b/interfaces/cxx/include/dace/AlgebraicVector_t.h @@ -45,42 +45,42 @@ #include "dace/AlgebraicMatrix_t.h" #endif /* WITH_ALGEBRAICMATRIX */ -namespace DACE{ +namespace DACE { /*********************************************************************************** * Constructors ************************************************************************************/ -template AlgebraicVector::AlgebraicVector() : std::vector(){ +template AlgebraicVector::AlgebraicVector() : std::vector() { /*! Default Constructor to create empty AlgebraicVector */ } -template AlgebraicVector::AlgebraicVector(size_t size) : std::vector(size){ +template AlgebraicVector::AlgebraicVector(size_t size) : std::vector(size) { /*! Constructor with size to allocate a vector of the given size with elements initialized using their default constructor. \param[in] size length of AlgebraicVector. */ } -template AlgebraicVector::AlgebraicVector(size_t size, const T &d) : std::vector(size, d){ +template AlgebraicVector::AlgebraicVector(size_t size, const T &d) : std::vector(size, d) { /*! Constructor with size and elements value to allocate a vector of the given size with elements initialized as copies of d. \param[in] size length of AlgebraicVector \param[in] d initial value for the elements */ } -template AlgebraicVector::AlgebraicVector(const std::vector &v) : std::vector(v){ +template AlgebraicVector::AlgebraicVector(const std::vector &v) : std::vector(v) { /*! Copy constructor to create a copy of any existing vector. \param[in] v vector to be copied into AlgebraicVector */ } -template AlgebraicVector::AlgebraicVector(std::initializer_list l) : std::vector(l){ +template AlgebraicVector::AlgebraicVector(std::initializer_list l) : std::vector(l) { /*! Constructor to create a vector from an initializer list. \param[in] l braced initializer list to be copied into the AlgebraicVector */ } -template AlgebraicVector::AlgebraicVector(const std::vector &v, size_t first, size_t last) : std::vector(v.begin()+first, v.begin()+last+1){ +template AlgebraicVector::AlgebraicVector(const std::vector &v, size_t first, size_t last) : std::vector(v.begin()+first, v.begin()+last+1) { /*! Extraction constructor to copy only a given range of elements from vector v. \param[in] v vector to be copied into AlgebraicVector \param[in] first index of the first element to be copied @@ -96,24 +96,21 @@ template AlgebraicVector::AlgebraicVector(const std::vector &v /*********************************************************************************** * Coefficient access routines ************************************************************************************/ -template AlgebraicVector AlgebraicVector::cons() const{ -/*! Return the constant part of a AlgebraicVector. - \return A AlgebraicVector of dimension 1 by size, where size is the - size of the AlgebraicVector. Each element contains the constant part of - the corresponding value included in the original AlgebraicVector. +template AlgebraicVector AlgebraicVector::cons() const { +/*! Return the constant parts of each element. + \return A AlgebraicVector containing the constant part of each element. */ using DACE::cons; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::extract(const size_t first, const size_t last) const -{ +template AlgebraicVector AlgebraicVector::extract(const size_t first, const size_t last) const { /*! Extracts elements from AlgebraicVector. \param[in] first index of first element to be extracted \param[in] last index of last element to be extracted @@ -127,7 +124,7 @@ template AlgebraicVector AlgebraicVector::extract(const size_t return AlgebraicVector(*this, first, last); } -template template AlgebraicVector::returnType> AlgebraicVector::concat(const std::vector &obj) const{ +template template AlgebraicVector::returnType> AlgebraicVector::concat(const std::vector &obj) const { /*! Append an AlgebraicVector to the end of the current one and return the new vector. \param[in] obj The AlgebraicVector to be appended. \return A new AlgebraicVector containing the elements of both vectors, @@ -148,14 +145,14 @@ template template AlgebraicVector AlgebraicVector AlgebraicVector::operator-() const{ +template AlgebraicVector AlgebraicVector::operator-() const { /*! Returns the additive inverse of the vector. \return A new AlgebraicVector, with the opposite sign. */ return -1.0*(*this); } -template template AlgebraicVector& AlgebraicVector::operator+=(const AlgebraicVector &obj){ +template template AlgebraicVector& AlgebraicVector::operator+=(const AlgebraicVector &obj) { /*! Add the given AlgebraicVector to ourselves. \param[in] obj An AlgebraicVector. \return A reference to ourselves. @@ -165,23 +162,23 @@ template template AlgebraicVector& AlgebraicVector if(size != obj.size()) throw std::runtime_error("DACE::AlgebraicVector::operator+=: Vectors must have the same length."); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator+=(const U &obj){ +template template AlgebraicVector& AlgebraicVector::operator+=(const U &obj) { /*! Add the given scalar to ourselves componentwise. \param[in] obj A scalar value. \return A reference to ourselves. */ const size_t size = this->size(); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator-=(const AlgebraicVector &obj){ +template template AlgebraicVector& AlgebraicVector::operator-=(const AlgebraicVector &obj) { /*! Subtract the given AlgebraicVector from ourselves. \param[in] obj An AlgebraicVector. \return A reference to ourselves. @@ -191,23 +188,23 @@ template template AlgebraicVector& AlgebraicVector if(size != obj.size()) throw std::runtime_error("DACE::AlgebraicVector::operator-=: Vectors must have the same length."); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator-=(const U &obj){ +template template AlgebraicVector& AlgebraicVector::operator-=(const U &obj) { /*! Subtract the given scalar from ourselves componentwise. \param[in] obj A scalar value. \return A reference to ourselves. */ const size_t size = this->size(); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator*=(const AlgebraicVector &obj){ +template template AlgebraicVector& AlgebraicVector::operator*=(const AlgebraicVector &obj) { /*! Multiply the given AlgebraicVector with ourselves componentwise. \param[in] obj An AlgebraicVector. \return A reference to ourselves. @@ -216,23 +213,23 @@ template template AlgebraicVector& AlgebraicVector if(size != obj.size()) throw std::runtime_error("DACE::AlgebraicVector::operator*=: Vectors must have the same length."); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator*=(const U &obj){ +template template AlgebraicVector& AlgebraicVector::operator*=(const U &obj) { /*! Multiply the given scalar with ourselves. \param[in] obj A scalar value. \return A reference to ourselves. */ const size_t size = this->size(); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator/=(const AlgebraicVector &obj){ +template template AlgebraicVector& AlgebraicVector::operator/=(const AlgebraicVector &obj) { /*! Divide ourselves by the given AlgebraicVector componentwise. \param[in] obj An AlgebraicVector. \return A reference to ourselves. @@ -242,38 +239,38 @@ template template AlgebraicVector& AlgebraicVector if(size != obj.size()) throw std::runtime_error("DACE::AlgebraicVector::operator/=: Vectors must have the same length."); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator/=(const U &obj){ +template template AlgebraicVector& AlgebraicVector::operator/=(const U &obj) { /*! Divide ourselves by the given scalar. \param[in] obj A scalar value. \return A reference to ourselves. */ const size_t size = this->size(); - for(size_t i=0; i template AlgebraicVector& AlgebraicVector::operator<<(const std::vector &obj){ +template template AlgebraicVector& AlgebraicVector::operator<<(const std::vector &obj) { /*! Append elements of vector obj to the end of ourself, converting the type to match ours if necessary. \param[in] obj Vector of elements to append. \return A reference to ourselves. */ const size_t size = obj.size(); - for(size_t i=0; i AlgebraicVector::returnType> operator+(const AlgebraicVector &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator+(const AlgebraicVector &obj1, const AlgebraicVector &obj2) { /*! Compute the addition between two AlgebraicVectors. \param[in] obj1 first AlgebraicVector. \param[in] obj2 second AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1+obj2). + \return A new AlgebraicVector, (obj1+obj2). \throw std::runtime_error */ if(obj1.size() != obj2.size()) @@ -281,42 +278,42 @@ template AlgebraicVector: const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator+(const AlgebraicVector &obj1, const V &obj2){ +template AlgebraicVector::returnType> operator+(const AlgebraicVector &obj1, const V &obj2) { /*! Compute the addition between a AlgebraicVector and a scalar value. \param[in] obj1 a AlgebraicVector. \param[in] obj2 a scalar value. - \return A new AlgebraicVector, containing the result of the operation (obj1+obj2). + \return A new AlgebraicVector, (obj1+obj2). */ const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator+(const U &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator+(const U &obj1, const AlgebraicVector &obj2) { /*! Compute the addition between a scalar value and a AlgebraicVector. \param[in] obj1 a scalar value. \param[in] obj2 a AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1+obj2). + \return A new AlgebraicVector, (obj1+obj2). */ const size_t size = obj2.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator-(const AlgebraicVector &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator-(const AlgebraicVector &obj1, const AlgebraicVector &obj2) { /*! Compute the subtraction between two AlgebraicVectors. \param[in] obj1 first AlgebraicVector. \param[in] obj2 second AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1-obj2). + \return A new AlgebraicVector, (obj1-obj2). \throw std::runtime_error */ if(obj1.size() != obj2.size()) @@ -324,42 +321,42 @@ template AlgebraicVector: const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator-(const AlgebraicVector &obj1, const V &obj2){ +template AlgebraicVector::returnType> operator-(const AlgebraicVector &obj1, const V &obj2) { /*! Compute the subtraction between a AlgebraicVector and a scalar value. \param[in] obj1 a AlgebraicVector. \param[in] obj2 a scalar value. - \return A new AlgebraicVector, containing the result of the operation (obj1-obj2). + \return A new AlgebraicVector, (obj1-obj2). */ const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator-(const U &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator-(const U &obj1, const AlgebraicVector &obj2) { /*! Compute the subtraction between a scalar value and a AlgebraicVector. \param[in] obj1 a scalar value. \param[in] obj2 a AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1-obj2). + \return A new AlgebraicVector, (obj1-obj2). */ const size_t size = obj2.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator*(const AlgebraicVector &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator*(const AlgebraicVector &obj1, const AlgebraicVector &obj2) { /*! Compute the element-wise multiplication between two AlgebraicVectors. \param[in] obj1 first AlgebraicVector. \param[in] obj2 second AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1*obj2). + \return A new AlgebraicVector, (obj1*obj2). \throw std::runtime_error */ if(obj1.size() != obj2.size()) @@ -367,42 +364,42 @@ template AlgebraicVector: const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator*(const AlgebraicVector &obj1, const V &obj2){ +template AlgebraicVector::returnType> operator*(const AlgebraicVector &obj1, const V &obj2) { /*! Compute the multiplication between a AlgebraicVector and a scalar value. \param[in] obj1 a AlgebraicVector. \param[in] obj2 a scalar value. - \return A new AlgebraicVector, containing the result of the operation (obj1*obj2). + \return A new AlgebraicVector, (obj1*obj2). */ const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator*(const U &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator*(const U &obj1, const AlgebraicVector &obj2) { /*! Compute the multiplication between a scalar value and a AlgebraicVector. \param[in] obj1 a scalar value. \param[in] obj2 a AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1*obj2). + \return A new AlgebraicVector, (obj1*obj2). */ const size_t size = obj2.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator/(const AlgebraicVector &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator/(const AlgebraicVector &obj1, const AlgebraicVector &obj2) { /*! Compute the element-wise division between two AlgebraicVectors. \param[in] obj1 first AlgebraicVector. \param[in] obj2 second AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1/obj2). + \return A new AlgebraicVector, (obj1/obj2). \throw std::runtime_error */ if(obj1.size() != obj2.size()) @@ -410,58 +407,58 @@ template AlgebraicVector: const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator/(const AlgebraicVector &obj1, const V &obj2){ +template AlgebraicVector::returnType> operator/(const AlgebraicVector &obj1, const V &obj2) { /*! Compute the division between a AlgebraicVector and a scalar value. \param[in] obj1 a AlgebraicVector. \param[in] obj2 a scalar value. - \return A new AlgebraicVector, containing the result of the operation (obj1/obj2). + \return A new AlgebraicVector, (obj1/obj2). */ const size_t size = obj1.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i AlgebraicVector::returnType> operator/(const U &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> operator/(const U &obj1, const AlgebraicVector &obj2) { /*! Compute the division between a scalar value and a AlgebraicVector. \param[in] obj1 a scalar value. \param[in] obj2 a AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation (obj1/obj2). + \return A new AlgebraicVector, (obj1/obj2). */ const size_t size = obj2.size(); AlgebraicVector::returnType> temp(size); - for(size_t i=0; i template typename PromotionTrait::returnType AlgebraicVector::dot(const AlgebraicVector &obj) const{ +template template typename PromotionTrait::returnType AlgebraicVector::dot(const AlgebraicVector &obj) const { /*! Compute the dot product with another AlgebraicVector. \param[in] obj the other AlgebraicVector. - \return A scalar value, containing the result of the operation. + \return A scalar value,. \throw std::runtime_error */ const size_t size = this->size(); if(size != obj.size()) throw std::runtime_error("DACE::AlgebraicVector::dot(): Vectors must have the same length."); - typename PromotionTrait::returnType temp = 0; - for(size_t i=0; i::returnType temp = 0.0; + for(size_t i=0; i template AlgebraicVector::returnType> AlgebraicVector::cross(const AlgebraicVector &obj) const{ +template template AlgebraicVector::returnType> AlgebraicVector::cross(const AlgebraicVector &obj) const { /*! Compute the cross product with another 3D AlgebraicVector. \param[in] obj The other AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation. + \return A new AlgebraicVector. \throw std::runtime_error */ if((this->size() != 3) || (obj.size() != 3)) @@ -479,340 +476,461 @@ template template AlgebraicVector AlgebraicVector AlgebraicVector::pow(const int p) const{ -/*! Elevate a AlgebraicVector to a given integer power. - The result is copied in a new AlgebraicVector. - \param[in] p power at which the AlgebraicVector is elevated. +template AlgebraicVector AlgebraicVector::absolute() const { +/*! Element-wise application of the absolute value function. + \return A new AlgebraicVector. + */ + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::trunc() const { +/*! Element-wise application of the truncation function. + \return A new AlgebraicVector. + */ + using std::trunc; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::round() const { +/*! Element-wise application of the round function. + \return A new AlgebraicVector. + */ + using std::round; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::mod() const { +/*! Element-wise application of the mod function. + \return A new AlgebraicVector. + */ + using std::trunc; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::pow(const int p) const { +/*! Element-wise application of the integer power function. + \param[in] p power to raise each element to. + \return A new AlgebraicVector. + */ + using std::pow; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::pow(const double p) const { +/*! Element-wise application of the double power function. + \param[in] p power to raise each element to. \return A new AlgebraicVector. */ using std::pow; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::sqrt() const{ -/*! Compute the square root of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::root(const int p) const { +/*! Element-wise application of the p-th root function. + \param[in] p root to be computed. + \return A new AlgebraicVector. + */ + using DACE::root; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::minv() const { +/*! Element-wise application of the multiplicative inverse function. + \return A new AlgebraicVector. + */ + using DACE::minv; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::sqr() const { +/*! Element-wise application of the square function. + \return A new AlgebraicVector. + */ + using DACE::sqr; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::sqrt() const { +/*! Element-wise application of the square root function. \return A new AlgebraicVector. */ using std::sqrt; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::exp() const{ -/*! Compute the exponent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::isrt() const { +/*! Element-wise application of the inverse square root function. \return A new AlgebraicVector. */ - using std::exp; + using DACE::isrt; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::log() const{ -/*! Compute the natural logarithm of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::cbrt() const { +/*! Element-wise application of the cube root function. \return A new AlgebraicVector. */ - using std::log; + using std::cbrt; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::sin() const{ -/*! Compute the sine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::icbrt() const { +/*! Element-wise application of the inverse cube root function. \return A new AlgebraicVector. */ - using std::sin; + using DACE::icbrt; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::cos() const{ -/*! Compute the cosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::hypot(const AlgebraicVector &obj) const { +/*! Element-wise application of the hypotenuse function hypot(x,y). + \param[in] obj AlgebraicVector second argument. \return A new AlgebraicVector. - */ - using std::cos; +*/ + using std::hypot; const size_t size = this->size(); + if(obj.size() != size) + throw std::runtime_error("DACE::AlgebraicVector::hypot(): Vectors must have the same length."); + AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::tan() const{ -/*! Compute the tangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::exp() const { +/*! Element-wise application of the exponential function. \return A new AlgebraicVector. */ - using std::tan; + using std::exp; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::asin() const{ -/*! Compute the arcsine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::log() const { +/*! Element-wise application of the natural logarithm function. \return A new AlgebraicVector. */ - using std::asin; + using std::log; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::acos() const{ -/*! Compute the arccosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::logb(const double b) const { +/*! Element-wise application of the logarithm function relative to given base. + \param[in] b base for the logarithm \return A new AlgebraicVector. */ - using std::acos; + using DACE::logb; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::atan() const{ -/*! Compute the arctangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::log10() const { +/*! Element-wise application of the decadic logarithm function. \return A new AlgebraicVector. */ - using std::atan; + using std::log10; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::atan2(const AlgebraicVector &obj) const{ -/*! Compute the four-quadrant arctangent of Y/X. Y is the current vector, - whereas X is the AlgebraicVector in input. - The result is copied in a new AlgebraicVector. - \param[in] obj AlgebraicVector - \return A new AlgebraicVector containing the result of the operation Y/X in [-pi, pi]. - \throw std::runtime_error -*/ - using std::atan2; +template AlgebraicVector AlgebraicVector::log2() const { +/*! Element-wise application of the binary logarithm function. + \return A new AlgebraicVector. + */ + using std::log2; const size_t size = this->size(); - if(obj.size() != size) - throw std::runtime_error("DACE::AlgebraicVector::atan2(): Vectors must have the same length."); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::sin() const { +/*! Element-wise application of the sine function. + \return A new AlgebraicVector. + */ + using std::sin; + + const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::sinh() const{ -/*! Compute the hyperbolic sine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::cos() const { +/*! Element-wise application of the cosine function. \return A new AlgebraicVector. */ - using std::sinh; + using std::cos; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::cosh() const{ -/*! Compute the hyperbolic cosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::tan() const { +/*! Element-wise application of the tangent function. \return A new AlgebraicVector. */ - using std::cosh; + using std::tan; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::tanh() const{ -/*! Compute the hyperbolic tangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::asin() const { +/*! Element-wise application of the arcsine function. \return A new AlgebraicVector. */ - using std::tanh; + using std::asin; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::logb(const double b) const{ -/*! Compute the logarithm of a AlgebraicVector with respect to a given base. - The result is copied in a new AlgebraicVector. - \param[in] b base with respect to which the logarithm is computed (default = 10). +template AlgebraicVector AlgebraicVector::acos() const { +/*! Element-wise application of the arccosine function. \return A new AlgebraicVector. */ - using DACE::logb; + using std::acos; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::isrt() const{ -/*! Compute the inverse square root of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::atan() const { +/*! Element-wise application of the arctangent function. \return A new AlgebraicVector. */ - using DACE::isrt; + using std::atan; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::sqr() const{ -/*! Compute the square of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::atan2(const AlgebraicVector &obj) const { +/*! Element-wise application of the four-quadrant arctangent of Y/X. + Y is the current object, whereas X is the argument. + \param[in] obj AlgebraicVector representing X + \return A new AlgebraicVector with elements in [-pi, pi]. + \throw std::runtime_error +*/ + using std::atan2; + + const size_t size = this->size(); + if(obj.size() != size) + throw std::runtime_error("DACE::AlgebraicVector::atan2(): Vectors must have the same length."); + + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector AlgebraicVector::sinh() const { +/*! Element-wise application of the hyperbolic sine function. \return A new AlgebraicVector. */ - using DACE::sqr; + using std::sinh; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::minv() const{ -/*! Compute the multiplicative inverse of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector AlgebraicVector::cosh() const { +/*! Element-wise application of the hyperbolic cosine function. \return A new AlgebraicVector. */ - using DACE::minv; + using std::cosh; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::root(const int p) const{ -/*! Compute the p-th root of a AlgebraicVector. - The result is copied in a new AlgebraicVector. - \param[in] p root to be computed (default = 2). +template AlgebraicVector AlgebraicVector::tanh() const { +/*! Element-wise application of the hyperbolic tangent function. \return A new AlgebraicVector. */ - using DACE::root; + using std::tanh; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::asinh() const{ -/*! Compute the hyperbolic arcsine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. +template AlgebraicVector AlgebraicVector::asinh() const { +/*! Element-wise application of the hyperbolic arcsine function. + \return A new AlgebraicVector. */ - using DACE::asinh; + using std::asinh; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::acosh() const{ -/*! Compute the hyperbolic arccosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. Currently not defined for double. - \return A new AlgebraicVector containing the result of the operation. +template AlgebraicVector AlgebraicVector::acosh() const { +/*! Element-wise application of the hyperbolic arccosine function. + \return A new AlgebraicVector. */ - using DACE::acosh; + using std::acosh; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::atanh() const{ -/*! Compute the hyperbolic arctangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. Currently not defined for double. - \return A new AlgebraicVector containing the result of the operation. +template AlgebraicVector AlgebraicVector::atanh() const { +/*! Element-wise application of the hyperbolic arctangent function. + \return A new AlgebraicVector. */ - using DACE::atanh; + using std::atanh; const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::atanh() const{ /*********************************************************************************** * Vector norm routines ************************************************************************************/ -template T AlgebraicVector::vnorm() const{ -/*! Compute the Euclidean vector norm (length) for a AlgebraicVector. - \return A scalar value corresponding to the result of the operation. +template T AlgebraicVector::length() const { +/*! Compute the length (Euclidean vector norm). + \return Length of the vector. */ using std::sqrt; using DACE::sqr; // Implementational note: these using statements are very subtle and absolutely needed. // They force the compiler to perform argument dependent lookup (ADL) which then finds - // the correct root() and pow() functions even if they are not in DACE:: or std::! + // the correct sqrt() and sqr() functions even if they are not in DACE:: or std::! const size_t size = this->size(); T norm = 0.0; - for(size_t i=0; i AlgebraicVector AlgebraicVector::normalize() const{ +template AlgebraicVector AlgebraicVector::normalize() const { /*! Normalize the vector. \return An AlgebraicVector of unit length. */ - const size_t size = this->size(); - AlgebraicVector temp(size); - T norm = 1.0/this->vnorm(); - - for(size_t i=0; i temp(*this); + temp *= minv(this->length()); return temp; } /*********************************************************************************** * Polynomial evaluation routines ************************************************************************************/ -template<> template V AlgebraicVector::eval(const V &args) const{ +template<> template V AlgebraicVector::eval(const V &args) const { /*! Evaluate a vector of polynomials with any vector type V with arguments and return a vector of results of the same type V. \param[in] args vector (e.g. AlgebraicVector<>) of arguments @@ -867,7 +982,7 @@ template<> template V AlgebraicVector::eval(const V &args) const return compiledDA(*this).eval(args); } -template<> template AlgebraicVector AlgebraicVector::eval(const std::initializer_list l) const{ +template<> template AlgebraicVector AlgebraicVector::eval(const std::initializer_list l) const { /*! Evaluate a vector of polynomials with an braced initializer list of type U and return an AlgebraicVector of type U with the results. \param[in] l Braced initializer list containing the arguments. @@ -879,7 +994,7 @@ template<> template AlgebraicVector AlgebraicVector::eval(con return compiledDA(*this).eval(l); } -template<> template AlgebraicVector AlgebraicVector::evalScalar(const U &arg) const{ +template<> template AlgebraicVector AlgebraicVector::evalScalar(const U &arg) const { /*! Evaluate a vector of polynomials with a single arithmetic type U argument. \param[in] arg single variable of arithmetic type T of the first independent DA variable. \return The result of the evaluation. @@ -898,7 +1013,7 @@ template<> template AlgebraicVector AlgebraicVector::evalScal /*********************************************************************************** * Input/Output routines ************************************************************************************/ -template std::ostream& operator<<(std::ostream &out, const AlgebraicVector &obj){ +template std::ostream& operator<<(std::ostream &out, const AlgebraicVector &obj) { /*! Output a vector to a C++ output stream. \param[in] out standard output stream. \param[in] obj AlgebraicVector to be written to the stream @@ -907,15 +1022,15 @@ template std::ostream& operator<<(std::ostream &out, const Algebraic const size_t size = obj.size(); out << "[[[ " << size << " vector" << std::endl; - for(size_t i=0; i std::istream& operator>>(std::istream &in, AlgebraicVector &obj){ -/*! Input a vector from a C++ input stream. +template std::istream& operator>>(std::istream &in, AlgebraicVector &obj) { +/*! Read a vector from a C++ input stream. \param[in] in standard input stream. \param[in] obj AlgebraicVector to be read from the stream \return Reference to input stream in. @@ -925,7 +1040,7 @@ template std::istream& operator>>(std::istream &in, AlgebraicVector< // try to read the first line getline(in, init_line); - if(in.good()){ + if(in.good()) { // retrieve the size of the vector to be read std::size_t found = init_line.find_first_of(' '); std::string size_str(init_line,4,found-4); @@ -950,8 +1065,8 @@ template std::istream& operator>>(std::istream &in, AlgebraicVector< return in; } -template std::string AlgebraicVector::toString() const{ -/*! Convert the current AlgebraicVector to string. +template std::string AlgebraicVector::toString() const { +/*! Convert to string. \return A string. */ std::ostringstream strs; @@ -963,266 +1078,335 @@ template std::string AlgebraicVector::toString() const{ /*********************************************************************************** * Non-member functions ************************************************************************************/ -template AlgebraicVector cons(const AlgebraicVector &obj){ -/*! Return the constant part of a AlgebraicVector. +template AlgebraicVector cons(const AlgebraicVector &obj) { +/*! Return the constant parts. \return An AlgebraicVector. \sa AlgebraicVector::cons */ return obj.cons(); } -template AlgebraicVector pow(const AlgebraicVector &obj, const int p){ -/*! Elevate a AlgebraicVector to a given integer power. - The result is copied in a new AlgebraicVector. +template AlgebraicVector absolute(const AlgebraicVector &obj) { +/*! Element-wise application of the absolute value function. + \param[in] obj AlgebraicVector. + \return A new AlgebraicVector. + \sa AlgebraicVector::absolute + */ + return obj.absolute(); +} + +template AlgebraicVector trunc(const AlgebraicVector &obj) { +/*! Element-wise application of the truncation function. + \param[in] obj AlgebraicVector. + \return A new AlgebraicVector. + \sa AlgebraicVector::trunc + */ + return obj.trunc(); +} + +template AlgebraicVector round(const AlgebraicVector &obj) { +/*! Element-wise application of the round function. \param[in] obj AlgebraicVector. - \param[in] p power at which the AlgebraicVector is elevated. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. + \sa AlgebraicVector::round + */ + return obj.round(); +} + +template AlgebraicVector mod(const AlgebraicVector &obj) { +/*! Element-wise application of the mod function. + \param[in] obj AlgebraicVector. + \return A new AlgebraicVector. + \sa AlgebraicVector::mod + */ + return obj.mod(); +} + +template AlgebraicVector pow(const AlgebraicVector &obj, const int p) { +/*! Element-wise application of the integer power function. + \param[in] obj AlgebraicVector. + \param[in] p power to raise each element to. + \return A new AlgebraicVector. \sa AlgebraicVector::pow */ return obj.pow(p); } -template AlgebraicVector root(const AlgebraicVector &obj, const int p){ -/*! Compute the p-th root of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector pow(const AlgebraicVector &obj, const double p) { +/*! Element-wise application of the double power function. \param[in] obj AlgebraicVector. - \param[in] p root to be computed (default = 2). - \return A new AlgebraicVector containing the result of the operation. + \param[in] p power to raise each element to. + \return A new AlgebraicVector. + \sa AlgebraicVector::pow + */ + return obj.pow(p); +} + +template AlgebraicVector root(const AlgebraicVector &obj, const int p) { +/*! Element-wise application of the root function. + \param[in] obj AlgebraicVector. + \param[in] p root to be computed. + \return A new AlgebraicVector. \sa AlgebraicVector::root */ return obj.root(p); } -template AlgebraicVector minv(const AlgebraicVector &obj){ -/*! Compute the multiplicative inverse of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector minv(const AlgebraicVector &obj) { +/*! Element-wise application of the multiplicative inverse function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::minv */ return obj.minv(); } -template AlgebraicVector sqr(const AlgebraicVector &obj){ -/*! Compute the square of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector sqr(const AlgebraicVector &obj) { +/*! Element-wise application of the square function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::sqr */ return obj.sqr(); } -template AlgebraicVector sqrt(const AlgebraicVector &obj){ -/*! Compute the square root of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector sqrt(const AlgebraicVector &obj) { +/*! Element-wise application of the square root function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::sqrt */ return obj.sqrt(); } -template AlgebraicVector isrt(const AlgebraicVector &obj){ -/*! Compute the inverse square root of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector isrt(const AlgebraicVector &obj) { +/*! Element-wise application of the inverse square root function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::isrt */ return obj.isrt(); } -template AlgebraicVector exp(const AlgebraicVector &obj){ -/*! Compute the exponential of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector cbrt(const AlgebraicVector &obj) { +/*! Element-wise application of the cube root function. + \param[in] obj AlgebraicVector. + \return A new AlgebraicVector. + \sa AlgebraicVector::cbrt + */ + return obj.cbrt(); +} + +template AlgebraicVector icbrt(const AlgebraicVector &obj) { +/*! Element-wise application of the inverse cube root function. + \param[in] obj AlgebraicVector. + \return A new AlgebraicVector. + \sa AlgebraicVector::icbrt + */ + return obj.icbrt(); +} + +template AlgebraicVector hypot(const AlgebraicVector &X, const AlgebraicVector &Y) { +/*! Element-wise application of the hypotenuse function. + \param[in] X AlgebraicVector containing X. + \param[in] Y AlgebraicVector containing Y. + \return A new AlgebraicVector. + \sa AlgebraicVector::hypot + */ + return X.hypot(Y); +} + +template AlgebraicVector exp(const AlgebraicVector &obj) { +/*! Element-wise application of the exponential function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::exp */ return obj.exp(); } -template AlgebraicVector log(const AlgebraicVector &obj){ -/*! Compute the natural logarithm of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector log(const AlgebraicVector &obj) { +/*! Element-wise application of the natural logarithm function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::log */ return obj.log(); } -template AlgebraicVector logb(const AlgebraicVector &obj, const double b){ -/*! Compute the logarithm of a AlgebraicVector with respect to a given base. - The result is copied in a new AlgebraicVector. +template AlgebraicVector logb(const AlgebraicVector &obj, const double b) { +/*! Element-wise application of the logarithm function relative to given base. \param[in] obj AlgebraicVector. - \param[in] b base with respect to which the logarithm is computed (default = 10). - \return A new AlgebraicVector containing the result of the operation. + \param[in] b base for the logarithm + \return A new AlgebraicVector. \sa AlgebraicVector::logb */ return obj.logb(b); } -template AlgebraicVector sin(const AlgebraicVector &obj){ -/*! Compute the sine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector log10(const AlgebraicVector &obj) { +/*! Element-wise application of the decadic logarithm function. + \param[in] obj AlgebraicVector. + \return A new AlgebraicVector. + \sa AlgebraicVector::log10 + */ + return obj.log10(); +} + +template AlgebraicVector log2(const AlgebraicVector &obj) { +/*! Element-wise application of the binary logarithm function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. + \sa AlgebraicVector::log2 + */ + return obj.log2(); +} + +template AlgebraicVector sin(const AlgebraicVector &obj) { +/*! Element-wise application of the sine function. + \param[in] obj AlgebraicVector. + \return A new AlgebraicVector. \sa AlgebraicVector::sin */ return obj.sin(); } -template AlgebraicVector cos(const AlgebraicVector &obj){ -/*! Compute the cosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector cos(const AlgebraicVector &obj) { +/*! Element-wise application of the cosine function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::cos */ return obj.cos(); } -template AlgebraicVector tan(const AlgebraicVector &obj){ -/*! Compute the tangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector tan(const AlgebraicVector &obj) { +/*! Element-wise application of the tangent function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::tan */ return obj.tan(); } -template AlgebraicVector asin(const AlgebraicVector &obj){ -/*! Compute the arcsine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector asin(const AlgebraicVector &obj) { +/*! Element-wise application of the arcsine function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::asin */ return obj.asin(); } -template AlgebraicVector acos(const AlgebraicVector &obj){ -/*! Compute the arccosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector acos(const AlgebraicVector &obj) { +/*! Element-wise application of the arccosine function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::acos */ return obj.acos(); } -template AlgebraicVector atan(const AlgebraicVector &obj){ -/*! Compute the arctangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector atan(const AlgebraicVector &obj) { +/*! Element-wise application of the arctangent function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::atan */ return obj.atan(); } -template AlgebraicVector atan2(const AlgebraicVector &obj1, const AlgebraicVector &obj2){ -/*! Compute the four-quadrant tangent of obj1/obj2. - The result is copied in a new AlgebraicVector. - \param[in] obj1 AlgebraicVector. - \param[in] obj2 AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation in [-pi, pi]. +template AlgebraicVector atan2(const AlgebraicVector &Y, const AlgebraicVector &X) { +/*! Element-wise application of the four-quadrant arctangent of Y/X. + \param[in] Y AlgebraicVector containing Y. + \param[in] X AlgebraicVector containing X. + \return A new AlgebraicVector with elements in [-pi, pi]. \sa AlgebraicVector::atan2 */ - return obj1.atan(obj2); + return Y.atan2(X); } -template AlgebraicVector sinh(const AlgebraicVector &obj){ -/*! Compute the hyperbolic sine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector sinh(const AlgebraicVector &obj) { +/*! Element-wise application of the hyperbolic sine function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::sinh */ return obj.sinh(); } -template AlgebraicVector cosh(const AlgebraicVector &obj){ -/*! Compute the hyperbolic cosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector cosh(const AlgebraicVector &obj) { +/*! Element-wise application of the hyperbolic cosine function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::cosh */ return obj.cosh(); } -template AlgebraicVector tanh(const AlgebraicVector &obj){ -/*! Compute the hyperbolic tangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector tanh(const AlgebraicVector &obj) { +/*! Element-wise application of the hyperbolic tangent function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::tanh */ return obj.tanh(); } -template AlgebraicVector asinh(const AlgebraicVector &obj){ -/*! Compute the hyperbolic arcsine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector asinh(const AlgebraicVector &obj) { +/*! Element-wise application of the hyperbolic arcsine function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::asinh */ return obj.asinh(); } -template AlgebraicVector acosh(const AlgebraicVector &obj){ -/*! Compute the hyperbolic arccosine of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector acosh(const AlgebraicVector &obj) { +/*! Element-wise application of the hyperbolic arccosine function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::acosh */ return obj.acosh(); } -template AlgebraicVector atanh(const AlgebraicVector &obj){ -/*! Compute the hyperbolic arctangent of a AlgebraicVector. - The result is copied in a new AlgebraicVector. +template AlgebraicVector atanh(const AlgebraicVector &obj) { +/*! Element-wise application of the hyperbolic arctangent function. \param[in] obj AlgebraicVector. - \return A new AlgebraicVector containing the result of the operation. + \return A new AlgebraicVector. \sa AlgebraicVector::atanh */ return obj.atanh(); } -template typename PromotionTrait::returnType dot(const AlgebraicVector &obj1, const AlgebraicVector &obj2){ +template typename PromotionTrait::returnType dot(const AlgebraicVector &obj1, const AlgebraicVector &obj2) { /*! Compute the dot product between two AlgebraicVectors. \param[in] obj1 a AlgebraicVector. \param[in] obj2 a AlgebraicVector. - \return A scalar value, containing the result of the operation. + \return A scalar value,. */ return obj1.dot(obj2); } -template AlgebraicVector::returnType> cross(const AlgebraicVector &obj1, const AlgebraicVector &obj2){ +template AlgebraicVector::returnType> cross(const AlgebraicVector &obj1, const AlgebraicVector &obj2) { /*! Compute the cross product between two 3D AlgebraicVectors. \param[in] obj1 a AlgebraicVector. \param[in] obj2 a AlgebraicVector. - \return A new AlgebraicVector, containing the result of the operation. + \return A new AlgebraicVector,. */ return obj1.cross(obj2); } -template T vnorm(const AlgebraicVector &obj){ -/*! Compute the Euclidean vector norm (length) of an AlgebraicVector. +template T length(const AlgebraicVector &obj) { +/*! Compute the length (Euclidean vector norm). \param[in] obj AlgebraicVector. - \return A scalar value containing the result of the operation. - \sa AlgebraicVector::norm + \return Length of the vector. */ - return obj.vnorm(); + return obj.length(); } -template AlgebraicVector normalize(const AlgebraicVector &obj){ +template AlgebraicVector normalize(const AlgebraicVector &obj) { /*! Normalize an AlgebraicVector. \param[in] obj An AlgebraicVector to normalize. \return An AlgebraicVector of unit length. @@ -1231,7 +1415,7 @@ template AlgebraicVector normalize(const AlgebraicVector &obj) return obj.normalize(); } -template V eval(const AlgebraicVector &obj, const V &args){ +template V eval(const AlgebraicVector &obj, const V &args) { /*! Evaluate an AlgebraicVector with a vector type V of arguments and return a vector of type V with the results. \param[in] obj An AlgebraicVector. @@ -1242,7 +1426,7 @@ template V eval(const AlgebraicVector &obj, const V &args){ return obj.eval(args); } -template AlgebraicVector eval(const AlgebraicVector &obj, const std::initializer_list l){ +template AlgebraicVector eval(const AlgebraicVector &obj, const std::initializer_list l) { /*! Evaluate an AlgebraicVector with an braced initializer list of type T and return an AlgebraicVector of type T with the results. \param[in] obj An AlgebraicVector. @@ -1256,7 +1440,7 @@ template AlgebraicVector eval(const AlgebraicVector &obj, con return obj.eval(l); } -template AlgebraicVector evalScalar(const AlgebraicVector &obj, const U &arg){ +template AlgebraicVector evalScalar(const AlgebraicVector &obj, const U &arg) { /*! Evaluate an AlgebraicVector with a single scalar argument of type U and return an AlgebraicVector containing the results. \param[in] obj The AlgebraicVector to evaluate. diff --git a/interfaces/cxx/include/dace/DA.h b/interfaces/cxx/include/dace/DA.h index d511c54..3fc243a 100644 --- a/interfaces/cxx/include/dace/DA.h +++ b/interfaces/cxx/include/dace/DA.h @@ -65,6 +65,8 @@ class DACE_API DA DACEDA m_index; //!< Index to the DA vector public: + typedef double value_type; //!< Underlying type of DA coefficients (for compatibility with C++ std lib, boost, etc) + /******************************************************************************** * DACE Setup *********************************************************************************/ @@ -99,15 +101,15 @@ class DACE_API DA /******************************************************************************** * Coefficient access and extraction routines *********************************************************************************/ - int isnan() const; - int isinf() const; + int isnan() const; //!< Check if any coefficients are NaNs + int isinf() const; //!< Check if any coefficients are Inf double cons() const; //!< Get constant part of a DA AlgebraicVector linear() const; //!< Get linear part of a DA AlgebraicVector gradient() const; //!< Gradient vector with respect to all independent DA variables double getCoefficient(const std::vector &jj) const; //!< Get specific coefficient void setCoefficient(const std::vector &jj, const double coeff); //!< Set specific coefficient Monomial getMonomial(const unsigned int npos) const; //!< Get the Monomial at given position - void getMonomial(const unsigned int npos, Monomial &m) const; //!< Extract the Monomial at given position + void getMonomial(const unsigned int npos, Monomial &m) const; //!< Extract the Monomial at given position std::vector getMonomials() const; //!< Get std::vector of all non-zero Monomials /******************************************************************************** @@ -208,7 +210,6 @@ class DACE_API DA * Norm and estimation routines *********************************************************************************/ unsigned int size() const; //!< Number of non-zero coefficients - double abs() const; //!< Maximum absolute value of all coefficients double norm(const unsigned int type = 0) const; //!< Different types of norms over all coefficients std::vector orderNorm(const unsigned int var = 0, const unsigned int type = 0) const; //!< Different types of norms over coefficients of each order separately @@ -320,7 +321,6 @@ DACE_API DA LogGammaFunction(const DA &da); DACE_API DA PsiFunction(const unsigned int n, const DA &da); DACE_API unsigned int size(const DA &da); -DACE_API double abs(const DA &da); DACE_API double norm(const DA &da, unsigned int type = 0); DACE_API std::vector orderNorm(const DA &da, unsigned int var = 0, unsigned int type = 0); DACE_API std::vector estimNorm(const DA &da, unsigned int var = 0, unsigned int type = 0, unsigned int nc = DA::getMaxOrder()); @@ -340,6 +340,17 @@ DACE_API DA translateVariable(const DA &da, const unsigned int var = 0, const do DACE_API std::string toString(const DA &da); DACE_API void write(const DA &da, std::ostream &os); +namespace abs_cons { + double abs(const DA &da); //!< Absolute value of constant part. +} + +namespace abs_max { + double abs(const DA &da); //!< Largest coefficient in absolute value. +} + +namespace abs_sum { + double abs(const DA &da); //!< Sum of absolute values of all coefficients. +} /*! Stored DA class representing a DA vector in a binary, setup independent format. */ @@ -349,7 +360,7 @@ class DACE_API storedDA : std::vector static const unsigned int headerSize; public: - storedDA(const DA &da); //!< Constructor from a DA. + storedDA(const DA &da); //!< Constructor from a DA. storedDA(const std::vector &data); //!< Constructor from binary data. storedDA(std::istream &is); //!< Constructor from stream. diff --git a/interfaces/cxx/include/dace/MathExtension.h b/interfaces/cxx/include/dace/MathExtension.h index f2bd00a..4dea2f8 100644 --- a/interfaces/cxx/include/dace/MathExtension.h +++ b/interfaces/cxx/include/dace/MathExtension.h @@ -35,6 +35,7 @@ DACE_API double absolute(const double x); //!< Absolute v DACE_API double cons(const double x); //!< Constant part (i.e. the value x) DACE_API double logb(const double x, const double b = 10.0); //!< Logarithm relative to base b DACE_API double isrt(const double x); //!< Inverse square root +DACE_API double icbrt(const double x); //!< Inverse cube root DACE_API double sqr(const double x); //!< Square DACE_API double minv(const double x); //!< Multiplicative inverse DACE_API double root(const double x, const int p = 2); //!< p-th root From 6530704adcb6d7f63d1a848a6b9632a3d5d32afa Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Fri, 25 Apr 2025 17:17:23 +0100 Subject: [PATCH 5/7] boost::odeint compatibility header --- .../cxx/include/dace/compat_boost_odeint.h | 74 +++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 interfaces/cxx/include/dace/compat_boost_odeint.h diff --git a/interfaces/cxx/include/dace/compat_boost_odeint.h b/interfaces/cxx/include/dace/compat_boost_odeint.h new file mode 100644 index 0000000..1af189c --- /dev/null +++ b/interfaces/cxx/include/dace/compat_boost_odeint.h @@ -0,0 +1,74 @@ +/****************************************************************************** +* * +* DIFFERENTIAL ALGEBRA CORE ENGINE * +* * +******************************************************************************* +* * +* Copyright 2025 Alexander Wittig * +* Licensed under the Apache License, Version 2.0 (the "License"); * +* you may not use this file except in compliance with the License. * +* You may obtain a copy of the License at * +* * +* http://www.apache.org/licenses/LICENSE-2.0 * +* * +* Unless required by applicable law or agreed to in writing, software * +* distributed under the License is distributed on an "AS IS" BASIS, * +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * +* See the License for the specific language governing permissions and * +* limitations under the License. * +* * +*******************************************************************************/ + +/* + * compat_boost_odeint.h + * + * Created on: Apr 25, 2025 + * Author: Alexander Wittig + */ + +#ifndef DINAMICA_COMPAT_BOOST_ODEINT_H_ +#define DINAMICA_COMPAT_BOOST_ODEINT_H_ + +/*! This file needs to be included after boost/numeric/odeint.hpp and dace/dace.h to provide + a compatibility shim to allow odeint to work with DACE::AlgebraicVector<> as a state type. + It can be used both with DA and double data types. + + Additionally, after including this file, an implementation of abs(DACE::DA) must be provided + within the DACE namespace. + For most practical uses, this can be done by selecting one of the predefined implementations + provided in DACE::abs_cons, DACE::abs_max, or DACE::abs_sum. + + Example: + \code + #include + #include + #include + + // select implementation of abs(DA) + namespace DACE { using DACE::abs_max::abs; } + + + using namespace boost::numeric::odeint; + using namespace DACE; + + // define RHS, call DA::init, ... + + typedef state_type AlgebraicVector; + + state_type x(3); + x[0] = 1.0; x[1] = 2.0; x[3] = 3.0; + x += 0.01 * x.identity(3); + integrate_adaptive( make_controlled( 1e-8, 1e-8, runge_kutta_dopri5() ), RHS, x, 0.0, 10.0, 0.1 ); + \endcode + */ + + namespace boost { namespace numeric { namespace odeint { + // mark AlgebraicVectors as resizable (using the standard container interface inherited from std::vector) + template struct is_resizeable> + { + typedef boost::true_type type; + static const bool value = type::value; + }; +} } } + +#endif /* DINAMICA_COMPAT_BOOST_ODEINT_H_ */ From 0a922d9ff3dd93609d5e857cce9a834108732155 Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Fri, 25 Apr 2025 18:17:31 +0100 Subject: [PATCH 6/7] Add AlgebraicVector::norm --- interfaces/cxx/AlgebraicVector.cpp | 59 +++++----- interfaces/cxx/MathExtension.cpp | 8 ++ interfaces/cxx/include/dace/AlgebraicVector.h | 14 ++- .../cxx/include/dace/AlgebraicVector_t.h | 103 ++++++++++++------ interfaces/cxx/include/dace/MathExtension.h | 1 + 5 files changed, 114 insertions(+), 71 deletions(-) diff --git a/interfaces/cxx/AlgebraicVector.cpp b/interfaces/cxx/AlgebraicVector.cpp index 3fbf82f..6b21eeb 100644 --- a/interfaces/cxx/AlgebraicVector.cpp +++ b/interfaces/cxx/AlgebraicVector.cpp @@ -43,7 +43,7 @@ namespace DACE { * Coefficient access routines ************************************************************************************/ #ifdef WITH_ALGEBRAICMATRIX -template<> AlgebraicMatrix AlgebraicVector::linear() const{ +template<> AlgebraicMatrix AlgebraicVector::linear() const { /*! Return the linear part of a AlgebraicVector. NOT DEFINED FOR TYPES OTHER THAN DA. \return A AlgebraicMatrix of dimension size by nvar, where size is the size of the AlgebraicVector considered and nvar is the number of variables defined @@ -57,13 +57,13 @@ template<> AlgebraicMatrix AlgebraicVector::linear() const{ const int nvar = DA::getMaxVariables(); AlgebraicMatrix out(size, nvar); - for(size_t i=0; i std::vector< std::vector > AlgebraicVector::linear() const{ +template<> std::vector< std::vector > AlgebraicVector::linear() const { /*! Return the linear part of a AlgebraicVector. NOT DEFINED FOR TYPES OTHER THAN DA. \return A std::vector< std::vector >, where each std::vector contains the linear part of the corresponding DA included in the original AlgebraicVector. @@ -74,15 +74,14 @@ template<> std::vector< std::vector > AlgebraicVector::linear() cons const size_t size = this->size(); std::vector< std::vector > out(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::trim(const unsigned int min, const unsigned int max) const -{ +template<> AlgebraicVector AlgebraicVector::trim(const unsigned int min, const unsigned int max) const { /*! Returns an AlgebraicVector with all monomials of order less than min and greater than max removed (trimmed). The result is copied in a new AlgebraicVector. \param[in] min minimum order to be preserved. \param[in] max maximum order to be preserved. @@ -102,7 +101,7 @@ template<> AlgebraicVector AlgebraicVector::trim(const unsigned int min, /*********************************************************************************** * Math routines ************************************************************************************/ -template<> AlgebraicVector AlgebraicVector::deriv(const unsigned int p) const{ +template<> AlgebraicVector AlgebraicVector::deriv(const unsigned int p) const { /*! Compute the derivative of a AlgebraicVector with respect to variable p. The result is copied in a new AlgebraicVector. NOT DEFINED FOR TYPES OTHER THAN DA. \param[in] p variable with respect to which the derivative is calculated. @@ -113,13 +112,13 @@ template<> AlgebraicVector AlgebraicVector::deriv(const unsigned int p) */ const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::integ(const unsigned int p) const{ +template<> AlgebraicVector AlgebraicVector::integ(const unsigned int p) const { /*! Compute the integral of a AlgebraicVector with respect to variable p. The result is copied in a new AlgebraicVector. NOT DEFINED FOR TYPES OTHER THAN DA. \param[in] p variable with respect to which the integral is calculated. @@ -130,7 +129,7 @@ template<> AlgebraicVector AlgebraicVector::integ(const unsigned int p) */ const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::integ(const unsigned int p) /*********************************************************************************** * Polynomial evaluation routines ************************************************************************************/ -template<> compiledDA AlgebraicVector::compile() const{ +template<> compiledDA AlgebraicVector::compile() const { /*! Compile vector of polynomials and create a compiledDA object. \return The compiled DA object. \note This DA specific function is only available in AlgebraicVector. @@ -149,7 +148,7 @@ template<> compiledDA AlgebraicVector::compile() const{ return compiledDA(*this); } -template<> AlgebraicVector AlgebraicVector::plug(const unsigned int var, const double val) const{ +template<> AlgebraicVector AlgebraicVector::plug(const unsigned int var, const double val) const { /*! Partial evaluation of vector of polynomials. In each element of the vector, variable var is replaced by the value val. The resulting vector of DAs is returned. @@ -162,7 +161,7 @@ template<> AlgebraicVector AlgebraicVector::plug(const unsigned int var, */ const size_t size = this->size(); AlgebraicVector temp(size); - for(size_t i=0; i AlgebraicVector AlgebraicVector::plug(const unsigned int var, This is NOT intended for public use. Limited error checking is performed in accordance with the exclusive use of this routine in map inversion. */ -template<> void AlgebraicVector::matrix_inverse(std::vector< std::vector > &A){ +template<> void AlgebraicVector::matrix_inverse(std::vector< std::vector > &A) { using std::abs; const size_t n = A.size(); std::vector indexc(n), indexr(n), ipiv(n, 0); - for (size_t i=0; i= big){ + if (abs(A[j][k]) >= big) { big = abs(A[j][k]); irow = j; icol = k;} @@ -204,7 +203,7 @@ template<> void AlgebraicVector::matrix_inverse(std::vector< std::vector void AlgebraicVector::matrix_inverse(std::vector< std::vector AlgebraicVector AlgebraicVector::invert() const{ +template<> AlgebraicVector AlgebraicVector::invert() const { /*! Invert the polynomials map given by the AlgebraicVector. \return the inverted polynomials \throw std::runtime_error @@ -258,13 +257,13 @@ template<> AlgebraicVector AlgebraicVector::invert() const{ // Compute DA representation of the inverse of the linear part of the map and its composition with non-linear part AN AlgebraicVector Linv(nvar); // Linv = AI*AN - for(size_t i=0; i AlgebraicVector AlgebraicVector::invert() const{ // Iterate to obtain the inverse map AlgebraicVector MI = Linv; - for(unsigned int i=1; i AlgebraicVector AlgebraicVector::invert() const{ /******************************************************************************** * Static factory routines *********************************************************************************/ -template<> AlgebraicVector AlgebraicVector::identity(const size_t n){ +template<> AlgebraicVector AlgebraicVector::identity(const size_t n) { /*! Return the DA identity of dimension n. \param[in] n The dimendion of the identity. \return AlgebraicVector containing the DA identity in n dimensions @@ -291,7 +290,7 @@ template<> AlgebraicVector AlgebraicVector::identity(const size_t n){ error will be the result. */ AlgebraicVector temp(n); - for(size_t i=0; i < n; i++){ + for(size_t i=0; i < n; i++) { temp[i] = DA((int)(i+1));} return temp; @@ -301,7 +300,7 @@ template<> AlgebraicVector AlgebraicVector::identity(const size_t n){ * Non-member functions ************************************************************************************/ #ifdef WITH_ALGEBRAICMATRIX -template<> AlgebraicMatrix linear(const AlgebraicVector &obj){ +template<> AlgebraicMatrix linear(const AlgebraicVector &obj) { /*! Return the linear part of a AlgebraicVector. \param[in] obj AlgebraicVector to extract linear part from \return An AlgebraicMatrix of dimensions size by nvar, where @@ -313,7 +312,7 @@ template<> AlgebraicMatrix linear(const AlgebraicVector &obj){ \sa AlgebraicVector::linear */ #else -template<> std::vector< std::vector > linear(const AlgebraicVector &obj){ +template<> std::vector< std::vector > linear(const AlgebraicVector &obj) { /*! Return the linear part of a AlgebraicVector. Only defined for AlgebraicVector. \param[in] obj AlgebraicVector to extract linear part from \return A std::vector< std::vector > containing the linear parts of @@ -327,7 +326,7 @@ template<> std::vector< std::vector > linear(const AlgebraicVector & return obj.linear(); } -template<> AlgebraicVector trim(const AlgebraicVector &obj, unsigned int min, unsigned int max){ +template<> AlgebraicVector trim(const AlgebraicVector &obj, unsigned int min, unsigned int max) { /*! Returns an AlgebraicVector with all monomials of order less than min and greater than max removed (trimmed). The result is copied in a new AlgebraicVector. \param[in] obj the AlgebraicVector to be trimmed. \param[in] min minimum order to be preserved. @@ -341,7 +340,7 @@ template<> AlgebraicVector trim(const AlgebraicVector &obj, unsigned int return obj.trim(min, max); } -template<> AlgebraicVector deriv(const AlgebraicVector &obj, const unsigned int p){ +template<> AlgebraicVector deriv(const AlgebraicVector &obj, const unsigned int p) { /*! Compute the derivative of a AlgebraicVector with respect to variable p. The result is copied in a new AlgebraicVector. \param[in] obj AlgebraicVector. @@ -355,7 +354,7 @@ template<> AlgebraicVector deriv(const AlgebraicVector &obj, const unsig return obj.deriv(p); } -template<> AlgebraicVector integ(const AlgebraicVector &obj, const unsigned int p){ +template<> AlgebraicVector integ(const AlgebraicVector &obj, const unsigned int p) { /*! Compute the integral of a AlgebraicVector with respect to variable p. The result is copied in a new AlgebraicVector. \param[in] obj AlgebraicVector. @@ -369,7 +368,7 @@ template<> AlgebraicVector integ(const AlgebraicVector &obj, const unsig return obj.integ(p); } -template<> compiledDA compile(const AlgebraicVector &obj){ +template<> compiledDA compile(const AlgebraicVector &obj) { /*! Compile vector of polynomials and create a compiledDA object. \param[in] obj The AlgebraicVector to compile. \return The compiled DA object. @@ -380,7 +379,7 @@ template<> compiledDA compile(const AlgebraicVector &obj){ return obj.compile(); } -template<> AlgebraicVector plug(const AlgebraicVector &obj, const unsigned int var, const double val){ +template<> AlgebraicVector plug(const AlgebraicVector &obj, const unsigned int var, const double val) { /*! Partial evaluation of vector of polynomials. In each element of the vector, variable var is replaced by the value val. The resulting vector of DAs is returned. diff --git a/interfaces/cxx/MathExtension.cpp b/interfaces/cxx/MathExtension.cpp index e381542..5325763 100644 --- a/interfaces/cxx/MathExtension.cpp +++ b/interfaces/cxx/MathExtension.cpp @@ -93,4 +93,12 @@ double root(const double x, const int p){ return std::pow(x, 1.0/p); } +double norm(const double x, const int p){ +/*! norm of x. + \param[in] x Function argument. + \param[in] type Type of norm (ignored for double). + */ + return std::abs(x); +} + } diff --git a/interfaces/cxx/include/dace/AlgebraicVector.h b/interfaces/cxx/include/dace/AlgebraicVector.h index 546f670..e0b3931 100644 --- a/interfaces/cxx/include/dace/AlgebraicVector.h +++ b/interfaces/cxx/include/dace/AlgebraicVector.h @@ -128,8 +128,9 @@ template class AlgebraicVector : public std::vector //!< Dot product (scalar product, inner product) of two vectors. template AlgebraicVector::returnType> cross(const AlgebraicVector &obj) const; //!< Cross product of two vectors of length 3. - T length() const; //!< Length of the vector in Euclidean norm. + T length() const; //!< Length of the vector in Euclidean norm. AlgebraicVector normalize() const; //!< Normalized vector of unit length along this vector. + // XXX: various Jacobians, gradients, curls, etc? /*********************************************************************************** * Special routines (DA related) @@ -144,10 +145,12 @@ template class AlgebraicVector : public std::vector AlgebraicVector trim(const unsigned int min, const unsigned int max = DA::getMaxOrder()) const; //!< Trim the coefficients of each components to particular orders. DA only. AlgebraicVector invert() const; //!< Inverse function of the AlgebraicVector. DA only. - // XXX: define and add the norm estimation routines from DA including norm(p), convergence radius estimation - // XXX: add jacobian routine (returns a DA matrix containing the jacobian) - /* - double norm(const unsigned int type = 0) const; //!< Different types of norms over all coefficients + + /******************************************************************************** + * DA norm routines + *********************************************************************************/ + AlgebraicVector norm(const unsigned int type = 0) const; //!< Element-wise DA norm + /* XXX: define and add the norm estimation routines from DA including convergence radius estimation std::vector orderNorm(const unsigned int var = 0, const unsigned int type = 0) const; //!< Different types of norms over coefficients of each order separately std::vector estimNorm(const unsigned int var = 0, const unsigned int type = 0, const unsigned int nc = DA::getMaxOrder()) const; @@ -233,6 +236,7 @@ template AlgebraicVector eval(const AlgebraicVector AlgebraicVector evalScalar(const AlgebraicVector &obj, const U &arg); template compiledDA compile(const AlgebraicVector &obj); template AlgebraicVector plug(const AlgebraicVector &obj, const unsigned int var, const double val = 0.0); +template AlgebraicVector norm(const AlgebraicVector &obj, const unsigned int type = 0); // specializations for various DA specific routines implemented and instantiated directly in the library instead of in a template #ifdef WITH_ALGEBRAICMATRIX diff --git a/interfaces/cxx/include/dace/AlgebraicVector_t.h b/interfaces/cxx/include/dace/AlgebraicVector_t.h index 138228c..48ac612 100644 --- a/interfaces/cxx/include/dace/AlgebraicVector_t.h +++ b/interfaces/cxx/include/dace/AlgebraicVector_t.h @@ -438,41 +438,6 @@ template AlgebraicVector: return temp; } -template template typename PromotionTrait::returnType AlgebraicVector::dot(const AlgebraicVector &obj) const { -/*! Compute the dot product with another AlgebraicVector. - \param[in] obj the other AlgebraicVector. - \return A scalar value,. - \throw std::runtime_error - */ - const size_t size = this->size(); - if(size != obj.size()) - throw std::runtime_error("DACE::AlgebraicVector::dot(): Vectors must have the same length."); - - typename PromotionTrait::returnType temp = 0.0; - for(size_t i=0; i template AlgebraicVector::returnType> AlgebraicVector::cross(const AlgebraicVector &obj) const { -/*! Compute the cross product with another 3D AlgebraicVector. - \param[in] obj The other AlgebraicVector. - \return A new AlgebraicVector. - \throw std::runtime_error - */ - if((this->size() != 3) || (obj.size() != 3)) - throw std::runtime_error("DACE::AlgebraicVector::cross(): Inputs must be 3 element AlgebraicVectors."); - - AlgebraicVector::returnType> temp(3); - - temp[0] = ((*this)[1] * obj[2]) - ((*this)[2] * obj[1]); - temp[1] = ((*this)[2] * obj[0]) - ((*this)[0] * obj[2]); - temp[2] = ((*this)[0] * obj[1]) - ((*this)[1] * obj[0]); - - return temp; -} - /*********************************************************************************** * Math routines ************************************************************************************/ @@ -937,8 +902,43 @@ template AlgebraicVector AlgebraicVector::atanh() const { } /*********************************************************************************** -* Vector norm routines +* Vector routines ************************************************************************************/ +template template typename PromotionTrait::returnType AlgebraicVector::dot(const AlgebraicVector &obj) const { +/*! Compute the dot product with another AlgebraicVector. + \param[in] obj the other AlgebraicVector. + \return A scalar value,. + \throw std::runtime_error + */ + const size_t size = this->size(); + if(size != obj.size()) + throw std::runtime_error("DACE::AlgebraicVector::dot(): Vectors must have the same length."); + + typename PromotionTrait::returnType temp = 0.0; + for(size_t i=0; i template AlgebraicVector::returnType> AlgebraicVector::cross(const AlgebraicVector &obj) const { +/*! Compute the cross product with another 3D AlgebraicVector. + \param[in] obj The other AlgebraicVector. + \return A new AlgebraicVector. + \throw std::runtime_error + */ + if((this->size() != 3) || (obj.size() != 3)) + throw std::runtime_error("DACE::AlgebraicVector::cross(): Inputs must be 3 element AlgebraicVectors."); + + AlgebraicVector::returnType> temp(3); + + temp[0] = ((*this)[1] * obj[2]) - ((*this)[2] * obj[1]); + temp[1] = ((*this)[2] * obj[0]) - ((*this)[0] * obj[2]); + temp[2] = ((*this)[0] * obj[1]) - ((*this)[1] * obj[0]); + + return temp; +} + template T AlgebraicVector::length() const { /*! Compute the length (Euclidean vector norm). \return Length of the vector. @@ -1010,6 +1010,25 @@ template<> template AlgebraicVector AlgebraicVector::evalScal return compiledDA(*this).evalScalar(arg); } +/*********************************************************************************** +* DA norm routines +************************************************************************************/ +template AlgebraicVector AlgebraicVector::norm(const unsigned int type) const { +/*! Element-wise application of the norm function. + \param[in] type type of norm to be computed. See DA::norm. + \return A new AlgebraicVector. + \sa DA::norm() + */ + using DACE::norm; + + const size_t size = this->size(); + AlgebraicVector temp(size); + for(size_t i=0; i AlgebraicVector evalScalar(const AlgebraicVector &ob return obj.evalScalar(arg); } +template AlgebraicVector norm(const AlgebraicVector &obj, const unsigned int type) { +/*! Element-wise application of the norm function. + \param[in] obj AlgebraicVector. + \param[in] type type of norm to be computed. See DA::norm. + \return A new AlgebraicVector. + \sa AlgebraicVector::norm + \sa DA::norm + */ + return obj.norm(type); +} + + } #endif /* DINAMICA_ALGEBRAICVECTOR_T_H_ */ diff --git a/interfaces/cxx/include/dace/MathExtension.h b/interfaces/cxx/include/dace/MathExtension.h index 4dea2f8..99829ec 100644 --- a/interfaces/cxx/include/dace/MathExtension.h +++ b/interfaces/cxx/include/dace/MathExtension.h @@ -39,6 +39,7 @@ DACE_API double icbrt(const double x); //!< Inverse cu DACE_API double sqr(const double x); //!< Square DACE_API double minv(const double x); //!< Multiplicative inverse DACE_API double root(const double x, const int p = 2); //!< p-th root +DACE_API double norm(const double x, const int type = 0); //!< norm (same as abs) } From 8f7a843d8acc748e8bba69074a7cc4b8062f15a9 Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Sat, 26 Apr 2025 23:13:58 +0100 Subject: [PATCH 7/7] add compatibility shims for boost::odeint vector_space_algebra --- .../cxx/include/dace/compat_boost_odeint.h | 52 +++++++++++++++++-- 1 file changed, 48 insertions(+), 4 deletions(-) diff --git a/interfaces/cxx/include/dace/compat_boost_odeint.h b/interfaces/cxx/include/dace/compat_boost_odeint.h index 1af189c..387f70d 100644 --- a/interfaces/cxx/include/dace/compat_boost_odeint.h +++ b/interfaces/cxx/include/dace/compat_boost_odeint.h @@ -30,10 +30,13 @@ #define DINAMICA_COMPAT_BOOST_ODEINT_H_ /*! This file needs to be included after boost/numeric/odeint.hpp and dace/dace.h to provide - a compatibility shim to allow odeint to work with DACE::AlgebraicVector<> as a state type. + a compatibility shim to allow odeint to work with DACE::AlgebraicVector as a state type. It can be used both with DA and double data types. - Additionally, after including this file, an implementation of abs(DACE::DA) must be provided + Use of the default range_algebra is recommended, but DACE::AlgebraicVector also works with + the vector_space_algebra. + + Additionally, before including this file, an implementation of abs(DACE::DA) must be provided within the DACE namespace. For most practical uses, this can be done by selecting one of the predefined implementations provided in DACE::abs_cons, DACE::abs_max, or DACE::abs_sum. @@ -42,10 +45,10 @@ \code #include #include - #include - // select implementation of abs(DA) + // select implementation of abs(DA) before compat_boost_odeint.h namespace DACE { using DACE::abs_max::abs; } + #include using namespace boost::numeric::odeint; @@ -62,6 +65,8 @@ \endcode */ +// range_algebra helpers (default) + namespace boost { namespace numeric { namespace odeint { // mark AlgebraicVectors as resizable (using the standard container interface inherited from std::vector) template struct is_resizeable> @@ -71,4 +76,43 @@ }; } } } +// vector_space_algebra helpers + +namespace boost { namespace numeric { namespace odeint { + // specialization to compute double max norm of AlgebraicVector + template struct vector_space_norm_inf> + { + typedef double result_type; + double operator()(const DACE::AlgebraicVector &x) const + { + using DACE::abs; + using std::abs; + + double res = 0.0; + for(unsigned int i = 0; i < x.size(); i++) + { + const double temp = abs(x[i]); + if(temp > res) res = temp; + } + return res; + } + }; +} } } + +namespace DACE { + // AlgebraicVector component-wise abs() with same output type. Required by vector_space_algebra in boost. + template DACE::AlgebraicVector abs(const DACE::AlgebraicVector &x) + { + using DACE::abs; + using std::abs; + + DACE::AlgebraicVector res(x.size()); + for(unsigned int i = 0; i < x.size(); i++) + { + res[i] = abs(x[i]); + } + return res; + } +} + #endif /* DINAMICA_COMPAT_BOOST_ODEINT_H_ */