tensor: bugfixes and hermiticity enforcement

la_traits.h

@@ -306,6 +306,8 @@ 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());};
 static inline C realpart(const std::complex<C> &x) {return x.real();}
 static inline C imagpart(const std::complex<C> &x) {return x.imag();}
+static inline void setrealpart(std::complex<C> &x, const C &y) {reinterpret_cast<C(&)[2]>(x)[0]=y;}
+static inline void setimagpart(std::complex<C> &x, const C &y) {reinterpret_cast<C(&)[2]>(x)[1]=y;}
 };
 
 
@@ -364,6 +366,8 @@ static void deallocate(C &x) {};
 static inline C conjugate(const C &x) {return x;};
 static inline C realpart(const C &x) {return x;}
 static inline C imagpart(const C &x) {return 0;}
+static inline void setrealpart(C &x, const C &y) {x=y;}
+static inline void setimagpart(C &x, const C &y) {}
 };
 
 

t.cc

@@ -4263,7 +4263,7 @@ for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int
 cout <<"Error = "<<(z-zz).norm()<<endl;
 }
 
-if(1)
+if(0)
 {
 int nn;
 cin>>nn;
@@ -4309,5 +4309,88 @@ cout <<"Error = "<<(z-zzz).norm()<<endl;
 
 }
 
+if(0)
+{
+//check symmetrizer/antisymmetrizer in general case
+int r,n,sym;
+cin>>r>>n>>sym;
+NRVec<INDEXGROUP> shape(3);
+	shape[0].number=2; 
+        shape[0].symmetry=0;
+        shape[0].range=n+1;
+        shape[0].offset=0;
+
+        shape[1].number=r;
+        shape[1].symmetry= sym;
+        shape[1].range=n;
+        shape[1].offset=0;
+
+        shape[2].number=2;
+        shape[2].symmetry=0;
+        shape[2].range=n+2;
+        shape[2].offset=0;
+
+
+Tensor<complex<double> > x(shape); x.randomize(1.);
+x.defaultnames();
+cout <<"x= "<<x.shape << " "<<x.names<<endl;
+Tensor<complex<double> > xf=x.flatten(1);
+cout <<"xf= "<<xf.shape << " "<<xf.names<<endl;
+Tensor<complex<double> > xxx=x.unwind_index_group(1);
+cout <<"xxx= "<<xxx.shape<<" "<<xxx.names<<endl;
+
+INDEXLIST il(r);
+for(int i=0; i<r; ++i) il[i]= {1+i,0};
+Tensor<complex<double> > xx = xf.merge_indices(il,sym);
+cout <<"xx = "<<xx.shape<< " "<<xx.names<<endl;
+cout <<"Error = "<<(xx-xxx).norm()<<endl;
+}
+
+if(1)
+{
+int nn;
+cin>>nn;
+int r=4;
+NRVec<INDEXGROUP> s(r);
+for(int i=0; i<r; ++i)
+        {
+        s[i].number=1;
+        s[i].symmetry=0;
+        s[i].range=nn;
+        s[i].offset=0;
+        }
+Tensor<complex<double> > x(s); x.randomize(1.);
+INDEXNAME xlist[] = {"i","j","k","l"};
+x.names=NRVec<INDEXNAME>(xlist);
+Tensor<complex<double> > y(s); y.randomize(1.);
+INDEXNAME ylist[] = {"n","m","j","i"};
+y.names=NRVec<INDEXNAME>(ylist);
+
+Tensor<complex<double> >z(s); z.clear();
+INDEXNAME zlist[] = {"k","l","m","n"};
+z.names=NRVec<INDEXNAME>(zlist);
+Tensor<complex<double> >zz(z);
+
+
+z.add_permuted_contractions("k/l,m",x,y,1,0,false,false);
+
+zz.copyonwrite();
+for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int n=0; n<nn; ++n)
+	{
+	complex<double>  s=0;
+	for(int i=0; i<nn; ++i) for(int j=0; j<nn; ++j) s += x(i,j,k,l)*y(n,m,j,i);
+	zz.lhs(k,l,m,n) = s;
+	}
+
+Tensor<complex<double> > zzz(s);
+for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int n=0; n<nn; ++n)
+	{
+	zzz.lhs(k,l,m,n) = (zz(k,l,m,n)-zz(l,k,m,n)+zz(m,k,l,n));
+	}
+
+cout <<"Error = "<<(z-zzz).norm()<<endl;
+
+}
+
 
 }//main

tensor.cc

@@ -24,6 +24,7 @@
 #include <complex>
 #include "bitvector.h"
 #include <string.h>
