jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

working on add_permuted_contractions

Browse Source
This commit is contained in:
Jiri Pittner
2025-11-14 16:16:22 +01:00
parent 89cc0c5b1e
commit 3f7586378d
3 changed files with 127 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
t.cc
View File

@@ -3568,7 +3568,7 @@ v.printsorted(cout,1,false);
}
if(1)
if(0)
{
//grassmann product of n identical rank=2 tensors in m-dim space
int n,m;
@@ -4164,5 +4164,52 @@ cout<<"resulting index order = "<<y.names<<endl;
}
if(1)
{
int nn;
cin>>nn;
int r=4;
NRVec<INDEXGROUP> s(r);
for(int i=0; i<r; ++i)
{
s[i].number=1;
s[i].symmetry=0;
s[i].range=nn;
s[i].offset=0;
}
Tensor<double> x(s); x.randomize(1.);
INDEXNAME xlist[] = {"i","j","k","l"};
x.names=NRVec<INDEXNAME>(xlist);
Tensor<double> y(s); y.randomize(1.);
INDEXNAME ylist[] = {"n","m","j","i"};
y.names=NRVec<INDEXNAME>(ylist);
Tensor<double>z(s); z.clear();
INDEXNAME zlist[] = {"k","l","m","n"};
z.names=NRVec<INDEXNAME>(zlist);
Tensor<double>zz(z);
INDEX i1[]={{0,0},{1,0}};
INDEX i2[]={{3,0},{2,0}};
NRVec<INDEX> il1(i1);
NRVec<INDEX> il2(i2);
Tensor<double> zzz = x.contractions(il1,y,il2,1,false,false);
//cout <<"zzz names = "<<zzz.names<<endl;
z.add_permuted_contractions("",x,y,1,0,false,false);
zz.copyonwrite();
for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int n=0; n<nn; ++n)
{
double s=0;
for(int i=0; i<nn; ++i) for(int j=0; j<nn; ++j) s += x(i,j,k,l)*y(n,m,j,i);
zz.lhs(k,l,m,n) = s;
//cout <<"test "<<k<<" "<<l<<" "<<m<<" "<<n<<" "<<zz(k,l,m,n)<<" "<<zzz(n,m,k,l)<< " : "<<z(k,l,m,n) <<endl;
}
cout <<"Error = "<<(z-zz).norm()<<endl;
}
}//main

76
tensor.cc
View File

