smat: reverted to former operator() just accessing the storage

smat.h

@@ -31,15 +31,6 @@ namespace LA {
 template<typename T> class NRPerm;
 template<typename T> class NRSMat_from1;
 
-//return type of the const element access NRSMat<T>::operator(): a reference, except for complex types,
-//where the upper triangle has to be conjugated on the fly
-template<typename T> struct NRSMat_const_element {typedef const T& type;};
-template<typename C> struct NRSMat_const_element<std::complex<C> > {typedef const std::complex<C> type;};
-
-//conjugation of the stored lower triangle element when the upper triangle element is requested
-template<typename T> inline const T& NRSMat_conj_upper(const T &x) {return x;}
-template<typename C> inline const std::complex<C> NRSMat_conj_upper(const std::complex<C> &x) {return std::conj(x);}
-
 
 /***************************************************************************//**
  * This class implements a general symmetric or hermitian matrix the elements
@@ -164,10 +155,9 @@ public:
 	inline T& operator[](const size_t ij);
 
 	//NOTE: complex NRSMat is hermitian, only the lower triangle is stored, the upper one is its complex conjugate
-	//const operator() returns conj(a(j,i)) for i<j (by value for complex types, by reference otherwise)
+	//operator() (both const and non-const) just accesses the storage: (i,j) and (j,i) refer to the same stored element,
-	//PITFALL: the non-const operator() cannot conjugate (it would need a smart pointer for the l-value), for i<j it returns
+	//NO conjugation is performed for i<j; the hermitian upper element a(i,j), i<j, is conj(a(j,i))
-	//a reference to the stored element a(j,i) WITHOUT conjugation - this also applies to reading from a non-const complex matrix
+	inline const T& operator()(const int i, const int j) const;
-	inline typename NRSMat_const_element<T>::type operator()(const int i, const int j) const;
 	inline T& operator()(const int i, const int j);
 
 	inline int nrows() const;
@@ -784,9 +774,8 @@ return (i>j)? (i-2)*(i-1)/2+j-1 : (j-2)*(j-1)/2+i-1;
 
 /***************************************************************************//**
  * determine matrix element of this symmetric matrix of general type <code>T</code>
- * PITFALL for complex (hermitian) matrices: for \f$i<j\f$ the reference to the stored element \f$a_{ji}\f$
+ * for complex (hermitian) matrices: for \f$i<j\f$ the reference to the stored element \f$a_{ji}\f$
- * is returned WITHOUT complex conjugation (an l-value cannot be conjugated); this also
+ * is returned WITHOUT complex conjugation, operator() just accesses the storage
- * applies when only reading from a non-const matrix, use the const version or \f$i\geq{}j\f$
  * @param[in] i row index running from 0
  * @param[in] j column index running from 0
  * @return reference to the corresponding matrix element
@@ -810,15 +799,16 @@ inline T & NRSMat<T>::operator()(const int i, const int j) {
 }
 
 /***************************************************************************//**
- * determine matrix element of this symmetric or hermitian matrix of general type <code>T</code>
+ * determine matrix element of this symmetric matrix of general type <code>T</code>
+ * for complex (hermitian) matrices: for \f$i<j\f$ the stored element \f$a_{ji}\f$ is returned
+ * WITHOUT complex conjugation, operator() just accesses the storage
  * @param[in] i row index running from 0
  * @param[in] j column index running from 0
- * @return constant reference to the corresponding matrix element; for complex types the value,
+ * @return constant reference to the corresponding matrix element
- * complex conjugated for \f$i<j\f$ (hermitian matrix)
  * @see count, SMat_index, NRSMat<T>::operator[]
  ******************************************************************************/
 template <typename T>
-inline typename NRSMat_const_element<T>::type NRSMat<T>::operator()(const int i, const int j) const {
+inline const T & NRSMat<T>::operator()(const int i, const int j) const {
 #ifdef DEBUG
 	if(i<0||i>=nn||j<0||j>=nn) std::cout<<"INDEX PROBLEM "<<0<<" "<<i<<" "<<nn-1<<" "<<0<<" "<<j<<" "<<nn-1<<std::endl;
         if(i<0) laerror("T & NRSMat<T>::operator()(const int, const int) first index out of range - low");
@@ -830,7 +820,7 @@ inline typename NRSMat_const_element<T>::type NRSMat<T>::operator()(const int i,
 #endif
 	NOT_GPU(*this);
 
-	return (i>=j) ? v[SMat_index_igej(i,j)] : NRSMat_conj_upper(v[SMat_index_igej(j,i)]);
+	return v[SMat_index(i,j)];
 }
 
 /***************************************************************************//**
@@ -1357,7 +1347,7 @@ public:
 	explicit NRSMat_from1(const NRMat<T> &rhs): NRSMat<T>(rhs) {};
 	explicit NRSMat_from1(const NRVec<T> &rhs, const int n): NRSMat<T>(rhs,n) {};
 
-	inline typename NRSMat_const_element<T>::type operator() (const int i, const int j) const { //for complex conjugated for i<j, see NRSMat<T>::operator()
+	inline const T& operator() (const int i, const int j) const { //storage access, no conjugation for complex i<j, see NRSMat<T>::operator()
 	#ifdef DEBUG
 		if(i<=0||i>NRSMat<T>::nn||j<=0||j>NRSMat<T>::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<<i<<" "<<NRSMat<T>::nn<<" "<<1<<" "<<j<<" "<<NRSMat<T>::nn<<std::endl;
 		if(i<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) first indexout of range - low");
@@ -1365,10 +1355,10 @@ public:
 		if(j<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) second index out of range - low");
 		if(j>NRSMat<T>::nn) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) second index out of range - high");
 	#endif
-		return (i>=j) ? NRSMat<T>::v[SMat_index_1igej(i,j)] : NRSMat_conj_upper(NRSMat<T>::v[SMat_index_1igej(j,i)]);
+		return NRSMat<T>::v[SMat_index_1(i,j)];
 	}
 	
-	inline  T& operator() (const int i, const int j){ //PITFALL: no conjugation for complex i<j, see NRSMat<T>::operator()
+	inline  T& operator() (const int i, const int j){ //storage access, no conjugation for complex i<j, see NRSMat<T>::operator()
 	#ifdef DEBUG
 		if(i<=0||i>NRSMat<T>::nn||j<=0||j>NRSMat<T>::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<<i<<" "<<NRSMat<T>::nn<<" "<<1<<" "<<j<<" "<<NRSMat<T>::nn<<std::endl;
                 if(i<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) first indexout of range - low");