From 2922c85a002fc3439f54c1d3d21ed4e80eeee6fb Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Mon, 21 Apr 2025 17:41:00 +0100 Subject: [PATCH 1/4] Add some LIKELY keywords to basic loops, whitespace --- core/dacemath.c | 66 ++++++++++++++++++++++++------------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/core/dacemath.c b/core/dacemath.c index a56c89f..195f29d 100755 --- a/core/dacemath.c +++ b/core/dacemath.c @@ -222,14 +222,14 @@ void daceMultiplyMonomials(const DACEDA *ina, const DACEDA *inb, DACEDA *inc) monomial *ib = ipob, *ic = ipoc; monomial *const ibmax = ipob + ilmb, *const icmax = ipoc + ilmc; - for (monomial *i = ipoa; i < ipoa + illa; i++) + for(monomial *i = ipoa; i < ipoa + illa; i++) { - while (ib->ii < i->ii && ib < ibmax) + while(ib->ii < i->ii && ib < ibmax) ib++; - if (ib == ibmax) break; - if (ib->ii == i->ii) + if(ib == ibmax) break; + if(ib->ii == i->ii) { - if (ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -324,7 +324,7 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb) { for(monomial *ia = ipoa; ia < ipoa+illa; ia++) { - if(DACECom.ieo[ia->ii] > DACECom_t.nocut) + if(UNLIKELY(DACECom.ieo[ia->ii] > DACECom_t.nocut)) continue; const double c = ia->cc*ckon; @@ -341,13 +341,13 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb) monomial *const ibmax = ipob+ilmb; for(monomial *ia = ipoa; ia < ipoa+illa; ia++) { - if(DACECom.ieo[ia->ii] > DACECom_t.nocut) + if(UNLIKELY(DACECom.ieo[ia->ii] > DACECom_t.nocut)) continue; const double c = ia->cc*ckon; if(!(fabs(c) <= DACECom_t.eps)) { - if(ib >= ibmax) + if(UNLIKELY(ib >= ibmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -451,13 +451,13 @@ void daceDivideByVariable(const DACEDA *ina, const unsigned int var, const unsig const unsigned int ic1 = DACECom.ie1[i->ii]; const unsigned int ic2 = DACECom.ie2[i->ii]; const unsigned int ipow = (ic2/idiv)%ibase; - if(ipow < p) + if(UNLIKELY(ipow < p)) { daceSetError(__func__, DACE_ERROR, 42); daceCreateConstant(inc, 0.0); return; } - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -474,13 +474,13 @@ void daceDivideByVariable(const DACEDA *ina, const unsigned int var, const unsig const unsigned int ic1 = DACECom.ie1[i->ii]; const unsigned int ic2 = DACECom.ie2[i->ii]; const unsigned int ipow = (ic1/idiv)%ibase; - if(ipow < p) + if(UNLIKELY(ipow < p)) { daceSetError(__func__, DACE_ERROR, 42); daceCreateConstant(inc, 0.0); return; } - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -534,7 +534,7 @@ void daceDifferentiate(const unsigned int idif, const DACEDA *ina, DACEDA *inc) const unsigned int ipow = (ic2/idiv)%ibase; if(ipow == 0 || DACECom.ieo[i->ii] > DACECom_t.nocut+1) continue; - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -553,7 +553,7 @@ void daceDifferentiate(const unsigned int idif, const DACEDA *ina, DACEDA *inc) const unsigned int ipow = (ic1/idiv)%ibase; if(ipow == 0 || DACECom.ieo[i->ii] > DACECom_t.nocut+1) continue; - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -602,15 +602,15 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc) { for(monomial *i = ipoa; i < ipoa+illa; i++) { - if(DACECom.ieo[i->ii] >= DACECom_t.nocut) + if(UNLIKELY(DACECom.ieo[i->ii] >= DACECom_t.nocut)) continue; const unsigned int ic1 = DACECom.ie1[i->ii]; const unsigned int ic2 = DACECom.ie2[i->ii]; const unsigned int ipow = (ic2/idiv)%ibase; const double ccc = i->cc/(ipow+1); - if(!(fabs(ccc) <= DACECom_t.eps)) + if(LIKELY(!(fabs(ccc) <= DACECom_t.eps))) { - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -625,15 +625,15 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc) { for(monomial *i = ipoa; i < ipoa+illa; i++) { - if(DACECom.ieo[i->ii] >= DACECom_t.nocut) + if(UNLIKELY(DACECom.ieo[i->ii] >= DACECom_t.nocut)) continue; const unsigned int ic1 = DACECom.ie1[i->ii]; const unsigned int ic2 = DACECom.ie2[i->ii]; const unsigned int ipow = (ic1/idiv)%ibase; const double ccc = i->cc/(ipow+1); - if(!(fabs(ccc) <= DACECom_t.eps)) + if(LIKELY(!(fabs(ccc) <= DACECom_t.eps))) { - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); break; @@ -2047,7 +2047,7 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co monomial *ia = ipoa, *ib = ipob, *ic = ipoc; monomial *const iamax = ipoa+illa, *const ibmax = ipob+illb, *const icmax = ipoc+ilmc; - if(illa > 0 && illb > 0) + if(LIKELY(illa > 0 && illb > 0)) { // both polynomials have coefficients, merge until one runs out unsigned int ja = ia->ii; @@ -2057,12 +2057,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co if(ja == jb) { // add the two terms - if(DACECom.ieo[ja] <= DACECom_t.nocut) + if(LIKELY(DACECom.ieo[ja] <= DACECom_t.nocut)) { const double ccc = ia->cc*afac + ib->cc*bfac; - if(!(fabs(ccc) <= DACECom_t.eps)) + if(LIKELY(!(fabs(ccc) <= DACECom_t.eps))) { - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); daceSetLength(inc, ilmc); @@ -2081,12 +2081,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co else if(ja < jb) { // store term a - if(DACECom.ieo[ja] <= DACECom_t.nocut) + if(LIKELY(DACECom.ieo[ja] <= DACECom_t.nocut)) { const double ccc = ia->cc*afac; - if(!(fabs(ccc) <= DACECom_t.eps)) + if(LIKELY(!(fabs(ccc) <= DACECom_t.eps))) { - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); daceSetLength(inc, ilmc); @@ -2104,12 +2104,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co else { // store term b - if(DACECom.ieo[jb] <= DACECom_t.nocut) + if(LIKELY(DACECom.ieo[jb] <= DACECom_t.nocut)) { const double ccc = ib->cc*bfac; - if(!(fabs(ccc) <= DACECom_t.eps)) + if(LIKELY(!(fabs(ccc) <= DACECom_t.eps))) { - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); daceSetLength(inc, ilmc); @@ -2145,12 +2145,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co for(monomial *is = ismin; is < ismax; is++) { - if(DACECom.ieo[is->ii] <= DACECom_t.nocut) + if(LIKELY(DACECom.ieo[is->ii] <= DACECom_t.nocut)) { const double ccc = is->cc*fac; - if(!(fabs(ccc) <= DACECom_t.eps)) + if(LIKELY(!(fabs(ccc) <= DACECom_t.eps))) { - if(ic >= icmax) + if(UNLIKELY(ic >= icmax)) { daceSetError(__func__, DACE_ERROR, 21); daceSetLength(inc, ilmc); From 72f015c73d95a36bacd0f9e5d5342b1f6f22d92c Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Mon, 21 Apr 2025 18:17:29 +0100 Subject: [PATCH 2/4] isNan and isInf functions --- core/dacebasic.c | 38 +++++++++++++++++++++++++++++++++++- core/include/dace/dacebase.h | 2 ++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/core/dacebasic.c b/core/dacebasic.c index 25bc75a..905c61c 100644 --- a/core/dacebasic.c +++ b/core/dacebasic.c @@ -569,7 +569,7 @@ void daceTrim(const DACEDA *ina, const unsigned int imin, const unsigned int ima /*! Copy monomials from a DA object ina to DA object inb if the same monomial is non-zero in DA object inc, while filtering out terms below the current - cutoff + cutoff. \param[in] ina Pointer to DA object to filter \param[in] inb Pointer to DA object to store the filtered result in \param[in] inc Pointer to DA object providing the filter template @@ -627,4 +627,40 @@ void daceFilter(const DACEDA *ina, DACEDA *inb, const DACEDA *inc) daceSetLength(inb, ib-ipob); } + +/*! Check each coefficient of DA object ina to see if any of them are NANs (not a number). + \param[in] ina Pointer to DA object to check + \return True (non-zero) if any of the coefficients of ina is NAN +*/ +unsigned int daceIsNan(const DACEDA *ina) +{ + monomial *ipoa; unsigned int ilma, illa; + + daceVariableInformation(ina, &ipoa, &ilma, &illa); + + for(monomial *ia = ipoa; ia < ipoa+illa; ia++) + { + if(isnan(ia->cc)) + return true; + } + return false; +} + +/*! Check each coefficient of DA object ina to see if any of them are INF (infinity). + \param[in] ina Pointer to DA object to check + \return True (non-zero) if any of the coefficients of ina is INF +*/ +unsigned int daceIsInf(const DACEDA *ina) +{ + monomial *ipoa; unsigned int ilma, illa; + + daceVariableInformation(ina, &ipoa, &ilma, &illa); + + for(monomial *ia = ipoa; ia < ipoa+illa; ia++) + { + if(isinf(ia->cc)) + return true; + } + return false; +} /** @}*/ diff --git a/core/include/dace/dacebase.h b/core/include/dace/dacebase.h index 5062d92..ae6c6c3 100644 --- a/core/include/dace/dacebase.h +++ b/core/include/dace/dacebase.h @@ -144,6 +144,8 @@ DACE_API void daceCopy(const DACEDA REF(ina), DACEDA REF(inb)); DACE_API void daceCopyFiltering(const DACEDA REF(ina), DACEDA REF(inb)); DACE_API void daceFilter(const DACEDA REF(ina), DACEDA REF(inb), const DACEDA REF(inc)); DACE_API void daceTrim(const DACEDA REF(ina), const unsigned int imin, const unsigned int imax, DACEDA REF(inc)); +DACE_API unsigned int daceIsNan(const DACEDA REF(ina)); +DACE_API unsigned int daceIsInf(const DACEDA REF(ina)); /******************************************************************************** * Basic DACE arithmetic operations From 3e50d2cd0dc9919202ec7c9179b8c70929cc1747 Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Mon, 21 Apr 2025 18:17:37 +0100 Subject: [PATCH 3/4] white space --- core/daceeval.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/core/daceeval.c b/core/daceeval.c index 831db19..a885351 100644 --- a/core/daceeval.c +++ b/core/daceeval.c @@ -61,12 +61,12 @@ double daceEvalMonomials(const DACEDA *ina, const DACEDA *inb) monomial *const ibmax = ipob + ilmb; double res = 0.0; - for (monomial *i = ipoa; i < ipoa + illa; i++) + for(monomial *i = ipoa; i < ipoa + illa; i++) { - while (ib->ii < i->ii && ib < ibmax) + while(ib->ii < i->ii && ib < ibmax) ib++; - if (ib == ibmax) break; - if (ib->ii == i->ii) res += ib->cc*i->cc; + if(ib == ibmax) break; + if(ib->ii == i->ii) res += ib->cc*i->cc; } return res; From 570b72e25863243c8eaaa3215f8983e1bd7ad999 Mon Sep 17 00:00:00 2001 From: Alexander Wittig Date: Mon, 21 Apr 2025 18:34:37 +0100 Subject: [PATCH 4/4] expose isnan and isinf in C++ interface --- interfaces/cxx/DA.cpp | 42 +++++++++++++++++++++++++++++--- interfaces/cxx/include/dace/DA.h | 5 +++- 2 files changed, 42 insertions(+), 5 deletions(-) diff --git a/interfaces/cxx/DA.cpp b/interfaces/cxx/DA.cpp index aa6dbf7..e5b4a4b 100644 --- a/interfaces/cxx/DA.cpp +++ b/interfaces/cxx/DA.cpp @@ -317,6 +317,28 @@ DA::~DA() throw(){ /******************************************************************************** * Coefficient access routines *********************************************************************************/ +int DA::isnan() const{ +/*! Check if a DA object has any NAN coefficients. + \return True is any coefficients of the DA object are NAN. + \throw DACE::DACEException +*/ + const int temp = daceIsNan(m_index); + if(daceGetError()) DACEException(); + + return temp; +} + +int DA::isinf() const{ +/*! Check if a DA object has any INF coefficients. + \return True is any coefficients of the DA object are INF. + \throw DACE::DACEException +*/ + const int temp = daceIsInf(m_index); + if(daceGetError()) DACEException(); + + return temp; +} + double DA::cons() const{ /*! Return the constant part of a DA object. \return A double corresponding to the constant part of the DA object. @@ -1826,9 +1848,6 @@ std::istream& operator>>(std::istream &in, DA &da){ // convert string vector to DA da = DA::fromString(strs); } - - - } return in; @@ -1932,6 +1951,22 @@ void DA::memdump(){ /******************************************************************************** * DACE non-member functions *********************************************************************************/ +int isnan(const DA &da) { +/*! Check if a DA object has any NAN coefficients. + \param[in] da a given DA object. + \return True if any coefficients of the DA object are NAN. + \throw DACE::DACEException +*/ + return da.isnan();} + +int isinf(const DA &da) { +/*! Check if a DA object has any INF coefficients. + \param[in] da a given DA object. + \return True if any coefficients of the DA object are INF. + \throw DACE::DACEException +*/ + return da.isinf();} + double cons(const DA &da) { /*! Return the constant part of a DA object. \param[in] da a given DA object. @@ -1959,7 +1994,6 @@ AlgebraicVector gradient(const DA &da) { return da.gradient();} - DA divide(const DA &da, const unsigned int var, const unsigned int p){ /*! Divide by independent variable var raised to power p. 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 6bbf723..20d96c0 100644 --- a/interfaces/cxx/include/dace/DA.h +++ b/interfaces/cxx/include/dace/DA.h @@ -62,7 +62,6 @@ class DACE_API DA private: static bool initialized; //!< Indicates if DA::init() was called static std::stack TOstack; //!< Truncation order stack - // XXX: Mauro, MSVC is bitching around (also for other templated classes, e.g. std::string): Warning C4251 'DACE::DA::TOstack': class 'std::stack>>' needs to have dll - interface to be used by clients of class 'DACE::DA' DACEDA m_index; //!< Index to the DA vector public: @@ -100,6 +99,8 @@ class DACE_API DA /******************************************************************************** * Coefficient access and extraction routines *********************************************************************************/ + int isnan() const; + int isinf() const; 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 @@ -259,6 +260,8 @@ class DACE_API DA /******************************************************************************** * DACE non-member functions *********************************************************************************/ +DACE_API int isnan(const DA &da); +DACE_API int isinf(const DA &da); DACE_API double cons(const DA &da); DACE_API AlgebraicVector linear(const DA &da); DACE_API AlgebraicVector gradient(const DA &da);