From 3f7586378db76f15168e2238e443b7b974654836 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Fri, 14 Nov 2025 16:16:22 +0100
Subject: [PATCH] working on add_permuted_contractions

---
 t.cc      | 49 ++++++++++++++++++++++++++++++++++-
 tensor.cc | 76 ++++++++++++++++++++++++++++++++++++++++++++++++++++---
 tensor.h  |  9 ++++---
 3 files changed, 127 insertions(+), 7 deletions(-)

diff --git a/t.cc b/t.cc
index 78e50db..d31c91b 100644
--- a/t.cc
+++ b/t.cc
@@ -3568,7 +3568,7 @@ v.printsorted(cout,1,false);
 }
 
 
-if(1)
+if(0)
 {
 //grassmann product of n identical rank=2 tensors in m-dim space
 int n,m;
@@ -4164,5 +4164,52 @@ cout<<"resulting index order = "<<y.names<<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<double> x(s); x.randomize(1.);
+INDEXNAME xlist[] = {"i","j","k","l"};
+x.names=NRVec<INDEXNAME>(xlist);
+Tensor<double> y(s); y.randomize(1.);
+INDEXNAME ylist[] = {"n","m","j","i"};
+y.names=NRVec<INDEXNAME>(ylist);
+
+Tensor<double>z(s); z.clear();
+INDEXNAME zlist[] = {"k","l","m","n"};
+z.names=NRVec<INDEXNAME>(zlist);
+Tensor<double>zz(z);
+
+INDEX i1[]={{0,0},{1,0}};
+INDEX i2[]={{3,0},{2,0}};
+NRVec<INDEX> il1(i1);
+NRVec<INDEX> il2(i2);
+Tensor<double> zzz = x.contractions(il1,y,il2,1,false,false);
+//cout <<"zzz names = "<<zzz.names<<endl;
+
+z.add_permuted_contractions("",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)
+	{
+	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;
+	//cout <<"test "<<k<<" "<<l<<" "<<m<<" "<<n<<" "<<zz(k,l,m,n)<<" "<<zzz(n,m,k,l)<< " : "<<z(k,l,m,n) <<endl;
+	}
+
+cout <<"Error = "<<(z-zz).norm()<<endl;
+
+}
+
 
 }//main
diff --git a/tensor.cc b/tensor.cc
index 588c709..12280a7 100644
--- a/tensor.cc
+++ b/tensor.cc
@@ -1198,7 +1198,7 @@ if(rhsindex<0||rhsindex>=rhs.shape[rhsgroup].number)  laerror("wrong index numbe
 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");
 #ifdef LA_TENSOR_INDEXPOSITION
-if(rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex == false) laerror("can contact only upper with lower index");
+if(rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex == false) laerror("can contract only upper with lower index");
 #endif
 
 const Tensor<T> u = conjugate1? (rhs1.unwind_index(group,index)).conjugate() : rhs1.unwind_index(group,index);
@@ -1458,9 +1458,11 @@ if(rank()!=rhs.rank()) laerror("rank mismatch in apply_permutation_algebra");
 
 if(rhs.is_named())
 	{
-	NRVec<INDEXNAME> namperm = rhs.names.permuted(pa[0].perm,inverse);
+	NRVec<INDEXNAME> namperm = rhs.names.permuted(pa[0].perm,!inverse);
 	if(is_named())
 		{
+		std::cout <<"LHS names = "<<names<<std::endl;
+		std::cout <<"permuted RHS names = "<<namperm<<std::endl;
 		if(names!=namperm) laerror("inconsistent index names in apply_permutation_algebra");
 		}
 	else	
@@ -1492,7 +1494,7 @@ if(allnamed)
 	{
 	NRVec<INDEXNAME> allrhsnames=rhsvec[0].names;
 	for(int i=1; i<rhsvec.size(); ++i) allrhsnames.concatme(rhsvec[i].names);
-	NRVec<INDEXNAME> namperm = allrhsnames.permuted(pa[0].perm,inverse);
+	NRVec<INDEXNAME> namperm = allrhsnames.permuted(pa[0].perm,!inverse);
         if(is_named())
                 {
                 if(names!=namperm) laerror("inconsistent index names in apply_permutation_algebra");
@@ -1905,6 +1907,74 @@ std::istream & operator>>(std::istream &s, INDEX &x)
 s>>x.group>>x.index;
 return s;
 }
+
+
+template<typename T>
+void Tensor<T>::add_permuted_contractions(const char *antisymmetrizer, const Tensor &rhs1, const Tensor &rhs2, T alpha, T beta, bool conjugate1, bool conjugate2)
+{
+if(!rhs1.is_uniquely_named()||!rhs2.is_uniquely_named()|| !is_uniquely_named()) laerror("tensors must have unique named indices in add_permuted_contractions");
+
+//find contraction indices
+int nc=0;
+for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs1.names.size(); ++j) 
+	if(rhs1.names[i]==rhs2.names[j])
+		{
+		//std::cout << "found contraction "<<nc<<" th. index = "<<rhs1.names[i]<<std::endl;
+#ifdef LA_TENSOR_INDEXPOSITION
+		int group1=indexposition(i,rhs1.shape);
+		int group2=indexposition(j,rhs2.shape);
+		if(rhs1.shape[group1].upperindex ^ rhs2.shape[group2].upperindex == false) laerror("can contract only upper with lower index");
+#endif
+		++nc;
+		}
+
+INDEXLIST il1(nc),il2(nc);
+int ii=0;
+for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs1.names.size(); ++j)
+        if(rhs1.names[i]==rhs2.names[j])
+		{
+		il1[ii] = indexposition(i,rhs1.shape);
+		il2[ii] = indexposition(j,rhs2.shape);
+		++ii;
+		}
+
+//std::cout<<"contraction list1 = "<<il1<<std::endl;
+//std::cout<<"contraction list2 = "<<il2<<std::endl;
+
+Tensor<T> tmp=rhs1.contractions(il1,rhs2,il2,alpha,conjugate1,conjugate2);
+if(rank()!=tmp.rank()) laerror("rank mismatch in add_permuted_contractions");
+
+//generate the antisymmetrizer, adding also indices not involved as a constant subpermutation
+//@@@
+PermutationAlgebra<int,T> pa(1);
+pa[0].weight=1;
+pa[0].perm.resize(rank());
+pa[0].perm.identity();
+
+//std::cout <<"LHS names = "<<names<<std::endl;
+//std::cout <<"TMP names = "<<tmp.names<<std::endl;
+
+//find permutation which will bring indices of tmp to the order as in *this: this->names[i] == tmp.names[p[i]]
+NRPerm<int> basicperm(rank());
+basicperm.clear();
+for(int i=0; i<rank(); ++i) 
+	{
+	for(int j=0; j<tmp.rank(); ++j)
+             if((* const_cast<const NRVec<INDEXNAME> *>(&names))[i] == tmp.names[j])
+		{
+		basicperm[i+1]=j+1;
+		break;
+		}
+	}
+if(!basicperm.is_valid()) laerror("indices mismatch between lhs and rhs in add_permuted_contractions");
+std::cout <<"Basic permutation = "<<basicperm<<std::endl;
+
+PermutationAlgebra<int,T> pb = basicperm*pa;
+apply_permutation_algebra(tmp,pb,false,(T)1/(T)pb.size(),beta);
+//equivalently possible 
+//PermutationAlgebra<int,T> pb = basicperm.inverse()*pa;
+//apply_permutation_algebra(tmp,pb,true,(T)1/(T)pb.size(),beta);
+}
         
 
 
diff --git a/tensor.h b/tensor.h
index e5792ec..4056cb0 100644
--- a/tensor.h
+++ b/tensor.h
@@ -230,7 +230,7 @@ public:
 	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;};
 	void clear() {data.clear();};
