From b1fea2f1c2f0baccf8952cedf4698b4b66a61e2f Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Sat, 13 Nov 2021 19:00:46 +0100
Subject: [PATCH] diffabs return normtype entity, elementwise abs() implemented

---
 mat.h  | 23 +++++++++++++++--------
 smat.h | 20 ++++++++++++++++----
 vec.h  | 16 +++++++++++++---
 3 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/mat.h b/mat.h
index 512843f..1e0ff9b 100644
--- a/mat.h
+++ b/mat.h
@@ -409,7 +409,8 @@ public:
 #endif
 
 		}
-	NRMat diffabs(const NRMat &rhs) const; //difference of absolute values
+	NRMat<typename LA_traits<T>::normtype> abs() const; //elementwise absolute values
+	NRMat<typename LA_traits<T>::normtype> diffabs(const NRMat &rhs) const; //difference of absolute values
 };
 
 }//namespace
@@ -1505,19 +1506,25 @@ NRMat<T> & NRMat<T>::operator-=(const NRMat<T>  &rhs) {
 
 /*difference of absolute values*/
 template <typename T>
-NRMat<T>  NRMat<T>::diffabs(const NRMat<T> &rhs) const {
+NRMat<typename LA_traits<T>::normtype>  NRMat<T>::diffabs(const NRMat<T> &rhs) const {
 #ifdef DEBUG
         if (nn != rhs.nn ||mm!=rhs.mm) laerror("incompatible dimensions");
 #endif
         NOT_GPU(*this);
         NOT_GPU(rhs);
 
-NRMat<T> 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
+	NRMat<typename LA_traits<T>::normtype> r(nn,mm);
+        for(size_t i=0; i< nn; i++)  for(size_t j=0; j< mm; j++) r(i,j) = MYABS((*this)(i,j)) - MYABS(rhs(i,j));
+        return r;
+}
+
+/*elementwise  absolute values*/
+template <typename T>
+NRMat<typename LA_traits<T>::normtype>  NRMat<T>::abs() const {
+        NOT_GPU(*this);
+
+        NRMat<typename LA_traits<T>::normtype> r(nn,mm);
+        for(size_t i=0; i< nn; i++)  for(size_t j=0; j< mm; j++) r(i,j) = MYABS((*this)(i,j));
         return r;
 }
 
diff --git a/smat.h b/smat.h
index 975e5b4..621e33a 100644
--- a/smat.h
+++ b/smat.h
@@ -186,7 +186,8 @@ public:
 #endif
 
 		}
-	NRSMat diffabs(const NRSMat &rhs) const; //difference of absolute values
+	NRSMat<typename LA_traits<T>::normtype> diffabs(const NRSMat &rhs) const; //difference of absolute values
+	NRSMat<typename LA_traits<T>::normtype> abs() const; //elementwise absolute values
 };
 
 }//namespace
@@ -395,18 +396,29 @@ inline NRSMat<T> & NRSMat<T>::operator-=(const T &a) {
 
 /*difference of absolute values*/
 template <typename T>
-NRSMat<T>  NRSMat<T>::diffabs(const NRSMat<T> &rhs) const {
+NRSMat<typename LA_traits<T>::normtype>  NRSMat<T>::diffabs(const NRSMat<T> &rhs) const {
 #ifdef DEBUG
         if (nn != rhs.nn) laerror("incompatible dimensions");
 #endif
         NOT_GPU(*this);
         NOT_GPU(rhs);
 
-NRSMat<T> r(nn);
-        for(int i=0; i<NN2; ++i) r.v[i] = MYABS(v[i]) - MYABS(rhs.v[i]);
+NRSMat<typename LA_traits<T>::normtype> r(nn);
+        for(int i=0; i<nn; ++i) for(int j=0; j<=i; ++j) r(i,j) = MYABS((*this)(i,j)) - MYABS(rhs(i,j));
         return r;
 }
 
+/*elementwise absolute values*/
+template <typename T>
+NRSMat<typename LA_traits<T>::normtype>  NRSMat<T>::abs() const {
+        NOT_GPU(*this);
+
+NRSMat<typename LA_traits<T>::normtype> r(nn);
+        for(int i=0; i<nn; ++i) for(int j=0; j<=i; ++j) r(i,j) = MYABS((*this)(i,j));
+        return r;
+}
+
+
 
 /***************************************************************************//**
  * add up this real symmetric matrix with given symmetric matrix
diff --git a/vec.h b/vec.h
index 84ca3cc..4debbab 100644
--- a/vec.h
+++ b/vec.h
@@ -426,7 +426,8 @@ public:
 
 		}
 
-	NRVec diffabs(const NRVec &rhs) const; //difference of absolute values
+	NRVec<typename LA_traits<T>::normtype> diffabs(const NRVec &rhs) const; //difference of absolute values
+	NRVec<typename LA_traits<T>::normtype> abs() const; //element-wise absolute values
 };
 
 
@@ -834,18 +835,27 @@ inline NRVec<T> & NRVec<T>::operator-=(const NRVec<T> &rhs) {
  * @return reference to the modified vector
  ******************************************************************************/
 template <typename T>
-NRVec<T>  NRVec<T>::diffabs(const NRVec<T> &rhs) const {
+NRVec<typename LA_traits<T>::normtype>  NRVec<T>::diffabs(const NRVec<T> &rhs) const {
 #ifdef DEBUG
 	if (nn != rhs.nn) laerror("incompatible dimensions");
 #endif
 	NOT_GPU(*this);
 	NOT_GPU(rhs);
 
-NRVec<T> r(nn);
+NRVec<typename LA_traits<T>::normtype> r(nn);
 	for(int i=0; i<nn; ++i) r[i] = MYABS(v[i]) - MYABS(rhs.v[i]);
         return r;
 }
 
+template <typename T>
+NRVec<typename LA_traits<T>::normtype>  NRVec<T>::abs() const {
+        NOT_GPU(*this);
+
+NRVec<typename LA_traits<T>::normtype> r(nn);
+        for(int i=0; i<nn; ++i) r[i] = MYABS(v[i]);
+        return r;
+}
+
 
 /***************************************************************************//**
  * multiply this general vector \f$\vec{x}\f$ by scalar value \f$\lambda\f$