diff --git a/la_traits.h b/la_traits.h index c804d8b..6fe3abf 100644 --- a/la_traits.h +++ b/la_traits.h @@ -300,6 +300,7 @@ static void copy(std::complex *dest, std::complex *src, size_t n) {memcpy( static void clear(std::complex *dest, size_t n) {memset(dest,0,n*sizeof(std::complex));} static void copyonwrite(std::complex &x) {}; static bool is_plaindata() {return true;} +static bool is_complex() {return true;} static void clearme(std::complex &x) {x=0;}; static void deallocate(std::complex &x) {}; static inline std::complex conjugate(const std::complex &x) {return std::complex(x.real(),-x.imag());}; @@ -357,6 +358,7 @@ static void copy(C *dest, C *src, size_t n) {memcpy(dest,src,n*sizeof(C));} static void clear(C *dest, size_t n) {memset(dest,0,n*sizeof(C));} static void copyonwrite(C &x) {}; static bool is_plaindata() {return true;} +static bool is_complex() {return false;} static void clearme(C &x) {x=0;}; static void deallocate(C &x) {}; static inline C conjugate(const C &x) {return x;}; @@ -394,6 +396,7 @@ static void copy(X *dest, X *src, size_t n) {for(size_t i=0; i *dest, size_t n) {for(size_t i=0; i &x) {x.copyonwrite();}\ static bool is_plaindata() {return false;}\ +static bool is_complex() {return LA_traits::is_complex();}\ static void clearme(X &x) {x.clear();}\ static void deallocate(X &x) {x.dealloc();}\ }; @@ -436,6 +439,7 @@ static void copy(C *dest, C *src, size_t n) {for(size_t i=0; i &x) {x.copyonwrite();} \ static bool is_plaindata() {return false;}\ +static bool is_complex() {return LA_traits::is_complex();}\ static void clearme(X &x) {x.clear();} \ static void deallocate(X &x) {x.dealloc();} \ }; diff --git a/smat.h b/smat.h index 0cab729..96eb517 100644 --- a/smat.h +++ b/smat.h @@ -154,6 +154,8 @@ 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; inline T& operator()(const int i, const int j); diff --git a/t.cc b/t.cc index 72c72a2..5c3109c 100644 --- a/t.cc +++ b/t.cc @@ -4030,6 +4030,7 @@ cout <<"Error = "<<(xx-x).norm()<::is_complex()) *sign = signmult(*sign,gsign); else *sign *= gsign; if(groupindex == -1) return -1; r += groupindex * cumsizes[g]; } @@ -276,7 +293,7 @@ for(int g=0; g::is_complex()) *sign = signmult(*sign,gsign); else *sign *= gsign; if(groupindex == -1) return -1; r += groupindex * cumsizes[g]; } @@ -408,6 +425,8 @@ switch(sh->symmetry) istart= sh->offset; iend= sh->offset+sh->range-1; break; + case 2: + case -2: case 1: istart= sh->offset; if(igroup==sh->number-1) iend= sh->offset+sh->range-1; @@ -473,6 +492,8 @@ switch(sh->symmetry) istart= sh->offset; iend= sh->offset+sh->range-1; break; + case 2: + case -2: case 1: istart= sh->offset; if(igroup==sh->number-1) iend= sh->offset+sh->range-1; @@ -1135,7 +1156,7 @@ if(rhsgroup<0||rhsgroup>=rhs.shape.size()) laerror("wrong rhsgroup number in con if(rhs1.shape[group].offset != rhs.shape[rhsgroup].offset) laerror("incompatible index offset in contraction"); if(rhs1.shape[group].range != rhs.shape[rhsgroup].range) laerror("incompatible index range in contraction"); if(rhs1.shape[group].symmetry != rhs.shape[rhsgroup].symmetry) laerror("incompatible index symmetry in addgroupcontraction"); -if(rhs1.shape[group].symmetry == 1) laerror("addgroupcontraction not implemented for symmetric index groups"); +if(rhs1.shape[group].symmetry !=0 && rhs1.shape[group].symmetry != -1) laerror("addgroupcontraction only implemented for nonsymmetric and antisymmetric index groups"); #ifdef LA_TENSOR_INDEXPOSITION if(rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex == false) laerror("can contact only upper with lower index"); #endif @@ -1179,6 +1200,8 @@ if(kk!=rhsu.groupsizes[0]) laerror("internal error in addgroupcontraction"); T factor=alpha; if(u.shape[0].symmetry== -1) factor=alpha*(T)factorial(u.shape[0].number); if(u.shape[0].symmetry== 1) laerror("addgroupcontraction not implemented for symmetric index groups"); +if(u.shape[0].symmetry== 2) laerror("addgroupcontraction not implemented for hermitean index groups"); +if(u.shape[0].symmetry== -2) laerror("addgroupcontraction not implemented for antihermitean index groups"); nn=1; for(int i=1; i *>(&shape))[0]; for(int i=0; i1 ) {shape.copyonwrite(); shape[i].symmetry=1;} - if(sh[i].symmetry<-1) {shape.copyonwrite(); shape[i].symmetry= -1;} + int maxlegal = LA_traits::is_complex() ? 2 : 1; + if(sh[i].symmetry> maxlegal || sh[i].symmetry< -maxlegal) laerror("illegal index group symmetry specified"); } } diff --git a/tensor.h b/tensor.h index c8081ff..3ae3129 100644 --- a/tensor.h +++ b/tensor.h @@ -19,8 +19,9 @@ //a simple tensor class with arbitrary symmetry of index subgroups //stored in an efficient way -//indices can optionally have names and by handled by name -//each index group has a specific symmetry (nosym,sym,antisym) +//indices can optionally have names and be handled by name +//each index group has a specific symmetry (antihermitean= -2, antisym= -1, nosymmetry= 0, symmetric= 1,hermitean=2) +//NOTE: diagonal elements of antihermitean and hermitean matrices are stored including the zero imag/real part and the zeroness is NOT checked and similarly for higher rank tensors //additional symmetry between index groups (like in 2-electron integrals) is not supported directly, you would need to nest the class to Tensor > //leftmost index is least significant (changing fastest) in the storage order //presently only a rudimentary implementation @@ -48,6 +49,45 @@ namespace LA { +template +inline T signeddata(const int sgn, const T data, const bool lhs=false) +{ +if(LA_traits::is_complex()) //condition known at compile time + { + switch(sgn) + { + case 2: + return LA_traits::conjugate(data); + break; + case 1: + return data; + break; + case -1: + return -data; + break; + case -2: + return -LA_traits::conjugate(data); + break; + case 0: +#ifdef DEBUG + if(lhs) laerror("dereferencing lhs Signedpointer to nonexistent tensor element"); +#endif + return 0; + break; + } + return 0; + } + else // for real + { + if(sgn>0) return data; + if(sgn<0) return -data; +#ifdef DEBUG + if(sgn==0 && lhs) laerror("dereferencing lhs Signedpointer to nonexistent tensor element"); +#endif + return 0; + } + +} template class Signedpointer @@ -57,19 +97,11 @@ int sgn; public: Signedpointer(T *p, int s) : ptr(p),sgn(s) {}; //dereferencing *ptr should be ignored for sgn==0 - const T operator=(const T rhs) - { - if(sgn>0) *ptr = rhs; - if(sgn<0) *ptr = -rhs; -#ifdef DEBUG - if(sgn==0) laerror("dereferencing lhs Signedpointer to nonexistent tensor element"); -#endif - return rhs; - } - T& operator*=(const T rhs) {*ptr *= rhs; return *ptr;} - T& operator/=(const T rhs) {*ptr /= rhs; return *ptr;} - T& operator+=(const T rhs) {if(sgn>0) *ptr += rhs; else *ptr -= rhs; return *ptr;} - T& operator-=(const T rhs) {if(sgn>0) *ptr -= rhs; else *ptr += rhs; return *ptr;} + const T operator=(const T rhs) {*ptr = signeddata(sgn,rhs); return rhs;} + void operator*=(const T rhs) {*ptr *= rhs;} + void operator/=(const T rhs) {*ptr /= rhs;} + void operator+=(T rhs) {*ptr += signeddata(sgn,rhs);} + void operator-=(T rhs) {*ptr -= signeddata(sgn,rhs);} }; @@ -104,7 +136,7 @@ class LA_traits { typedef class INDEXGROUP { public: int number; //number of indices -int symmetry; //-1 0 or 1, later 2 for hermitian and -2 for antihermitian? - would need change in operator() and Signedpointer +int symmetry; //-1 0 or 1, later 2 for hermitian and -2 for antihermitian #ifdef LA_TENSOR_ZERO_OFFSET static const LA_index offset = 0; //compiler can optimize away some computations #else @@ -229,6 +261,7 @@ public: bool is_flat() const {for(int i=0; i1) return false; return true;}; bool is_compressed() const {for(int i=0; i1&&shape[i].symmetry!=0) return true; return false;}; bool has_symmetry() const {for(int i=0; i::is_complex()) return false; for(int i=0; i 1) return true; return false;}; void clear() {data.clear();}; void defaultnames(const char *basename="i") {names.resize(rank()); for(int i=0; i &s) {shape=s; data.resize(calcsize()); calcrank(); names.clear();}; void deallocate() {data.resize(0); shape.resize(0); groupsizes.resize(0); cumsizes.resize(0); names.resize(0);}; + inline Signedpointer lhs(const SUPERINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer(&data[i],sign);}; - inline T operator()(const SUPERINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); if(sign==0) return 0; return sign>0 ?data[i] : -data[i];}; + inline T operator()(const SUPERINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); return signeddata(sign,data[i]);}; inline Signedpointer lhs(const FLATINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer(&data[i],sign);}; - inline T operator()(const FLATINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); if(sign==0) return 0; return sign>0 ?data[i] : -data[i];}; + inline T operator()(const FLATINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); return signeddata(sign,data[i]);}; inline Signedpointer lhs(LA_index i1...) {va_list args; int sign; LA_largeindex i; va_start(args,i1); i= vindex(&sign, i1,args); return Signedpointer(&data[i],sign); }; - inline T operator()(LA_index i1...) const {va_list args; ; int sign; LA_largeindex i; va_start(args,i1); i= vindex(&sign, i1,args); if(sign==0) return 0; return sign>0 ?data[i] : -data[i];}; + inline T operator()(LA_index i1...) const {va_list args; ; int sign; LA_largeindex i; va_start(args,i1); i= vindex(&sign, i1,args); return signeddata(sign,data[i]);}; inline Tensor& operator=(const Tensor &rhs) {myrank=rhs.myrank; shape=rhs.shape; groupsizes=rhs.groupsizes; cumsizes=rhs.cumsizes; data=rhs.data; names=rhs.names; return *this;}; @@ -298,7 +332,7 @@ public: void put(int fd, bool with_names=false) const; void get(int fd, bool with_names=false); - inline void randomize(const typename LA_traits::normtype &x) {data.randomize(x);}; + inline void randomize(const typename LA_traits::normtype &x) {if(has_hermiticity()) laerror("randomization does not support correct treatment of hermitean/antihermitean index groups"); data.randomize(x);}; void loopover(void (*callback)(const SUPERINDEX &, T *)); //loop over all elements void constloopover(void (*callback)(const SUPERINDEX &, const T *)) const; //loop over all elements