@@ -1198,7 +1198,7 @@ if(rhsindex<0||rhsindex>=rhs.shape[rhsgroup].number) laerror("wrong index numbe
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");
#ifdef LA_TENSOR_INDEXPOSITION
if(rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex == false) laerror("can contact only upper with lower index");
if(rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex == false) laerror("can contract only upper with lower index");
#endif
const Tensor<T> u = conjugate1? (rhs1.unwind_index(group,index)).conjugate() : rhs1.unwind_index(group,index);
@@ -1458,9 +1458,11 @@ if(rank()!=rhs.rank()) laerror("rank mismatch in apply_permutation_algebra");
if(rhs.is_named())
{
NRVec<INDEXNAME> namperm = rhs.names.permuted(pa[0].perm,inverse);
NRVec<INDEXNAME> namperm = rhs.names.permuted(pa[0].perm,!inverse);
if(is_named())
{
std::cout <<"LHS names = "<<names<<std::endl;
std::cout <<"permuted RHS names = "<<namperm<<std::endl;
if(names!=namperm) laerror("inconsistent index names in apply_permutation_algebra");
}
else
@@ -1492,7 +1494,7 @@ if(allnamed)
{
NRVec<INDEXNAME> allrhsnames=rhsvec[0].names;
for(int i=1; i<rhsvec.size(); ++i) allrhsnames.concatme(rhsvec[i].names);
NRVec<INDEXNAME> namperm = allrhsnames.permuted(pa[0].perm,inverse);
NRVec<INDEXNAME> namperm = allrhsnames.permuted(pa[0].perm,!inverse);
if(is_named())
{
if(names!=namperm) laerror("inconsistent index names in apply_permutation_algebra");
@@ -1907,6 +1909,74 @@ return s;
}
template<typename T>
void Tensor<T>::add_permuted_contractions(const char *antisymmetrizer, const Tensor &rhs1, const Tensor &rhs2, T alpha, T beta, bool conjugate1, bool conjugate2)
{
if(!rhs1.is_uniquely_named()||!rhs2.is_uniquely_named()|| !is_uniquely_named()) laerror("tensors must have unique named indices in add_permuted_contractions");
//find contraction indices
int nc=0;
for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs1.names.size(); ++j)
if(rhs1.names[i]==rhs2.names[j])
{
//std::cout << "found contraction "<<nc<<" th. index = "<<rhs1.names[i]<<std::endl;
#ifdef LA_TENSOR_INDEXPOSITION
int group1=indexposition(i,rhs1.shape);
int group2=indexposition(j,rhs2.shape);
if(rhs1.shape[group1].upperindex ^ rhs2.shape[group2].upperindex == false) laerror("can contract only upper with lower index");
#endif
++nc;
}
INDEXLIST il1(nc),il2(nc);
int ii=0;
for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs1.names.size(); ++j)
if(rhs1.names[i]==rhs2.names[j])
{
il1[ii] = indexposition(i,rhs1.shape);
il2[ii] = indexposition(j,rhs2.shape);
++ii;
}
//std::cout<<"contraction list1 = "<<il1<<std::endl;
//std::cout<<"contraction list2 = "<<il2<<std::endl;
Tensor<T> tmp=rhs1.contractions(il1,rhs2,il2,alpha,conjugate1,conjugate2);
if(rank()!=tmp.rank()) laerror("rank mismatch in add_permuted_contractions");
//generate the antisymmetrizer, adding also indices not involved as a constant subpermutation
//@@@
PermutationAlgebra<int,T> pa(1);
pa[0].weight=1;
pa[0].perm.resize(rank());
pa[0].perm.identity();
//std::cout <<"LHS names = "<<names<<std::endl;
//std::cout <<"TMP names = "<<tmp.names<<std::endl;
//find permutation which will bring indices of tmp to the order as in *this: this->names[i] == tmp.names[p[i]]
NRPerm<int> basicperm(rank());
basicperm.clear();
for(int i=0; i<rank(); ++i)
{
for(int j=0; j<tmp.rank(); ++j)
if((* const_cast<const NRVec<INDEXNAME> *>(&names))[i] == tmp.names[j])
{
basicperm[i+1]=j+1;
break;
}
}
if(!basicperm.is_valid()) laerror("indices mismatch between lhs and rhs in add_permuted_contractions");
std::cout <<"Basic permutation = "<<basicperm<<std::endl;
PermutationAlgebra<int,T> pb = basicperm*pa;
apply_permutation_algebra(tmp,pb,false,(T)1/(T)pb.size(),beta);
//equivalently possible
//PermutationAlgebra<int,T> pb = basicperm.inverse()*pa;
//apply_permutation_algebra(tmp,pb,true,(T)1/(T)pb.size(),beta);
}
template class Tensor<double>;

5
tensor.h
View File

@@ -230,7 +230,7 @@ public:
bool is_compressed() const {for(int i=0; i<shape.size(); ++i) if(shape[i].number>1&&shape[i].symmetry!=0) return true; return false;};
bool has_symmetry() const {for(int i=0; i<shape.size(); ++i) if(shape[i].symmetry!=0) return true; return false;};
void clear() {data.clear();};
void defaultnames() {names.resize(rank()); for(int i=0; i<rank(); ++i) sprintf(names[i].name,"i%03d",i);}
void defaultnames(const char *basename="i") {names.resize(rank()); for(int i=0; i<rank(); ++i) sprintf(names[i].name,"%s%03d",basename,i);}
int rank() const {return myrank;};
int calcrank(); //is computed from shape
void canonicalize_shape();
@@ -325,6 +325,9 @@ public:
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<T> r; r.addgroupcontraction(*this,group,rhs,rhsgroup,alpha,0,true, conjugate1, conjugate); return r; };
void add_permuted_contractions(const char *antisymmetrizer, const Tensor &rhs1, const Tensor &rhs2, T alpha=1, T beta=1, bool conjugate1=false, bool conjugate2=false); //must be used on named tensors, finds repeated indices to contract and antisymmetrizes and permutes remaining indices to match *this
Tensor innercontraction(const INDEXLIST &il1, const INDEXLIST &il2) const; //contraction(s) inside this tensor
Tensor innercontraction(const INDEX &i1, const INDEX &i2) const {INDEXLIST il1(1); il1[0]=i1; INDEXLIST il2(1); il2[0]=i2; return innercontraction(il1,il2);};
Tensor innercontraction(const NRVec<INDEXNAME> &nl1, const NRVec<INDEXNAME> &nl2) const {return innercontraction(findindexlist(nl1),findindexlist(nl2));};