From 518c75fb20d883f7dc9d433a3072bec7e380c869 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Fri, 3 May 2024 13:57:11 +0200
Subject: [PATCH] class tensor - apply permutation algebra

---
 t.cc      |  2 ++
 tensor.cc | 33 +++++++++++++++++++++++++++++++++
 tensor.h  |  8 ++++++--
 3 files changed, 41 insertions(+), 2 deletions(-)

diff --git a/t.cc b/t.cc
index 22bc995..1f9c871 100644
--- a/t.cc
+++ b/t.cc
@@ -3350,4 +3350,6 @@ for(int i=0; i<n; ++i)
 //cout <<c;
 }
 
+//test Tensor apply_permutation_algebra
+
 }
diff --git a/tensor.cc b/tensor.cc
index 52b1785..dfed10d 100644
--- a/tensor.cc
+++ b/tensor.cc
@@ -715,6 +715,39 @@ auxmatmult<T>(nn,mm,kk,&data[0],&u.data[0], &rhsu.data[0],alpha,beta);
 }
 
 
+template<typename T>
+static const PermutationAlgebra<int,T> *help_pa;
+
+static bool help_inverse;
+
+template<typename T>
+static T help_alpha;
+
+template<typename T>
+static void permutationalgebra_callback(const SUPERINDEX &I, T *v)
+{
+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_t<T>)(Jp);
+	}
+}
+
+
+
+template<typename T>
+void Tensor<T>::apply_permutation_algebra(const Tensor<T> &rhs, const PermutationAlgebra<int,T> &pa, bool inverse, T alpha, T beta)
+{
+if(beta!=(T)0) *this *= beta; else clear();
+if(alpha==(T)0) return;
+
+help_t<T> = const_cast<Tensor<T> *>(&rhs);
+help_pa<T> = &pa;
+help_inverse = inverse;
+help_alpha<T> = alpha;
+loopover(permutationalgebra_callback);
+}
 
 
 template class Tensor<double>;
diff --git a/tensor.h b/tensor.h
index 6f727c0..9926dae 100644
--- a/tensor.h
+++ b/tensor.h
@@ -183,8 +183,12 @@ public:
 	void addcontraction(const Tensor &rhs1, int group, int index, const Tensor &rhs, int rhsgroup, int rhsindex, T alpha=1, T beta=1, bool doresize=false);
 	inline Tensor contraction(int group, int index, const Tensor &rhs, int rhsgroup, int rhsindex, T alpha=1) const {Tensor<T> r; r.addcontraction(*this,group,index,rhs,rhsgroup,rhsindex,alpha,0,true); return r; }
 
-	//@@@ general antisymmetrization operator Kucharski style - or that will be left to a code generator?
-	//@@@symmetrize a group, antisymmetrize a group, expand a (anti)symmetric group - obecne symmetry change krome +1 na -1 vse mozne
+	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:
+//								 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
+
+
+	//TODO perhaps implement application of a permutation algebra to a product of several tensors
 };