From cf86493a6ff91b8ffc7f9de7c17ea26c6464b7f1 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Wed, 30 Jun 2021 14:54:35 +0200 Subject: [PATCH] implemented diffabs() useful for checks of results up to a sign --- mat.cc | 46 ---------------------------------- mat.h | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++++ polynomial.h | 6 ----- smat.h | 16 ++++++++++++ vec.h | 35 ++++++++++++++++++++++++++ 5 files changed, 121 insertions(+), 52 deletions(-) diff --git a/mat.cc b/mat.cc index 0bbec47..a5f315c 100644 --- a/mat.cc +++ b/mat.cc @@ -1429,29 +1429,6 @@ NRMat >::operator+=(const NRMat< std::complex > &r return *this; } -/***************************************************************************//** - * add a given general matrix (type T) \f$A\f$ to the current complex matrix - * @param[in] rhs matrix \f$A\f$ of type T - * @return reference to the modified matrix - ******************************************************************************/ -template -NRMat & NRMat::operator+=(const NRMat &rhs) { -#ifdef DEBUG - if (nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices"); -#endif - SAME_LOC(*this, rhs); - NOT_GPU(*this); - - copyonwrite(); - - #ifdef MATPTR - for(size_t i=0; i< (size_t)nn*mm; i++) v[0][i] += rhs.v[0][i]; - #else - for(size_t i=0; i< (size_t)nn*mm; i++) v[i] += rhs.v[i]; - #endif - return *this; -} - /***************************************************************************//** * subtract a given real matrix \f$A\f$ from the current real matrix @@ -1505,29 +1482,6 @@ NRMat< std::complex >::operator-=(const NRMat< std::complex > & } -/***************************************************************************//** - * subtract a given general matrix (type T) \f$A\f$ from the current matrix - * @param[in] rhs matrix \f$A\f$ of type T - * @return reference to the modified matrix - ******************************************************************************/ -template -NRMat & NRMat::operator-=(const NRMat &rhs) { -#ifdef DEBUG - if (nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices"); -#endif - SAME_LOC(*this, rhs); - NOT_GPU(*this); - - copyonwrite(); - - #ifdef MATPTR - for(size_t i=0; i< (size_t)nn*mm; i++) v[0][i] += rhs.v[0][i]; - #else - for(size_t i=0; i<(size_t) nn*mm; i++) v[i] += rhs.v[i]; - #endif - return *this; -} - /***************************************************************************//** * add a given sparse real matrix \f$A\f$ stored in packed form to the current diff --git a/mat.h b/mat.h index 936e89c..1b96dd3 100644 --- a/mat.h +++ b/mat.h @@ -406,6 +406,7 @@ public: #endif } + NRMat diffabs(const NRMat &rhs) const; //difference of absolute values }; }//namespace @@ -1402,6 +1403,75 @@ void NRMat::moveto(const GPUID dest) { #endif +/***************************************************************************//** + * add a given general matrix (type T) \f$A\f$ to the current complex matrix + * @param[in] rhs matrix \f$A\f$ of type T + * @return reference to the modified matrix + ******************************************************************************/ +template +NRMat & NRMat::operator+=(const NRMat &rhs) { +#ifdef DEBUG + if (nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices"); +#endif + SAME_LOC(*this, rhs); + NOT_GPU(*this); + + copyonwrite(); + + #ifdef MATPTR + for(size_t i=0; i< (size_t)nn*mm; i++) v[0][i] += rhs.v[0][i]; + #else + for(size_t i=0; i< (size_t)nn*mm; i++) v[i] += rhs.v[i]; + #endif + return *this; +} + + +/***************************************************************************//** + * subtract a given general matrix (type T) \f$A\f$ from the current matrix + * @param[in] rhs matrix \f$A\f$ of type T + * @return reference to the modified matrix + ******************************************************************************/ +template +NRMat & NRMat::operator-=(const NRMat &rhs) { +#ifdef DEBUG + if (nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices"); +#endif + SAME_LOC(*this, rhs); + NOT_GPU(*this); + + copyonwrite(); + + #ifdef MATPTR + for(size_t i=0; i< (size_t)nn*mm; i++) v[0][i] -= rhs.v[0][i]; + #else + for(size_t i=0; i<(size_t) nn*mm; i++) v[i] -= rhs.v[i]; + #endif + return *this; +} + + +/*difference of absolute values*/ +template +NRMat NRMat::diffabs(const NRMat &rhs) const { +#ifdef DEBUG + if (nn != rhs.nn ||mm!=rhs.mm) laerror("incompatible dimensions"); +#endif + NOT_GPU(*this); + NOT_GPU(rhs); + +NRMat r(nn,mm); + #ifdef MATPTR + for(size_t i=0; i< (size_t)nn*mm; i++) r.v[0][i] = MYABS(v[0][i]) - MYABS(rhs.v[0][i]); + #else + for(size_t i=0; i<(size_t) nn*mm; i++) r.v[i] = MYABS(v[i]) - MYABS(rhs.v[i]); + #endif + return r; +} + + + + diff --git a/polynomial.h b/polynomial.h index 73e9678..a41598b 100644 --- a/polynomial.h +++ b/polynomial.h @@ -26,12 +26,6 @@ namespace LA { -template -inline typename LA_traits::normtype MYABS(const T &x) {return abs(x);} - -template <> -inline unsigned int MYABS(const unsigned int &x) {return x;} - template class Polynomial : public NRVec { diff --git a/smat.h b/smat.h index e2d34e0..052c7a1 100644 --- a/smat.h +++ b/smat.h @@ -186,6 +186,7 @@ public: #endif } + NRSMat diffabs(const NRSMat &rhs) const; //difference of absolute values }; }//namespace @@ -389,6 +390,21 @@ inline NRSMat & NRSMat::operator-=(const T &a) { return *this; } +/*difference of absolute values*/ +template +NRSMat NRSMat::diffabs(const NRSMat &rhs) const { +#ifdef DEBUG + if (nn != rhs.nn) laerror("incompatible dimensions"); +#endif + NOT_GPU(*this); + NOT_GPU(rhs); + +NRSMat r(nn); + for(int i=0; i void lawritemat(FILE *file, const T *a, int r, int c, template class NRPerm; template class CyclePerm; +/***************************************************************************//** + * auxiliary macro to avoid compilation errors for some types + ******************************************************************************/ +template +inline typename LA_traits::normtype MYABS(const T &x) {return abs(x);} + +template <> inline unsigned char MYABS(const unsigned char &x) {return x;} +template <> inline unsigned short MYABS(const unsigned short &x) {return x;} +template <> inline unsigned int MYABS(const unsigned int &x) {return x;} +template <> inline unsigned long MYABS(const unsigned long &x) {return x;} +template <> inline unsigned long long MYABS(const unsigned long long &x) {return x;} + + + /***************************************************************************//** * static constants used in several cblas-routines ******************************************************************************/ @@ -394,6 +408,8 @@ public: #endif } + + NRVec diffabs(const NRVec &rhs) const; //difference of absolute values }; @@ -759,6 +775,25 @@ inline NRVec & NRVec::operator-=(const NRVec &rhs) { return *this; } +/***************************************************************************//** + * difference of elements of two vectors in absolute values + * \f[\vec{z}_i = \vec{x}_i-\vec{y}_i\f] + * @param[in] rhs vector \f$\vec{y}\f$ + * @return reference to the modified vector + ******************************************************************************/ +template +NRVec NRVec::diffabs(const NRVec &rhs) const { +#ifdef DEBUG + if (nn != rhs.nn) laerror("incompatible dimensions"); +#endif + NOT_GPU(*this); + NOT_GPU(rhs); + +NRVec r(nn); + for(int i=0; i