From 417a7d1d1a6a933dd2181a6860b937969785cb61 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Tue, 18 Nov 2025 15:11:18 +0100
Subject: [PATCH] tensor: add_permuted_contractions works

---
 t.cc      | 46 ++++++++++++++++++++++++++++++++++++++++++++++
 tensor.cc | 21 +++++++--------------
 2 files changed, 53 insertions(+), 14 deletions(-)

diff --git a/t.cc b/t.cc
index f789a56..72c72a2 100644
--- a/t.cc
+++ b/t.cc
@@ -4169,6 +4169,52 @@ cout<<"resulting index order = "<<y.names<<endl;
 }
 
 
+if(0)
+{
+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[] = {"m","n","j","i"};
+y.names=NRVec<INDEXNAME>(ylist);
+
+Tensor<double>z(s); z.clear();
+INDEXNAME zlist[] = {"m","n","k","l"};
+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(m,n,j,i);
+	zz.lhs(m,n,k,l) = s;
+	//cout <<"test "<<k<<" "<<l<<" "<<m<<" "<<n<<" "<<zz(k,l,m,n)<<" "<<zzz(m,n,k,l)<< " : "<<z(m,n,k,l) <<endl;
+	}
+
+cout <<"Error = "<<(z-zz).norm()<<endl;
+}
+
 if(0)
 {
 int nn;
diff --git a/tensor.cc b/tensor.cc
index f92e9bd..a1c1f32 100644
--- a/tensor.cc
+++ b/tensor.cc
@@ -2107,44 +2107,37 @@ if(antinames.size()<tmpnames.size())
 			antinames[lastname] = tmpnames[j];
 			}
 	}
-std::cout <<"final antisymmmetrizer names and groups"<<antinames<<antigroups;
+//std::cout <<"final antisymmmetrizer names and groups"<<antinames<<antigroups;
 //std::cout <<"LHS names = "<<names<<std::endl; 
 //std::cout <<"TMP names = "<<tmpnames<<std::endl;
 
 //prepare the antisymmetrizer
 PermutationAlgebra<int,int> pa = general_antisymmetrizer(antigroups,-2,true);
-std::cout <<"initial antisymmetrizer = "<<pa;
+//std::cout <<"initial antisymmetrizer = "<<pa;
 
 //find permutation between antisym and TMP index order 
 NRPerm<int> antiperm=find_indexperm(antinames,tmpnames);
-std::cout<<"permutation from rhs to antisymmetrizer = "<<antiperm<<std::endl;
+//std::cout<<"permutation from rhs to antisymmetrizer = "<<antiperm<<std::endl;
 
 //conjugate the PA by antiperm or its inverse
-pa= pa.conjugated_by(antiperm,true);//@@@ or false?
+pa= pa.conjugated_by(antiperm,true);
 
 //find permutation which will bring indices of tmp to the order as in *this: this->names[i] == tmpnames[p[i]]
 NRPerm<int> basicperm=find_indexperm(names,tmpnames);
-std::cout <<"permutation from rhs to lhs = "<<basicperm<<std::endl;
-
-//@@@PermutationAlgebra<int,T> pb = basicperm*pa; //for apply_permutation_algebra(tmp,pb,false,(T)1,beta);
 PermutationAlgebra<int,T> pb = basicperm.inverse()*pa;
-//@@@ OR PermutationAlgebra<int,T> pb = basicperm.inverse()*pa; apply_permutation_algebra(tmp,pb,true,(T)1,beta);
+//std::cout <<"permutation from rhs to lhs = "<<basicperm<<std::endl;
 
-//save some work if the PA is trivial
+//save some work if the PA is trivial - contract only
 if(pb.is_identity())
 	{
-std::cout <<"simplified version\n";
 	addcontractions(rhs1,il1,rhs2,il2,alpha,beta,false,conjugate1,conjugate2);
 	return;
 	}
 
-std::cout <<"full version\n";
-//create an intermediate contracted tensir and the permute it
+//create an intermediate contracted tensor and the permute it
 Tensor<T> tmp=rhs1.contractions(il1,rhs2,il2,alpha,conjugate1,conjugate2);
 if(rank()!=tmp.rank()) laerror("rank mismatch in add_permuted_contractions");
 if(tmp.names!=tmpnames) laerror("internal error in add_permuted_contractions");
-
-//@@@apply_permutation_algebra(tmp,pb,false,(T)1,beta);
 apply_permutation_algebra(tmp,pb,true,(T)1,beta);
 }