-	void defaultnames() {names.resize(rank()); for(int i=0; i<rank(); ++i) sprintf(names[i].name,"i%03d",i);}
+	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;};
 	int calcrank(); //is computed from shape
 	void canonicalize_shape(); 
@@ -322,8 +322,11 @@ public:
 	inline Tensor contractions( const INDEXLIST &il1, const Tensor &rhs2, const INDEXLIST &il2,  T alpha=1, bool conjugate1=false, bool conjugate2=false) const {Tensor<T> r; r.addcontractions(*this,il1,rhs2,il2,alpha,0,true,conjugate1, conjugate2); return r; };
 	inline Tensor contractions(const Tensor &rhs2, const NRVec<INDEXNAME> &names,  T alpha=1, bool conjugate1=false, bool conjugate2=false) const {return contractions(findindexlist(names),rhs2,rhs2.findindexlist(names),alpha,conjugate1,conjugate2); };
 
-	 void addgroupcontraction(const Tensor &rhs1, int group, const Tensor &rhs2, int rhsgroup, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false); //over all indices in a group of same symmetry; rhs1 will have more significant non-contracted indices in the result than rhs2
-	 inline Tensor groupcontraction(int group, const Tensor &rhs, int rhsgroup, T alpha=1, bool conjugate1=false, bool conjugate=false) const {Tensor<T> r; r.addgroupcontraction(*this,group,rhs,rhsgroup,alpha,0,true, conjugate1, conjugate); return r; };
+	void addgroupcontraction(const Tensor &rhs1, int group, const Tensor &rhs2, int rhsgroup, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false); //over all indices in a group of same symmetry; rhs1 will have more significant non-contracted indices in the result than rhs2
+	inline Tensor groupcontraction(int group, const Tensor &rhs, int rhsgroup, T alpha=1, bool conjugate1=false, bool conjugate=false) const {Tensor<T> r; r.addgroupcontraction(*this,group,rhs,rhsgroup,alpha,0,true, conjugate1, conjugate); return r; };
+
+	void add_permuted_contractions(const char *antisymmetrizer, const Tensor &rhs1, const Tensor &rhs2, T alpha=1, T beta=1, bool conjugate1=false, bool conjugate2=false); //must be used on named tensors, finds repeated indices to contract and antisymmetrizes and permutes remaining indices to match *this
+
 
 	Tensor innercontraction(const INDEXLIST &il1, const INDEXLIST &il2) const; //contraction(s) inside this tensor
 	Tensor innercontraction(const INDEX &i1, const INDEX &i2) const {INDEXLIST il1(1); il1[0]=i1; INDEXLIST il2(1); il2[0]=i2; return innercontraction(il1,il2);};