diff --git a/mat.cc b/mat.cc index ace1db3..2a8c83c 100644 --- a/mat.cc +++ b/mat.cc @@ -1025,7 +1025,7 @@ std::complex tmp; #ifdef CUDALA if(location == cpu){ #endif - for(register int i=1; i (v[i][j].real()-v[j][i].real(),v[i][j].imag()+v[j][i].imag()); @@ -1094,9 +1094,9 @@ NRMat >::NRMat(const NRMat &rhs, bool imagpart): nn if(location == cpu){ #endif #ifdef MATPTR - v = new std::complex*[n]; + v = new std::complex*[nn]; v[0] = new std::complex[nn_mm]; - for(register int i=1; i)); cblas_dcopy(nn_mm, &rhs[0][0], 1, ((double *)v[0]) + (imagpart?1:0), 2); @@ -1135,9 +1135,9 @@ NRMat::NRMat(const NRMat > &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 > NRMat >::timestranspos for(i=0; i >::operator*(const NRSMat< std::complex > &r if(location == cpu){ #endif for(register int i=0; i& NRMat::swap_rows_cols(){ #ifdef CUDALA if(location == cpu){ #endif - for(register int i=0; i >& NRMat >::swap_rows_cols(){ #ifdef CUDALA if(location == cpu){ #endif - for(register int i=0;i& NRMat::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 &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 inline NRMat(const T (&a)[R][C]); @@ -110,7 +113,7 @@ public: //! explicit constructor converting vector into a NRMat object #ifdef MATPTR - explicit NRMat(const NRVec &rhs, const int n, const int m, const int offset = 0):NRMat(&rhs[0][0] + offset , n, m){ + explicit NRMat(const NRVec &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::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::is_plaindata()) memcpy(v[0], a, nn*mm*sizeof(T)); else for( int i=0; i::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::copyonwrite(v[0][i]); +#else for(int i=0; i::copyonwrite(v[i]); +#endif } if(*count > 1){ @@ -1184,7 +1191,7 @@ void NRMat::copyonwrite(bool detachonly, bool deep) { if(LA_traits::is_plaindata()) memcpy(newv[0], v[0], (size_t)mm*nn*sizeof(T)); else { - for(int i=0; i::copyonwrite(newv[i]);} + for(int i=0; i::copyonwrite(newv[0][i]);} } } v = newv; @@ -1419,7 +1426,7 @@ class NRMat_from1 : public NRMat { public: NRMat_from1(): NRMat() {}; template explicit NRMat_from1(const T (&a)[R][C]) : NRMat(a) {}; - explicit NRMat_from1(const int n): NRMat(n) {}; + explicit NRMat_from1(const int n) = delete; //single integer argument forbidden, use NRMat_from1(n,m) explicit NRMat_from1(const NRSMat_from1 &rhs) : NRMat(rhs) {}; NRMat_from1(const NRMat &rhs): NRMat(rhs) {};//!< be able to convert the parent class transparently to this NRMat_from1(const int n, const int m): NRMat(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, *) /***************************************************************************//** diff --git a/noncblas.cc b/noncblas.cc index c85fbba..e2471f7 100644 --- a/noncblas.cc +++ b/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 alphac = std::conj(*(const std::complex *)alpha); diff --git a/vecmat3.h b/vecmat3.h index d1a3333..78065a7 100644 --- a/vecmat3.h +++ b/vecmat3.h @@ -32,6 +32,7 @@ #include #ifndef NO_NUMERIC_LIMITS #include +#include #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::value,int>::type = 0> Mat3(const I x) = delete; //single integer argument forbidden (write e.g. Mat3(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));}