implemented diffabs() useful for checks of results up to a sign

mat.cc

@@ -1429,29 +1429,6 @@ NRMat<std::complex<double> >::operator+=(const NRMat< std::complex<double> >  &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 <typename T>
-NRMat<T> & NRMat<T>::operator+=(const NRMat<T>  &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<double> >::operator-=(const NRMat< std::complex<double> >  &
 }
 
 
-/***************************************************************************//**
- * 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 <typename T>
-NRMat<T> & NRMat<T>::operator-=(const NRMat<T>  &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

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<T>::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 <typename T>
+NRMat<T> & NRMat<T>::operator+=(const NRMat<T>  &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 <typename T>
+NRMat<T> & NRMat<T>::operator-=(const NRMat<T>  &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 <typename T>
+NRMat<T>  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
+        return r;
+}
+
+
+
+
 
  
 

polynomial.h

@@ -26,12 +26,6 @@
 
 namespace LA {
 
-template <typename T>
-inline  typename LA_traits<T>::normtype  MYABS(const T &x) {return abs(x);}
-
-template <>
-inline unsigned int MYABS(const  unsigned int &x) {return x;}
-
 
 template <typename T>
 class Polynomial : public NRVec<T> {

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<T> & NRSMat<T>::operator-=(const T &a) {
 	return *this;
 }
 
+/*difference of absolute values*/
+template <typename T>
+NRSMat<T>  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]);
+        return r;
+}
+
+
 /***************************************************************************//**
  * add up this real symmetric matrix with given symmetric matrix
  * @param[in] rhs real symmetric matrix to be added

vec.h

@@ -33,6 +33,20 @@ template <typename T> void lawritemat(FILE *file, const T *a, int r, int c,
 template <typename T> class NRPerm;
 template <typename T> class CyclePerm;
 
+/***************************************************************************//**
+ * auxiliary macro to avoid compilation errors for some types
+ ******************************************************************************/
+template <typename T>
+inline  typename LA_traits<T>::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<T> & NRVec<T>::operator-=(const NRVec<T> &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 <typename T>
+NRVec<T>  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);
+	for(int i=0; i<nn; ++i) r[i] = MYABS(v[i]) - MYABS(rhs.v[i]);
+        return r;
+}
+
 
 /***************************************************************************//**
  * multiply this general vector \f$\vec{x}\f$ by scalar value \f$\lambda\f$