class tensor - apply permutation algebra

t.cc

@@ -3350,4 +3350,6 @@ for(int i=0; i<n; ++i)
 //cout <<c;
 }
 
+//test Tensor apply_permutation_algebra
+
 }

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

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?
+	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:
-	//@@@symmetrize a group, antisymmetrize a group, expand a (anti)symmetric group - obecne symmetry change krome +1 na -1 vse mozne
+//								 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
 };