apply_permutation_algebra for tensors from product rhs

t.cc

@@ -3451,7 +3451,7 @@ cout <<mm;
 cout <<m.getcount()<<" "<<v.getcount()<<" "<<mm.getcount()<<endl;
 }
 
-if(1)
+if(0)
 {
 NRVec<double> v(12);
 v.randomize(1.);
@@ -3462,7 +3462,7 @@ cout<<"abssorted\n";
 v.printsorted(cout,1,true);
 }
 
-if(1)
+if(0)
 {
 NRSMat<double> v(4);
 v.randomize(1.);
@@ -3471,7 +3471,7 @@ cout<<"smat sorted\n";
 v.printsorted(cout,1,false);
 }
 
-if(1)
+if(0)
 {
 NRMat<double> v(4,5);
 v.randomize(1.);
@@ -3481,7 +3481,57 @@ v.printsorted(cout,1,false);
 }
 
 
+if(1)
+{
+//grassmann product of n identical rank=2 tensors in m-dim space
+int n,m;
+cin >>n>>m;
 
+//generate the source tensor
+INDEXGROUP g;
+g.number=2;
+g.symmetry= 0;
+g.offset=0;
+g.range=m;
+Tensor<double> x(g);
+x.randomize(1);
 
+cout <<x;
 
+//generate antisymmetrizer of even indices, with identity on odd indices
+NRVec<NRVec_from1<int> > indexclasses(1);
+indexclasses[0].resize(n);
+for(int i=1; i<=n; ++i) {indexclasses[0][i]= i;}
+PermutationAlgebra<int,int> a=general_antisymmetrizer(indexclasses,0,true);
+//antisymmetrize only in the even indices
+PermutationAlgebra<int,double> b(a.size());
+for(int i=0; i<a.size(); ++i) 
+	{
+	b[i].weight=a[i].weight;
+	b[i].perm.resize(2*n);
+	for(int j=1; j<=n; ++j)
+		{
+		b[i].perm[2*j-1] = 2*j-1;
+		b[i].perm[2*j] = 2*a[i].perm[j];
 		}
+	}
+cout <<b;
+
+//prepare output tensor (ignoring its symmetry)
+INDEXGROUP gg;
+gg.number= 2*n;
+gg.symmetry= 0;
+gg.offset=0;
+gg.range=m;
+Tensor<double> y(gg);
+
+NRVec<Tensor<double> > rhsvec(n);
+for(int i=0; i<n; ++i) rhsvec[i]=x;
+
+y.apply_permutation_algebra(rhsvec,b,false,1.,0.);
+
+cout <<y;
+}
+
+
+}//main

tensor.cc

@@ -924,6 +924,32 @@ for(int p=0; p<help_pa<T>->size(); ++p)
 }
 
 
+static int help_tn;
+template<typename T>
+static Tensor<T> *help_tv;
+
+template<typename T>
+static void permutationalgebra_callback2(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);
+	T product;
+	//build product from the vector of tensors using the Jp index
+	int rankshift=0;
+	for(int i=0; i<help_tn; ++i)
+		{
+		int rank = help_tv<T>[i].rank();
+		FLATINDEX indi = Jp.subvector(rankshift,rankshift+rank-1);
+		if(i==0) product= help_tv<T>[i](indi);
+		else product *= help_tv<T>[i](indi);
+		rankshift +=  rank;
+		}
+        *v += help_alpha<T> * (*help_pa<T>)[p].weight * product;
+        }
+}
+
 
 template<typename T>
 void Tensor<T>::apply_permutation_algebra(const Tensor<T> &rhs, const PermutationAlgebra<int,T> &pa, bool inverse, T alpha, T beta)
@@ -932,6 +958,8 @@ if(beta!=(T)0) {if(beta!=(T)1) *this *= beta;} else clear();
 if(alpha==(T)0) return;
 copyonwrite();
 
+if(rank()!=rhs.rank()) laerror("rank mismatch in apply_permutation_algebra");
+
 help_t<T> = const_cast<Tensor<T> *>(&rhs);
 help_pa<T> = &pa;
 help_inverse = inverse;
@@ -940,6 +968,29 @@ loopover(permutationalgebra_callback);
 }
 
 
+template<typename T>
+void Tensor<T>::apply_permutation_algebra(const NRVec<Tensor<T> > &rhsvec, const PermutationAlgebra<int,T> &pa, bool inverse, T alpha, T beta)
+{
+if(beta!=(T)0) {if(beta!=(T)1) *this *= beta;} else clear();
+if(alpha==(T)0) return;
+copyonwrite();
+
+int totrank=0;
+for(int i=0; i<rhsvec.size(); ++i) totrank+=rhsvec[i].rank();
+if(totrank!=rank()) laerror("rank mismatch in apply_permutation_algebra");
+
+help_tv<T> = const_cast<Tensor<T> *>(&rhsvec[0]);
+help_tn = rhsvec.size();
+help_pa<T> = &pa;
+help_inverse = inverse;
+help_alpha<T> = alpha;
+loopover(permutationalgebra_callback2);
+}
+
+
+
+
+
 template<typename T>
 void Tensor<T>::split_index_group(int group)
 {

tensor.h

@@ -41,8 +41,7 @@
 //@@@ will not be particularly efficient
 //
 //@@@permutation of individual indices - chech the indices in sym groups remain adjacent, calculate result's shape, loopover the result and permute using unwind_callback
-//@@@todo - implement index names - flat vector of names, and contraction by named index list
+//@@@todo!!! - implement index names - flat vector of names, and contraction by named index list
-//@@@todo grassmann product and support for RDM and cumulat stuff
 
 namespace LA {
 
@@ -217,8 +216,12 @@ 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; };
 
 	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  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
+//								 PermutationAlgebra can represent e.g. general_antisymmetrizer in Kucharski-Bartlett notation, or Grassmann products building RDM from cumulants
+//								 Note that *this tensor can be e.g. antisymmetric while rhs is not and is being antisymmetrized by the PermutationAlgebra
+//								 The efficiency is not optimal, even when avoiding the outer product, the calculation is done indexing element by element
+//								 More efficient would be applying permutation algebra symbolically and efficiently computing term by term
 	void split_index_group(int group); //formal split of a non-symmetric index group WITHOUT the need for data reorganization
 	void merge_adjacent_index_groups(int groupfrom, int groupto); //formal merge of non-symmetric index groups  WITHOUT the need for data reorganization
 	Tensor merge_index_groups(const NRVec<int> &groups) const;