+#include <bits/stdc++.h>
 
 
 namespace LA {
@@ -899,7 +900,7 @@ template<typename T>
 static void flatten_callback(const SUPERINDEX &I, T *v)
 {
 FLATINDEX J = superindex2flat(I);
-//std::cout <<"TEST flatten_callback: from "<<JP<<" TO "<<J<<std::endl;
+//std::cout <<"TEST flatten_callback: "<<J<<std::endl;
 *v = (*help_tt<T>)(J); //rhs operator() generates the redundant elements for the unwinded lhs tensor
 }       
 //
@@ -1432,6 +1433,7 @@ template<typename T>
 static const PermutationAlgebra<int,T> *help_pa;
 
 static bool help_inverse;
+static bool help_conjugate;
 
 template<typename T>
 static T help_alpha;
@@ -1443,7 +1445,10 @@ FLATINDEX J = superindex2flat(I);
 for(int p=0; p<help_pa<T>->size(); ++p)
 	{
 	FLATINDEX Jp = J.permuted((*help_pa<T>)[p].perm,help_inverse);
-	*v += help_alpha<T> * (*help_pa<T>)[p].weight * (*help_tt<T>)(Jp);
+	bool conj=false;
+	if(help_conjugate) conj= ((*help_pa<T>)[p].perm).parity()<0;
+	T tmp= (*help_tt<T>)(Jp);
+	*v += help_alpha<T> * (*help_pa<T>)[p].weight * (conj ? LA_traits<T>::conjugate(tmp) : tmp);
 	}
 }
 
@@ -1476,7 +1481,7 @@ for(int p=0; p<help_pa<T>->size(); ++p)
 
 
 template<typename T>
-void Tensor<T>::apply_permutation_algebra(const Tensor<T> &rhs, const PermutationAlgebra<int,T> &pa, bool inverse, T alpha, T beta)
+void Tensor<T>::apply_permutation_algebra(const Tensor<T> &rhs, const PermutationAlgebra<int,T> &pa, bool inverse, T alpha, T beta, bool conjugate)
 {
 if(beta!=(T)0) {if(beta!=(T)1) *this *= beta;} else clear();
 if(alpha==(T)0) return;
@@ -1500,6 +1505,7 @@ if(rhs.is_named())
 help_tt<T> = &rhs;
 help_pa<T> = &pa;
 help_inverse = inverse;
+help_conjugate= LA_traits<T>::is_complex() ? conjugate : false;
 help_alpha<T> = alpha;
 loopover(permutationalgebra_callback);
 }
@@ -1882,6 +1888,7 @@ if(!basicperm.is_valid()) laerror("internal error in merge_indices");
 
 //prepare permutation algebra
 PermutationAlgebra<int,T> pa;
+bool doconjugate=false;
 if(sym==0)
 	{
 	pa.resize(1);
@@ -1890,6 +1897,7 @@ if(sym==0)
 	}
 else
 	{
+	doconjugate = sym>=2||sym<= -2;
 	PermutationAlgebra<int,int> sa = sym>0 ? symmetrizer<int>(il.size()) : antisymmetrizer<int>(il.size());
 	//std::cout <<"SA = "<<sa<<std::endl;
 	pa.resize(sa.size());
@@ -1905,7 +1913,7 @@ else
 //std::cout <<"Use PA = "<<pa<<std::endl;
 
 Tensor<T> r(newshape);
-r.apply_permutation_algebra(*this,pa,false,(T)1/(T)pa.size(),0);
+r.apply_permutation_algebra(*this,pa,false,(T)1/(T)pa.size(),0,doconjugate);
 if(is_named()) r.names=names.permuted(basicperm,true);	
 return r;
 }
@@ -2200,6 +2208,82 @@ if(is_named())
 }
 
 
+static void checkindex(const NRVec<INDEXGROUP> &shape, const  SUPERINDEX &I, bool &zeroreal, bool &zeroimag)
+{
+#ifdef DEBUG
+if(shape.size()!=I.size()) laerror("inconsistent shape and index in checkindex");
+#endif
+zeroreal=zeroimag=false;
+for(int g=0; g<I.size(); ++g)
+	{
+#ifdef DEBUG
+	if(I[g].size()!=shape[g].number) laerror("inconsistent2 shape and index in checkindex");
+#endif
+	if(shape[g].symmetry == 2 || shape[g].symmetry == -2)
+		{
+		bool sameindex=false;
+		for(int i=1; i<shape[g].number; ++i) for(int j=0; j<i; ++j)
+			if(I[g][i]==I[g][j]) 
+				{
+				sameindex=true;
+				goto checkfailed; //for this group, but do the remaining ones for further restrictions
+				}
+		
+		checkfailed:
+		if(sameindex)
+			{
+			if(shape[g].symmetry == 2) zeroimag=true;
+			if(shape[g].symmetry == -2) zeroreal=true;
+			}
+		}
+	}
+}
+
+
+template<typename T>
+static void hermiticity_callback(const SUPERINDEX &I, T *v)
+{
+bool zeroimag;
+bool zeroreal;
+checkindex(help_t<T>->shape,I,zeroreal,zeroimag);
+if(zeroreal) LA_traits<T>::setrealpart(*v,0);
+if(zeroimag) LA_traits<T>::setimagpart(*v,0);
+}
+
+
+template<typename T>
+static void hermiticity_callback2(const SUPERINDEX &I, const T *v)
+{
+bool zeroimag;
+bool zeroreal;
+checkindex(help_tt<T>->shape,I,zeroreal,zeroimag);
+if(zeroreal && LA_traits<T>::realpart(*v)!=(T)0 || zeroimag  &&  LA_traits<T>::imagpart(*v)!=(T)0) throw true;
+}
+
+
+
+template<typename T>
+void Tensor<T>::enforce_hermiticity()
+{
+if(!has_hermiticity()) return;
+help_t<T> = this;
+loopover(hermiticity_callback);
+}
+
+
+template<typename T>
+bool Tensor<T>::fulfills_hermiticity() const
+{
+if(!has_hermiticity()) return true;
+help_tt<T> = this;
+try	{
+	constloopover(hermiticity_callback2);
+	}
+catch(bool failed) {return false;}
+return true;
+}
+
+
 
 template class Tensor<double>;
 template class Tensor<std::complex<double> >;

