diff --git a/mat.cc b/mat.cc index 2a8c83c..84d608f 100644 --- a/mat.cc +++ b/mat.cc @@ -1705,8 +1705,8 @@ NRMat< std::complex >::operator+=(const NRSMat< std::complex > & p += i + 1; } p = rhs; p++; - for(register int i=1; i & NRMat::operator+=(const NRSMat &rhs) { } p = rhs; p++; for(register int i=1; i::conjugate(p[j]); //hermitian for complex p += i+1; } return *this; @@ -1829,8 +1829,8 @@ NRMat >::operator-=(const NRSMat > &rh p += i + 1; } p = rhs; p++; - for(register int i=1; i & NRMat::operator-=(const NRSMat &rhs) { } p = rhs; p++; for(register int i=1; i::conjugate(p[j]); //hermitian for complex p += i+1; } return *this; diff --git a/mat.h b/mat.h index 1bcb7be..95fbcff 100644 --- a/mat.h +++ b/mat.h @@ -675,7 +675,7 @@ NRMat::NRMat(const NRMat &rhs) { } /***************************************************************************//** - * create matrix from a \c NRSMat object + * create matrix from a \c NRSMat object (for complex types hermitian, the upper triangle is conjugated) * @param[in] rhs \c NRSMat input object to be converted * @see count, v, vec.h, NRSMat ******************************************************************************/ @@ -702,13 +702,15 @@ NRMat::NRMat(const NRSMat &rhs) { #ifdef MATPTR for (i=0; i::conjugate(rhs[k++]); } } #else for (i=0; i::conjugate(rhs[k++]); } } #endif diff --git a/nonclass.cc b/nonclass.cc index 21620dc..9f845cc 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -102,7 +102,8 @@ void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, ii=j; jj=i; } - fprintf(file, f, ((std::complex)a[ii*(ii+1)/2+jj]).real(), ((std::complex)a[ii*(ii+1)/2+jj]).imag()); + const std::complex x = (std::complex)a[(ii-1)*ii/2+jj-1]; //packed lower triangle, 1-based ii,jj + fprintf(file, f, x.real(), i>=j ? x.imag() : -x.imag()); //hermitian for complex } else fprintf(file, f, ((std::complex)a[(i-1)*c+j-1]).real(), ((std::complex)a[(i-1)*c+j-1]).imag()); if (j < n2) fputc(' ',file); } @@ -122,7 +123,8 @@ void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, ii=j; jj=i; } - fprintf(file, f, ((std::complex)a[ii*(ii+1)/2+jj]).real(), ((std::complex)a[ii*(ii+1)/2+jj]).imag()); + const std::complex x = (std::complex)a[(ii-1)*ii/2+jj-1]; //packed lower triangle, 1-based ii,jj + fprintf(file, f, x.real(), i>=j ? x.imag() : -x.imag()); //hermitian for complex } else fprintf(file,f,((std::complex)a[(i-1)*c+j-1]).real(), ((std::complex)a[(i-1)*c+j-1]).imag()); putc(j realpart > >(const NRSMat< std #ifdef CUDALA if(a.location == cpu){ #endif - cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a(0,0)), 2, result, 1); + cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a[0]), 2, result, 1); #ifdef CUDALA }else{ laerror("not implemented for cuda yet"); @@ -1573,7 +1575,7 @@ const NRSMat imagpart > >(const NRSMat< std if(a.location == cpu){ #endif - cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a(0,0))+1, 2, result, 1); + cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a[0])+1, 2, result, 1); #ifdef CUDALA }else{ laerror("not implemented for cuda yet"); diff --git a/smat.cc b/smat.cc index bbb3b9e..a948a50 100644 --- a/smat.cc +++ b/smat.cc @@ -102,7 +102,8 @@ void NRSMat::get(int fd, bool dim, bool transp, bool orcaformat) { /***************************************************************************//** - * constructor symmetrizing given matrix \f$A\f$ of general type T yielding \f$(A+A^\mathrm{T})/2\f$ + * constructor symmetrizing given matrix \f$A\f$ of general type T yielding \f$(A+A^\mathrm{T})/2\f$, + * for complex types hermitizing it, i.e. yielding \f$(A+A^\dagger)/2\f$ * @param[in] rhs matrix \f$A\f$ ******************************************************************************/ template @@ -122,7 +123,7 @@ NRSMat::NRSMat(const NRMat &rhs) { int i, j, k(0); for(i=0; i::conjugate(rhs[j][i]))/((T)2); } } } @@ -759,7 +760,7 @@ void NRSMat >::axpy(const std::complex alpha, const #ifdef CUDALA if(location == cpu){ #endif - cblas_zaxpy(nn, &alpha, x.v, 1, v, 1); + cblas_zaxpy(NN2, &alpha, x.v, 1, v, 1); #ifdef CUDALA }else{ const cuDoubleComplex _alpha = make_cuDoubleComplex(alpha.real(), alpha.imag()); diff --git a/smat.h b/smat.h index 645bd75..fdd76fe 100644 --- a/smat.h +++ b/smat.h @@ -31,6 +31,15 @@ namespace LA { template class NRPerm; template class NRSMat_from1; +//return type of the const element access NRSMat::operator(): a reference, except for complex types, +//where the upper triangle has to be conjugated on the fly +template struct NRSMat_const_element {typedef const T& type;}; +template struct NRSMat_const_element > {typedef const std::complex type;}; + +//conjugation of the stored lower triangle element when the upper triangle element is requested +template inline const T& NRSMat_conj_upper(const T &x) {return x;} +template inline const std::complex NRSMat_conj_upper(const std::complex &x) {return std::conj(x);} + /***************************************************************************//** * This class implements a general symmetric or hermitian matrix the elements @@ -154,9 +163,11 @@ public: inline const T& operator[](const size_t ij) const; inline T& operator[](const size_t ij); - //NOTE: it stores the matrix as symemtric and operator() assumes it is symmetric, does not support complex hermitean as that would require a smart pointer for l-value version - //complex conjugation options are available for BLAS routines to facilitate hermitan matrices - inline const T& operator()(const int i, const int j) const; + //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::type operator()(const int i, const int j) const; inline T& operator()(const int i, const int j); inline int nrows() const; @@ -177,7 +188,7 @@ public: { T s = NRVec(*this).sum(); if(onlytriangle) return s; - s*=2; + s += LA_traits::conjugate(s); //upper triangle is conjugated for complex (hermitian) for(int i=0; i csum() const {return rsum();}; + const NRVec csum() const //differs from rsum() by conjugation for complex (hermitian) + { + NRVec r(nn); + for(int j=0; j::NRSMat(const T *a, const int n) : nn(n), count(new int(1)) { location = DEFAULT_LOC; if(location == cpu){ #endif + v = new T[NN2]; if(LA_traits::is_plaindata()) memcpy(v, a, NN2*sizeof(T)); else for( int i=0; ij)? (i-2)*(i-1)/2+j-1 : (j-2)*(j-1)/2+i-1; /***************************************************************************//** * determine matrix element of this symmetric matrix of general type T + * PITFALL for complex (hermitian) matrices: for \f$i::operator()(const int i, const int j) { } /***************************************************************************//** - * determine matrix element of this symmetric matrix of general type T + * determine matrix element of this symmetric or hermitian matrix of general type T * @param[in] i row index running from 0 * @param[in] j column index running from 0 - * @return constant reference to the corresponding matrix element + * @return constant reference to the corresponding matrix element; for complex types the value, + * complex conjugated for \f$i::operator[] ******************************************************************************/ template -inline const T & NRSMat::operator()(const int i, const int j) const { +inline typename NRSMat_const_element::type NRSMat::operator()(const int i, const int j) const { #ifdef DEBUG if(i<0||i>=nn||j<0||j>=nn) std::cout<<"INDEX PROBLEM "<<0<<" "<::operator()(const int, const int) first index out of range - low"); @@ -805,7 +830,7 @@ inline const T & NRSMat::operator()(const int i, const int j) const { #endif NOT_GPU(*this); - return v[SMat_index(i,j)]; + return (i>=j) ? v[SMat_index_igej(i,j)] : NRSMat_conj_upper(v[SMat_index_igej(j,i)]); } /***************************************************************************//** @@ -1284,7 +1309,7 @@ std::istream& operator>>(std::istream &s, NRSMat &x) { template void NRSMat_convert(NRSMat &a, const NRSMat &b) { -a.resize(b.nrows(),b.ncols()); +a.resize(b.nrows()); for(int i=0; i &rhs): NRSMat(rhs) {}; explicit NRSMat_from1(const NRVec &rhs, const int n): NRSMat(rhs,n) {}; - inline const T& operator() (const int i, const int j) const { + inline typename NRSMat_const_element::type operator() (const int i, const int j) const { //for complex conjugated for i::operator() #ifdef DEBUG if(i<=0||i>NRSMat::nn||j<=0||j>NRSMat::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<::nn<<" "<<1<<" "<::nn<::operator() (const int, const int) first indexout of range - low"); @@ -1340,10 +1365,10 @@ public: if(j<=0) laerror("index in const T& NRSMat_from1::operator() (const int, const int) second index out of range - low"); if(j>NRSMat::nn) laerror("index in const T& NRSMat_from1::operator() (const int, const int) second index out of range - high"); #endif - return NRSMat::v[SMat_index_1(i,j)]; + return (i>=j) ? NRSMat::v[SMat_index_1igej(i,j)] : NRSMat_conj_upper(NRSMat::v[SMat_index_1igej(j,i)]); } - inline T& operator() (const int i, const int j){ + inline T& operator() (const int i, const int j){ //PITFALL: no conjugation for complex i::operator() #ifdef DEBUG if(i<=0||i>NRSMat::nn||j<=0||j>NRSMat::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<::nn<<" "<<1<<" "<::nn<::operator() (const int, const int) first indexout of range - low"); diff --git a/vec.cc b/vec.cc index fd2f79e..ede85d1 100644 --- a/vec.cc +++ b/vec.cc @@ -655,7 +655,7 @@ void NRVec >::gemv(const double beta, const NRSMat * \f[\vec{y}\leftarrow \alpha\operatorname{op}(A)\cdot\vec{x}+\beta\vec{y}\f] * @param[in] beta complex parameter \f$\beta\f$ * @param[in] A complex Hermitian matrix \f$A\f$ stored in packed form - * @param[in] trans not used + * @param[in] trans 'n' or 'c' for \f$A\f$, 't' for \f$A^\mathrm{T}=A^*\f$ * @param[in] alpha complex parameter \f$\alpha\f$ * @param[in] x complex vector \f$\vec{x}\f$ * @see gemm, NRSMat @@ -674,7 +674,10 @@ void NRVec >::gemv(const std::complex beta, #ifdef CUDALA if(location == cpu){ #endif - cblas_zhpmv(CblasRowMajor, CblasLower, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1); + if(tolower(trans) == 't') //row-major lower packed H is column-major upper packed H^T + cblas_zhpmv(CblasColMajor, CblasUpper, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1); + else + cblas_zhpmv(CblasRowMajor, CblasLower, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1); #ifdef CUDALA }else{ const cuDoubleComplex _alpha = make_cuDoubleComplex(alpha.real(), alpha.imag());