bugfixes in mat.h mat.cc and related

mat.cc

@@ -1025,7 +1025,7 @@ std::complex<double> tmp;
 #ifdef CUDALA
 	if(location == cpu){
 #endif
-		for(register int i=1; i<nn; i++){
+		for(register int i=0; i<nn; i++){
 			for(register int j=0; j<=i; j++){
 			#ifdef MATPTR
 				tmp = std::complex<double> (v[i][j].real()-v[j][i].real(),v[i][j].imag()+v[j][i].imag()); 
@@ -1094,9 +1094,9 @@ NRMat<std::complex<double> >::NRMat(const NRMat<double> &rhs, bool imagpart): nn
 	if(location == cpu){
 #endif
 	#ifdef MATPTR
-	        v = new std::complex<double>*[n];
+	        v = new std::complex<double>*[nn];
 	        v[0] = new std::complex<double>[nn_mm];
-	        for(register int i=1; i<n; i++) v[i] = v[i-1] + m;
+	        for(register int i=1; i<nn; i++) v[i] = v[i-1] + mm;
 	
 		memset(v[0], 0, nn_mm*sizeof(std::complex<double>));
 	        cblas_dcopy(nn_mm, &rhs[0][0], 1, ((double *)v[0]) + (imagpart?1:0), 2);
@@ -1135,9 +1135,9 @@ NRMat<double>::NRMat(const NRMat<std::complex<double> > &rhs, bool imagpart): nn
 	if(location == cpu){
 #endif
 	#ifdef MATPTR
-	        v = new double*[n];
+	        v = new double*[nn];
 	        v[0] = new double[nn_mm];
-	        for(register int i=1; i<n; i++) v[i] = v[i-1] + m;
+	        for(register int i=1; i<nn; i++) v[i] = v[i-1] + mm;
 	
 	        cblas_dcopy(nn_mm, ((double *)&rhs[0][0]) + (imagpart?1:0), 2, v[0], 1);
 	#else
@@ -1357,9 +1357,9 @@ const NRSMat<std::complex<double> >  NRMat<std::complex<double> >::timestranspos
 		for(i=0; i<nn; ++i){
 			for(j=0; j<=i; ++j){
 			#ifdef MATPTR
-				cblas_zdotc_sub(nn, v[i], 1, v[j], 1, &r(i,j));
+				cblas_zdotc_sub(mm, v[j], 1, v[i], 1, &r(i,j)); //(A A^H)_ij = sum_k A_ik conj(A_jk), lower triangle of the hermitian result
 			#else
-				cblas_zdotc_sub(nn, v + i*(size_t)mm, 1, v + j*(size_t)mm, 1, &r(i,j));
+				cblas_zdotc_sub(mm, v + j*(size_t)mm, 1, v + i*(size_t)mm, 1, &r(i,j)); //(A A^H)_ij = sum_k A_ik conj(A_jk), lower triangle of the hermitian result
 			#endif
 			}
 		}
@@ -2162,7 +2162,8 @@ NRMat< std::complex<double> >::operator*(const NRSMat< std::complex<double> > &r
 	if(location == cpu){
 #endif
 		for(register int i=0; i<nn; i++){
-			cblas_zhpmv(CblasRowMajor, CblasLower, mm, &CONE, &rhs[0],
+			//row i of A H is H^T a_i; the C-order lower packed hermitian H read as fortran-order upper packed is H^T
+			cblas_zhpmv(CblasColMajor, CblasUpper, mm, &CONE, &rhs[0],
 					(*this)[i], 1, &CZERO, result[i], 1);
 		}
 #ifdef CUDALA
@@ -3143,24 +3144,9 @@ NRMat<double>& NRMat<double>::swap_rows_cols(){
 #ifdef CUDALA
 	if(location == cpu){
 #endif	
-		for(register int i=0; i<n_pul; i++){
+		//reversal of the whole C-order array; for a negative increment BLAS expects the pointer to the first element in memory
-			std::cout << "swapping row " << i << " and " << nn-i-1 << std::endl;
+		const size_t dim = (size_t)nn*mm;
-
+		cblas_dswap(dim/2, (*this)[0], 1, (*this)[0] + dim - dim/2, -1);
-			std::cout << "elements: " << *((*this)[i]) << std::endl;
-			std::cout << "elements: " << *((*this)[nn - i - 1] + mm - 1) << std::endl;
-
-			cblas_dswap(mm, (*this)[i], 1, (*this)[nn - i - 1] + mm - 1, -1);
-		}
-		return *this;
-
-		if(nn & 1){ // odd number of rows
-			for(register int i=0; i<=m_pul; i++){
-
-				tmp = (*this)(n_pul, i);
-				(*this)(n_pul, i) = (*this)(n_pul, mm-i-1);
-				(*this)(n_pul, mm-i-1) = tmp;
-			}
-		}
 #ifdef CUDALA
 	}else{
 		for(register int i=0; i<n_pul; i++){
@@ -3197,17 +3183,9 @@ NRMat<std::complex<double> >& NRMat<std::complex<double> >::swap_rows_cols(){
 #ifdef CUDALA
 	if(location == cpu){
 #endif	
-		for(register int i=0;i<n_pul;i++){
+		//reversal of the whole C-order array; for a negative increment BLAS expects the pointer to the first element in memory
-			cblas_zswap(mm, (*this)[i], 1, (*this)[nn - i - 1] + mm - 1, -1);
+		const size_t dim = (size_t)nn*mm;
-		}
+		cblas_zswap(dim/2, (*this)[0], 1, (*this)[0] + dim - dim/2, -1);
-
-		if(nn & 1){
-			for(register int i=0; i<=m_pul; i++){ // odd number of rows
-				tmp = (*this)(n_pul, i);
-				(*this)(n_pul, i) = (*this)(n_pul, mm-i-1);
-				(*this)(n_pul, mm-i-1) = tmp;
-			}
-		}
 #ifdef CUDALA
 	}else{
 		for(register int i=0;i<n_pul;i++){
@@ -3247,7 +3225,7 @@ NRMat<T>& NRMat<T>::swap_rows_cols(){
 #endif
 		data_ptr = (*this)[0];	
 		const int dim_pul = dim >> 1;
-		for(register int i=0; i<=dim_pul; i++){
+		for(register int i=0; i<dim_pul; i++){
 			tmp = data_ptr[i];
 			data_ptr[i] = data_ptr[dim - i - 1];
 			data_ptr[dim - i - 1] = tmp;

mat.h

@@ -93,6 +93,9 @@ public:
 	//! inlined constructor creating matrix of given size filled with data located at given memory location 
 	inline NRMat(const T *a, const int n, const int m);
 	inline NRMat(const Mat3<T> &rhs) : NRMat(&rhs(0,0),3,3) {};
+
+	//! a single integer argument is forbidden (it would silently create a 3x3 matrix via Mat3(const T)); use NRMat(n,m)
+	explicit NRMat(const int n) = delete;
 	
 	//! inlined constructor creating matrix of given size from an array
 	template<int R, int C> inline NRMat(const T (&a)[R][C]);
@@ -110,7 +113,7 @@ public:
 
 	//! explicit constructor converting vector into a <code>NRMat<T></code> object
 #ifdef MATPTR
-	explicit NRMat(const NRVec<T> &rhs, const int n, const int m, const int offset = 0):NRMat(&rhs[0][0] + offset , n, m){
+	explicit NRMat(const NRVec<T> &rhs, const int n, const int m, const int offset = 0):NRMat(&rhs[0] + offset , n, m){
 		if (offset < 0 || (size_t)n*m + offset > rhs.nn) laerror("matrix dimensions and offset incompatible with vector length");
 	};
 #else
@@ -561,7 +564,7 @@ inline NRMat<T>::NRMat(const T (&a)[R][C]) : count(new int) {
 	#ifdef MATPTR
 		v = new T*[nn];
 		v[0] = new T[nn*mm];
-		for (register int i=1; i<n; i++) v[i] = v[i-1] + m;
+		for (register int i=1; i<nn; i++) v[i] = v[i-1] + mm;
 		if(LA_traits<T>::is_plaindata()) memcpy(v[0], a, nn*mm*sizeof(T));
 		else for( int i=0; i<nn; i++) for(int j=0; j<mm; ++j) v[i][j]=a[i][j];
 	#else
@@ -1166,7 +1169,11 @@ void NRMat<T>::copyonwrite(bool detachonly, bool deep) {
 #ifdef CUDALA
                 if(location != cpu) laerror("nested types not supported on gpu memory");
 #endif
+#ifdef MATPTR
+		for(int i=0; i<nn*mm; ++i) LA_traits<T>::copyonwrite(v[0][i]);
+#else
 		for(int i=0; i<nn*mm; ++i) LA_traits<T>::copyonwrite(v[i]);
+#endif
                 }
 
 	if(*count > 1){
@@ -1184,7 +1191,7 @@ void NRMat<T>::copyonwrite(bool detachonly, bool deep) {
 				if(LA_traits<T>::is_plaindata()) memcpy(newv[0], v[0], (size_t)mm*nn*sizeof(T));
                                 else
                                         {
-                                        for(int i=0; i<nn*mm; ++i) {newv[i]=v[i]; LA_traits<T>::copyonwrite(newv[i]);}
+                                        for(int i=0; i<nn*mm; ++i) {newv[0][i]=v[0][i]; LA_traits<T>::copyonwrite(newv[0][i]);}
                                         }
 				}
 			v = newv;
@@ -1419,7 +1426,7 @@ class NRMat_from1 : public NRMat<T> {
 public:
 	NRMat_from1(): NRMat<T>() {};
 	template<int R, int C> explicit NRMat_from1(const T (&a)[R][C]) : NRMat<T>(a) {};
-	explicit NRMat_from1(const int n): NRMat<T>(n) {};
+	explicit NRMat_from1(const int n) = delete; //single integer argument forbidden, use NRMat_from1(n,m)
 	explicit NRMat_from1(const NRSMat_from1<T> &rhs) : NRMat<T>(rhs) {};
 	NRMat_from1(const NRMat<T> &rhs): NRMat<T>(rhs) {};//!< be able to convert the parent class transparently to this
 	NRMat_from1(const int n, const int m): NRMat<T>(n, m) {};
@@ -1658,7 +1665,7 @@ tmp.printsorted(s,direction,absvalue,stable,thr,topmax,&indexmap);
  * corresponding macro is defined in vec.h
  ******************************************************************************/
 NRVECMAT_OPER(Mat, +)
-NRVECMAT_OPER(Mat, -)
+NRVECMAT_OPER_NONCOMMUTATIVE(Mat, -)
 NRVECMAT_OPER(Mat, *)
 
 /***************************************************************************//**

noncblas.cc

@@ -347,9 +347,21 @@ void cblas_zhpmv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
 		const void *X, const int incX,
 		const void *beta, void *Y, const int incY)
 {
-	if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted");
-	if(Uplo!=CblasLower) laerror("CblasLower uplo asserted");
 	char U = BLAS_FORTRANCASE('u');
+	if(Order==CblasColMajor && Uplo==CblasUpper) //native fortran convention, no conjugation needed
+		{
+#ifdef FORINT
+		const FINT ntmp=N;
+		const FINT incxtmp=incX;
+		const FINT incytmp=incY;
+		FORNAME(zhpmv) (&U,&ntmp, alpha, Ap, X, &incxtmp, beta, Y, &incytmp);
+#else
+		FORNAME(zhpmv) (&U,&N, alpha, Ap, X, &incX, beta, Y, &incY);
+#endif
+		return;
+		}
+	if(Order!=CblasRowMajor) laerror("CblasRowMajor or CblasColMajor+CblasUpper order asserted");
+	if(Uplo!=CblasLower) laerror("CblasLower uplo asserted");
 //row-major lower packed hermitian H is column-major upper packed conj(H), therefore
 //conj(y) = conj(beta) conj(y) + conj(alpha) conj(H) conj(x) is computed by zhpmv on conjugated vectors and scalars
 std::complex<double> alphac = std::conj(*(const std::complex<double> *)alpha);

vecmat3.h

@@ -32,6 +32,7 @@
 #include <stdio.h>
 #ifndef NO_NUMERIC_LIMITS
 #include <limits>
+#include <type_traits>
 #endif
 
 namespace LA_Vecmat3 {
@@ -136,6 +137,7 @@ public:
 	Mat3(void) {};
 	Mat3(const T (&a)[3][3]) {memcpy(q,a,3*3*sizeof(T));}
 	Mat3(const T x) {memset(q,0,9*sizeof(T)); q[0][0]=q[1][1]=q[2][2]=x;}; //scalar matrix
+	template<typename I, typename std::enable_if<std::is_integral<I>::value,int>::type = 0> Mat3(const I x) = delete; //single integer argument forbidden (write e.g. Mat3<double>(1.))
 	Mat3& operator=(const T &x) {memset(q,0,9*sizeof(T)); q[0][0]=q[1][1]=q[2][2]=x; return *this;}; //scalar matrix
 	void indentity() {*this = (T)1;};
 	Mat3(const T* x) {memcpy(q,x,9*sizeof(T));}