From b3c7d212684befbf5ee062fb00199630800f27f9 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Thu, 6 Nov 2025 15:07:13 +0100 Subject: [PATCH] tensor: addgroupcontractions --- t.cc | 4 +- tensor.cc | 152 ++++++++++++++++++++++++++++++++++++++++++++++++------ tensor.h | 8 +++ 3 files changed, 146 insertions(+), 18 deletions(-) diff --git a/t.cc b/t.cc index 9d0a973..dd77fb3 100644 --- a/t.cc +++ b/t.cc @@ -3362,7 +3362,7 @@ for(int i=0; i +Tensor Tensor::unwind_index_group(int group) const +{ +if(group==0) return *this; //is already the least significant group +if(group<0||group>=shape.size()) laerror("wrong group number in unwind_index_group"); + +NRPerm p(shape.size()); +p[1]= 1+group; +int ii=1; +if(ii==1+group) ii++; //skip this +for(int i=2; i<=shape.size(); ++i) + { + p[i]=ii++; + if(ii==1+group) ii++; //skip this + } +if(!p.is_valid()) laerror("internal error in unwind_index"); +return permute_index_groups(p); +} + + + template Tensor Tensor::unwind_index(int group, int index) const { if(group<0||group>=shape.size()) laerror("wrong group number in unwind_index"); if(index<0||index>=shape[group].number) laerror("wrong index number in unwind_index"); -if(shape[group].number==1) //single index in the group + +if(group==0 && index==0 && shape[0].number == 1) return *this; +if(group==0 && index==0 && shape[0].symmetry==0) //formal split of 1 index without data rearrangement { - if(group==0) return *this; //is already the least significant group - NRPerm p(shape.size()); - p[1]= 1+group; - int ii=1; - if(ii==1+group) ii++; //skip this - for(int i=2; i<=shape.size(); ++i) - { - p[i]=ii++; - if(ii==1+group) ii++; //skip this - } - if(!p.is_valid()) laerror("internal error in unwind_index"); - return permute_index_groups(p); + Tensor r(*this); + r.split_index_group1(0); + return r; } +if(shape[group].number==1) return unwind_index_group(group); //single index in the group //general case - recalculate the shape and allocate the new tensor NRVec newshape(shape.size()+1); @@ -1088,6 +1103,65 @@ cblas_zgemm(CblasRowMajor, CblasNoTrans, (conjugate?CblasConjTrans:CblasTrans), +template +void Tensor::addgroupcontraction(const Tensor &rhs1, int group, const Tensor &rhs, int rhsgroup, T alpha, T beta, bool doresize, bool conjugate1, bool conjugate) +{ +if(group<0||group>=rhs1.shape.size()) laerror("wrong group number in contraction"); +if(rhsgroup<0||rhsgroup>=rhs.shape.size()) laerror("wrong rhsgroup number in contraction"); +if(rhs1.shape[group].offset != rhs.shape[rhsgroup].offset) laerror("incompatible index offset in contraction"); +if(rhs1.shape[group].range != rhs.shape[rhsgroup].range) laerror("incompatible index range in contraction"); +if(rhs1.shape[group].symmetry != rhs.shape[rhsgroup].symmetry) laerror("incompatible index symmetry in addgroupcontraction"); +if(rhs1.shape[group].symmetry == 1) laerror("addgroupcontraction not implemented for symmetric index groups"); +#ifdef LA_TENSOR_INDEXPOSITION +if(rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex == false) laerror("can contact only upper with lower index"); +#endif + +const Tensor u = conjugate1? (rhs1.unwind_index_group(group)).conjugate() : rhs1.unwind_index_group(group); +const Tensor rhsu = rhs.unwind_index_group(rhsgroup); + +NRVec newshape(u.shape.size()+rhsu.shape.size()-2); +int ii=0; +for(int i=1; i newnames; +if(u.is_named() && rhsu.is_named()) + { + for(int i=0; i(nn,mm,kk,&data[0],&u.data[0], &rhsu.data[0],factor,beta,conjugate); +} + + //Conntraction could be implemented without the temporary storage for unwinding, but then we would need @@ -1168,9 +1242,8 @@ for(int i=0; i u = rhs1.unwind_indices(il1); -if(conjugate1) u.conjugateme(); -Tensor rhsu = rhs2.unwind_indices(il2); +const Tensor u = conjugate1? (rhs1.unwind_indices(il1)).conjugateme() : rhs1.unwind_indices(il1); +const Tensor rhsu = rhs2.unwind_indices(il2); NRVec newshape(u.shape.size()+rhsu.shape.size()-2*il1.size()); @@ -1337,6 +1410,36 @@ if(data.size()!=newsize) laerror("internal error in split_index_group"); } +template +void Tensor::split_index_group1(int group) +{ +if(group<0||group >= shape.size()) laerror("illegal index group number"); +if(shape[group].number==1) return; //nothing to split +if(shape[group].symmetry!=0) laerror("only non-symmetric index group can be splitted, use flatten instead"); + +NRVec newshape(shape.size()+1); +int gg=0; +for(int g=0; g void Tensor:: merge_adjacent_index_groups(int groupfrom, int groupto) { @@ -1391,6 +1494,23 @@ return r; } + +template +T Tensor::dot(const Tensor &rhs) const +{ +if(shape!=rhs.shape) laerror("incompatible tensor shapes in dot"); +if(is_named() && rhs.is_named() && names!=rhs.names) laerror("incompatible tensor index names in dot"); +T factor=1; +for(int i=0; i::normtype norm() const {return data.norm();}; inline Tensor operator*(const Tensor &rhs) const {return Tensor(rhs.shape.concat(shape),data.otimes2vec(rhs.data),rhs.names.concat(names));} //outer product, rhs indices will be the less significant + T dot(const Tensor &rhs) const; //scalar product (full contraction), in complex case is automatically conjugated, not for symmetric group indices Tensor& conjugateme() {data.conjugateme(); return *this;}; inline Tensor conjugate() const {Tensor r(*this); r.conjugateme(); return r;}; @@ -275,6 +278,7 @@ public: Tensor permute_index_groups(const NRPerm &p) const; //rearrange the tensor storage permuting index groups as a whole + Tensor unwind_index_group(int group) const; //make the index group leftmost (least significant) Tensor unwind_index(int group, int index) const; //separate an index from a group and expand it to full range as the least significant one (the leftmost one) Tensor unwind_index(const INDEX &I) const {return unwind_index(I.group,I.index);}; Tensor unwind_index(const INDEXNAME &N) const {return unwind_index(findindex(N));}; @@ -295,6 +299,9 @@ 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; }; inline Tensor contractions(const Tensor &rhs2, const NRVec names, T alpha=1, bool conjugate1=false, bool conjugate2=false) const {return contractions(findindexlist(names),rhs2,rhs2.findindexlist(names),alpha,conjugate1,conjugate2); }; + void addgroupcontraction(const Tensor &rhs1, int group, const Tensor &rhs2, int rhsgroup, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false); //over all indices in a group of same symmetry; rhs1 will have more significant non-contracted indices in the result than rhs2 + inline Tensor groupcontraction(int group, const Tensor &rhs, int rhsgroup, T alpha=1, bool conjugate1=false, bool conjugate=false) const {Tensor r; r.addgroupcontraction(*this,group,rhs,rhsgroup,alpha,0,true, conjugate1, conjugate); 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)) @@ -304,6 +311,7 @@ public: // More efficient would be applying permutation algebra symbolically and efficiently computing term by term void split_index_group(int group); //formal in-place split of a non-symmetric index group WITHOUT the need for data reorganization or names rearrangement + void split_index_group1(int group); //formal in-place split of the leftmost index in a non-symmetric index group WITHOUT the need for data reorganization or names rearrangement void merge_adjacent_index_groups(int groupfrom, int groupto); //formal merge of non-symmetric index groups WITHOUT the need for data reorganization or names rearrangement Tensor merge_index_groups(const NRVec &groups) const;