tensor.h

@@ -50,7 +50,7 @@
 namespace LA {
 
 template<typename T>
-inline T signeddata(const int sgn, const T data, const bool lhs=false)
+static inline T signeddata(const int sgn, const T &data)
 {
 if(LA_traits<T>::is_complex()) //condition known at compile time
                                 {
@@ -69,9 +69,6 @@ if(LA_traits<T>::is_complex()) //condition known at compile time
                                                 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;
                                         }
@@ -81,9 +78,6 @@ if(LA_traits<T>::is_complex()) //condition known at compile time
                                 {
                                 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;
                                 }
 
@@ -273,12 +267,24 @@ public:
 	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 Signedpointer<T> lhs(const SUPERINDEX &I) {int sign; LA_largeindex i=index(&sign,I);
-	inline T operator()(const SUPERINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); return signeddata(sign,data[i]);};
+#ifdef DEBUG
-        inline Signedpointer<T> lhs(const FLATINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer<T>(&data[i],sign);};
+				if(sign==0) laerror("l-value pointer to nonexistent tensor element"); 
-        inline T operator()(const FLATINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); return signeddata(sign,data[i]);};
+#endif
-        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); };
+				 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); return signeddata(sign,data[i]);};
+	inline T operator()(const SUPERINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); if(sign==0) return 0; else return signeddata(sign,data[i]);};
+        inline Signedpointer<T> lhs(const FLATINDEX &I) {int sign; LA_largeindex i=index(&sign,I); 
+#ifdef DEBUG
+                                if(sign==0) laerror("l-value pointer to nonexistent tensor element");
+#endif
+				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; else 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);  
+#ifdef DEBUG
+                                if(sign==0) laerror("l-value pointer to nonexistent tensor element");
+#endif
+				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; else 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;};
 
@@ -332,7 +338,9 @@ 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) {if(has_hermiticity()) laerror("randomization does not support correct treatment of hermitean/antihermitean index groups"); data.randomize(x);};
+	void enforce_hermiticity(); //zero out real/imag parts for repeated indices appropriately
+	bool fulfills_hermiticity() const; //check it is so
+	inline void randomize(const typename LA_traits<T>::normtype &x) {data.randomize(x); enforce_hermiticity();};
 
 	void loopover(void (*callback)(const SUPERINDEX &, T *)); //loop over all elements
 	void constloopover(void (*callback)(const SUPERINDEX &, const T *)) const; //loop over all elements
@@ -373,7 +381,7 @@ public:
 	Tensor innercontraction(const NRVec<INDEXNAME> &nl1, const NRVec<INDEXNAME> &nl2) const {return innercontraction(findindexlist(nl1),findindexlist(nl2));};
 	Tensor innercontraction(const INDEXNAME &n1, const INDEXNAME &n2) const {return innercontraction(findindex(n1),findindex(n2));};
 
-	void apply_permutation_algebra(const  Tensor &rhs, const PermutationAlgebra<int,T> &pa, bool inverse=false, T alpha=1, T beta=0); //general (not optimally efficient) symmetrizers, antisymmetrizers etc. acting on the flattened index list:
+	void apply_permutation_algebra(const  Tensor &rhs, const PermutationAlgebra<int,T> &pa, bool inverse=false, T alpha=1, T beta=0, bool conjugate_by_parity=false); //general (not optimally efficient) symmetrizers, antisymmetrizers etc. acting on the flattened index list:
 	void apply_permutation_algebra(const  NRVec<Tensor> &rhsvec, const PermutationAlgebra<int,T> &pa, bool inverse=false, T alpha=1, T beta=0); //avoids explicit outer product but not vectorized, rather inefficient
 //								 this *=beta; for I over this: this(I) += alpha * sum_P c_P rhs(P(I))
 //								 PermutationAlgebra can represent e.g. general_antisymmetrizer in Kucharski-Bartlett notation, or Grassmann products building RDM from cumulants