jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

tensor: support for complex (anti)hermitian tensors

Browse Source
This commit is contained in:
Jiri Pittner
2025-11-18 17:30:58 +01:00
parent 417a7d1d1a
commit 20a61e2fb9
5 changed files with 101 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
la_traits.h
View File

@@ -300,6 +300,7 @@ static void copy(std::complex<C> *dest, std::complex<C> *src, size_t n) {memcpy(
static void clear(std::complex<C> *dest, size_t n) {memset(dest,0,n*sizeof(std::complex<C>));}
static void copyonwrite(std::complex<C> &x) {};
static bool is_plaindata() {return true;}
static bool is_complex() {return true;}
static void clearme(std::complex<C> &x) {x=0;};
static void deallocate(std::complex<C> &x) {};
static inline std::complex<C> conjugate(const std::complex<C> &x) {return std::complex<C>(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<C> *dest, X<C> *src, size_t n) {for(size_t i=0; i<n; ++i) des
static void clear(X<C> *dest, size_t n) {for(size_t i=0; i<n; ++i) dest[i].clear();}\
static void copyonwrite(X<C> &x) {x.copyonwrite();}\
static bool is_plaindata() {return false;}\
static bool is_complex() {return LA_traits<C>::is_complex();}\
static void clearme(X<C> &x) {x.clear();}\
static void deallocate(X<C> &x) {x.dealloc();}\
};
@@ -436,6 +439,7 @@ static void copy(C *dest, C *src, size_t n) {for(size_t i=0; i<n; ++i) dest[i]=s
static void clear(C *dest, size_t n) {for(size_t i=0; i<n; ++i) dest[i].clear();} \
static void copyonwrite(X<C> &x) {x.copyonwrite();} \
static bool is_plaindata() {return false;}\
static bool is_complex() {return LA_traits<C>::is_complex();}\
static void clearme(X<C> &x) {x.clear();} \
static void deallocate(X<C> &x) {x.dealloc();} \
};

2
smat.h
View File

@@ -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);

2
t.cc
View File

@@ -4030,6 +4030,7 @@ cout <<"Error = "<<(xx-x).norm()<<endl;
}
if(0)
{
//check symmetrizer/antisymmetrizer in general case
@@ -4067,6 +4068,7 @@ cout <<"xx = "<<xx.shape<< " "<<xx.names<<endl;
cout <<"Error = "<<(xx-xxx).norm()<<endl;
}
if(0)
{
int r=4;

55
tensor.cc
View File

@@ -62,6 +62,8 @@ for(int i=0; i<shape.size(); ++i)
case 0:
s *= groupsizes[i] = longpow(sh->range,sh->number);
break;
case 2:
case -2:
case 1:
s *= groupsizes[i] = simplicial(sh->number,sh->range);
break;
@@ -104,18 +106,17 @@ switch(I.size()) //a few special cases for efficiency
break;
case 2:
{
*sign=1;
if(g.symmetry==0) return (I[1]-g.offset)*g.range+I[0]-g.offset;
if(g.symmetry==0) {*sign=1; return (I[1]-g.offset)*g.range+I[0]-g.offset;};
LA_index i0,i1;
if(I[0]>I[1]) {i1=I[0]; i0=I[1]; if(g.symmetry<0) *sign = -1;} else {i1=I[1]; i0=I[0];}
if(I[0]>I[1]) {i1=I[0]; i0=I[1]; *sign=g.symmetry;} else {i1=I[1]; i0=I[0]; *sign=1;}
i0 -= g.offset;
i1 -= g.offset;
if(g.symmetry<0)
if(g.symmetry == -1) //antisymmetric
{
if(i0==i1) {*sign=0; return -1;}
return i1*(i1-1)/2+i0;
}
else
else //symmetric, hermitian, antihermitian
{
return i1*(i1+1)/2+i0;
}
@@ -124,10 +125,9 @@ switch(I.size()) //a few special cases for efficiency
default: //general case
{
*sign=1;
if(g.symmetry==0) //rectangular case
{
*sign=1;
LA_largeindex r=0;
for(int i=I.size()-1; i>=0; --i)
{
@@ -143,8 +143,8 @@ switch(I.size()) //a few special cases for efficiency
II.copyonwrite();
if(g.offset!=0) II -= g.offset;
int parity=netsort(II.size(),&II[0]);
if(g.symmetry<0 && (parity&1)) *sign= -1;
if(g.symmetry<0) //antisymmetric
*sign= (parity&1) ? g.symmetry : 1;
if(g.symmetry == -1) //antisymmetric - do not store zero diagonal
{
for(int i=0; i<I.size()-1; ++i)
if(II[i]==II[i+1])
@@ -154,7 +154,7 @@ switch(I.size()) //a few special cases for efficiency
for(int i=0; i<II.size(); ++i) r += simplicial(i+1,II[i]-i);
return r;
}
else //symmetric
else //symmetric, hermitian, antihermitian
{
LA_largeindex r=0;
for(int i=0; i<II.size(); ++i) r += simplicial(i+1,II[i]);
@@ -181,6 +181,8 @@ switch(g.symmetry)
s /= g.range;
}
break;
case 2:
case -2:
case 1:
for(int i=g.number; i>0; --i)
{
@@ -221,6 +223,21 @@ for(int g=shape.size()-1; g>=0; --g)
return I;
}
//group-like multiplication table to combine symmetry adjustments due to several index groups
static const int signmultab[5][5] = {
{1,2,0,-2,-1},
{2,1,0,-1,-2},
{0,0,0,0,0},
{-2,-1,0,1,2},
{-1,-2,0,2,1}
};
static inline int signmult(int s1, int s2)
{
return signmultab[s1+2][s2+2];
}
template<typename T>
@@ -250,7 +267,7 @@ for(int g=0; g<shape.size(); ++g) //loop over index groups
int gsign;
LA_largeindex groupindex = subindex(&gsign,shape[g],I[g]);
//std::cout <<"INDEX TEST group "<<g<<" cumsizes "<< cumsizes[g]<<" groupindex "<<groupindex<<std::endl;
*sign *= gsign;
if(LA_traits<T>::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<shape.size(); ++g) //loop over index groups
gstart=gend+1;
LA_largeindex groupindex = subindex(&gsign,shape[g],subI);
//std::cout <<"FLATINDEX TEST group "<<g<<" cumsizes "<< cumsizes[g]<<" groupindex "<<groupindex<<std::endl;
*sign *= gsign;
if(LA_traits<T>::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<u.shape.size(); ++i) nn*= u.groupsizes[i];
mm=1; for(int i=1; i<rhsu.shape.size(); ++i) mm*= rhsu.groupsizes[i];
data.copyonwrite();
@@ -1645,7 +1668,7 @@ if(is_named() && rhs.is_named() && names!=rhs.names) laerror("incompatible tenso
T factor=1;
for(int i=0; i<shape.size(); ++i)
{
if(shape[i].symmetry==1) laerror("unsupported index group symmetry in dot");
if(shape[i].symmetry==1||shape[i].symmetry==2||shape[i].symmetry== -2) laerror("unsupported index group symmetry in dot");
if(shape[i].symmetry== -1) factor *= (T)factorial(shape[i].number);
}
return factor * data.dot(rhs.data);
@@ -1897,8 +1920,8 @@ const INDEXGROUP *sh = &(* const_cast<const NRVec<INDEXGROUP> *>(&shape))[0];
for(int i=0; i<shape.size(); ++i)
{
if(sh[i].number==1 && sh[i].symmetry!=0) {shape.copyonwrite(); shape[i].symmetry=0;}
if(sh[i].symmetry>1 ) {shape.copyonwrite(); shape[i].symmetry=1;}
if(sh[i].symmetry<-1) {shape.copyonwrite(); shape[i].symmetry= -1;}
int maxlegal = LA_traits<T>::is_complex() ? 2 : 1;
if(sh[i].symmetry> maxlegal || sh[i].symmetry< -maxlegal) laerror("illegal index group symmetry specified");
}
}

74
tensor.h
View File

@@ -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<Tensor<T> >
//leftmost index is least significant (changing fastest) in the storage order
//presently only a rudimentary implementation
@@ -48,6 +49,45 @@
namespace LA {
template<typename T>
inline T signeddata(const int sgn, const T data, const bool lhs=false)
{
if(LA_traits<T>::is_complex()) //condition known at compile time
{
switch(sgn)
{
case 2:
return LA_traits<T>::conjugate(data);
break;
case 1:
return data;
break;
case -1:
return -data;
break;
case -2:
return -LA_traits<T>::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<typename T>
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<INDEXNAME> {
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; i<shape.size(); ++i) if(shape[i].number>1) return false; return true;};
bool is_compressed() const {for(int i=0; i<shape.size(); ++i) if(shape[i].number>1&&shape[i].symmetry!=0) return true; return false;};
bool has_symmetry() const {for(int i=0; i<shape.size(); ++i) if(shape[i].symmetry!=0) return true; return false;};
bool has_hermiticity() const {if(!LA_traits<T>::is_complex()) return false; for(int i=0; i<shape.size(); ++i) if(shape[i].symmetry < -1 || shape[i].symmetry > 1) return true; return false;};
void clear() {data.clear();};
void defaultnames(const char *basename="i") {names.resize(rank()); for(int i=0; i<rank(); ++i) sprintf(names[i].name,"%s%03d",basename,i);}
int rank() const {return myrank;};
@@ -239,12 +272,13 @@ public:
void copyonwrite() {shape.copyonwrite(); groupsizes.copyonwrite(); cumsizes.copyonwrite(); data.copyonwrite(); names.copyonwrite();};
void resize(const NRVec<INDEXGROUP> &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<T> lhs(const SUPERINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer<T>(&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<T> lhs(const FLATINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer<T>(&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<T> lhs(LA_index i1...) {va_list args; int sign; LA_largeindex i; va_start(args,i1); i= vindex(&sign, i1,args); return Signedpointer<T>(&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<T>::normtype &x) {data.randomize(x);};
inline void randomize(const typename LA_traits<T>::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