From b1f415fe047e72e122f8c2b7a2d2d3cf997f94e4 Mon Sep 17 00:00:00 2001
From: jiri <jiri>
Date: Wed, 16 Apr 2008 12:56:02 +0000
Subject: [PATCH] *** empty log message ***

---
 mat.cc | 45 +++++++--------------------------------------
 vec.h  | 14 ++++----------
 2 files changed, 11 insertions(+), 48 deletions(-)

diff --git a/mat.cc b/mat.cc
index e467085..f393596 100644
--- a/mat.cc
+++ b/mat.cc
@@ -913,24 +913,6 @@ NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
 
 
 
-// sum of rows
-template<>
-const NRVec<double> NRMat<double>::rsum() const
-{
-	NRVec<double> result(mm);
-	for (int i=0; i<mm; i++) result[i] = cblas_dasum(nn,(*this)[0]+i,mm);
-	return result;
-}
-
-// sum of columns
-template<>
-const NRVec<double> NRMat<double>::csum() const
-{
-	NRVec<double> result(nn);
-	for (int i=0; i<nn; i++) result[i] = cblas_dasum(mm, (*this)[i], 1);
-	return result;
-}
-
 // complex conjugate of Mat
 template<>
 NRMat<double> &NRMat<double>::conjugateme() {return *this;}
@@ -1082,37 +1064,24 @@ void NRMat< complex<double> >::axpy(const complex<double> alpha,
 
 
 // trace of Mat
-template<>
-const double NRMat<double>::trace() const
+template <typename T>
+const T NRMat<T>::trace() const
 {
 #ifdef DEBUG
 	if (nn != mm) laerror("no-square matrix in Mat::trace()");
 #endif
-	return cblas_dasum(nn, (*this)[0], nn+1);
-
-
-}
-
-
-
-template<>
-const complex<double> NRMat< complex<double> >::trace() const
-{
-#ifdef DEBUG
-	if (nn != mm) laerror("no-square matrix in Mat::trace()");
-#endif
-	register complex<double> sum = CZERO;
-	for (int i=0; i<nn*nn; i+=(nn+1))
+T sum=0;
 #ifdef MATPTR
-		sum += v[0][i];
+for (int i=0; i<nn; ++i) sum += v[i][i];
 #else
-		sum += v[i];
+for (int i=0; i<nn*nn; i+=(nn+1)) sum += v[i];
 #endif
-	return sum;
+return sum;
 }
 
 
 
+
 //get diagonal; for compatibility with large matrices do not return newly created object
 //for non-square get diagonal of A^TA, will be used as preconditioner 
 template<>
diff --git a/vec.h b/vec.h
index 33c4d86..4a61365 100644
--- a/vec.h
+++ b/vec.h
@@ -103,7 +103,8 @@ public:
 	const NRVec operator*(const NRSMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
 	const NRVec operator*(const SparseMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
 	const NRMat<T> operator|(const NRVec<T> &rhs) const;
-	inline const T sum() const; //sum of its elements
+	inline const T sum() const {T sum=0; for(int i=0; i<nn; i++) sum += v[i]; return sum;}; //sum of its elements
+	inline const T asum() const; //sum of its elements absolute values
 	inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector
 	inline T & operator[](const int i);
 	inline const T & operator[](const int i) const;
@@ -447,18 +448,11 @@ inline const T NRVec<T>::operator*(const NRVec<T> &rhs) const
 
 // Sum of elements
 template<>
-inline const double NRVec<double>::sum() const
+inline const double NRVec<double>::asum() const
 {
 	return cblas_dasum(nn, v, 1);
 }
-template<>
-inline const complex<double>
-NRVec< complex<double> >::sum() const
-{
-	complex<double> sum = CZERO;
-	for (int i=0; i<nn; i++) sum += v[i];
-	return sum;
-}
+
 
 // Dot product: x * y
 template<>