From 08cdf7c97125062fc67d41fd826aee397ed53315 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Fri, 17 Oct 2025 15:07:43 +0200 Subject: [PATCH] apply_permutation_algebra for tensors from product rhs --- t.cc | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++--- tensor.cc | 51 ++++++++++++++++++++++++++++++++++++++++++++++++++ tensor.h | 9 ++++++--- 3 files changed, 110 insertions(+), 6 deletions(-) diff --git a/t.cc b/t.cc index 5d21575..cb31a96 100644 --- a/t.cc +++ b/t.cc @@ -3451,7 +3451,7 @@ cout < v(12); v.randomize(1.); @@ -3462,7 +3462,7 @@ cout<<"abssorted\n"; v.printsorted(cout,1,true); } -if(1) +if(0) { NRSMat v(4); v.randomize(1.); @@ -3471,7 +3471,7 @@ cout<<"smat sorted\n"; v.printsorted(cout,1,false); } -if(1) +if(0) { NRMat 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 x(g); +x.randomize(1); +cout < > indexclasses(1); +indexclasses[0].resize(n); +for(int i=1; i<=n; ++i) {indexclasses[0][i]= i;} +PermutationAlgebra a=general_antisymmetrizer(indexclasses,0,true); +//antisymmetrize only in the even indices +PermutationAlgebra b(a.size()); +for(int i=0; i y(gg); + +NRVec > rhsvec(n); +for(int i=0; i->size(); ++p) } +static int help_tn; +template +static Tensor *help_tv; + +template +static void permutationalgebra_callback2(const SUPERINDEX &I, T *v) +{ +FLATINDEX J = superindex2flat(I); +for(int p=0; p->size(); ++p) + { + FLATINDEX Jp = J.permuted((*help_pa)[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[i].rank(); + FLATINDEX indi = Jp.subvector(rankshift,rankshift+rank-1); + if(i==0) product= help_tv[i](indi); + else product *= help_tv[i](indi); + rankshift += rank; + } + *v += help_alpha * (*help_pa)[p].weight * product; + } +} + template void Tensor::apply_permutation_algebra(const Tensor &rhs, const PermutationAlgebra &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 = const_cast *>(&rhs); help_pa = &pa; help_inverse = inverse; @@ -940,6 +968,29 @@ loopover(permutationalgebra_callback); } +template +void Tensor::apply_permutation_algebra(const NRVec > &rhsvec, const PermutationAlgebra &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 = const_cast *>(&rhsvec[0]); +help_tn = rhsvec.size(); +help_pa = &pa; +help_inverse = inverse; +help_alpha = alpha; +loopover(permutationalgebra_callback2); +} + + + + + template void Tensor::split_index_group(int group) { diff --git a/tensor.h b/tensor.h index 8d2afef..a3aa56b 100644 --- a/tensor.h +++ b/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 grassmann product and support for RDM and cumulat stuff +//@@@todo!!! - implement index names - flat vector of names, and contraction by named index list 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 r; r.addcontractions(*this,il1,rhs2,il2,alpha,0,true,conjugate1, conjugate2); return r; }; void apply_permutation_algebra(const Tensor &rhs, const PermutationAlgebra &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 &rhsvec, const PermutationAlgebra &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 &groups) const;