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

Compare commits

base: jiri:5f74028ab65312a99d5fc355e6dff3825fd2fb17
Branches Tags
jiri:master
...
compare: jiri:71c890c39b5a6a693721adc4fc379de5b2e4a0e6
Branches Tags
jiri:master

2 Commits
5f74028ab6 ... 71c890c39b

Author SHA1 Message Date
Jiri Pittner
71c890c39b working on add_permuted_contractions 2025-11-16 14:56:51 +01:00
Jiri Pittner
ba5adcd5e6 rename conjugate_by to conjugated_by 2025-11-16 14:56:29 +01:00
5 changed files with 77 additions and 46 deletions
Expand all files Collapse all files

8
permutation.cc
View File

@@ -189,7 +189,7 @@ return r;
template <typename T>
NRPerm<T> NRPerm<T>::conjugate_by(const NRPerm<T> &q, bool inverse) const
NRPerm<T> NRPerm<T>::conjugated_by(const NRPerm<T> &q, bool inverse) const
{
#ifdef DEBUG
if(!this->is_valid() || !q.is_valid()) laerror("multiplication of invalid permutation");
@@ -224,7 +224,7 @@ return r;
template <typename T>
CyclePerm<T> CyclePerm<T>::conjugate_by(const CyclePerm<T> &q) const
CyclePerm<T> CyclePerm<T>::conjugated_by(const CyclePerm<T> &q) const
{
#ifdef DEBUG
if(!this->is_valid() || !q.is_valid()) laerror("multiplication of invalid permutation");
@@ -780,10 +780,10 @@ return res;
}
template <typename T, typename R>
PermutationAlgebra<T,R> PermutationAlgebra<T,R>::conjugate_by(const NRPerm<T> &q, bool reverse) const
PermutationAlgebra<T,R> PermutationAlgebra<T,R>::conjugated_by(const NRPerm<T> &q, bool reverse) const
{
PermutationAlgebra<T,R> res(this->size());
for(int i=0; i<this->size(); ++i) {res[i].perm = (*this)[i].perm.conjugate_by(q,reverse); res[i].weight=(*this)[i].weight;}
for(int i=0; i<this->size(); ++i) {res[i].perm = (*this)[i].perm.conjugated_by(q,reverse); res[i].weight=(*this)[i].weight;}
return res;
}

6
permutation.h
View File

@@ -73,7 +73,7 @@ public:
NRPerm operator*(const CyclePerm<T> &r) const;
template<typename R> PermutationAlgebra<T,R> operator*(const PermutationAlgebra<T,R> &pa) const;
NRPerm multiply(const NRPerm<T> &q, bool inverse) const; //multiplication but optionally q inversed
NRPerm conjugate_by(const NRPerm &q, bool reverse=false) const; //q^-1 p q or q p q^-1
NRPerm conjugated_by(const NRPerm &q, bool reverse=false) const; //q^-1 p q or q p q^-1
NRPerm commutator(const NRPerm &q, bool inverse=false) const; //p^-1 q^-1 p q or q^-1 p^-1 q p
int parity() const; //returns +/- 1
void randomize(void); //uniformly random by Fisher-Yates shuffle
@@ -199,7 +199,7 @@ public:
PermutationAlgebra operator*(const NRPerm<T> &rhs) const; //applied to all terms
PermutationAlgebra cutleft(int n) const; //applied to all terms
PermutationAlgebra cutright(int n) const; //applied to all terms
PermutationAlgebra conjugate_by(const NRPerm<T> &q, bool reverse=false) const; //q^-1 p q or q p q^-1 , applied to all terms
PermutationAlgebra conjugated_by(const NRPerm<T> &q, bool reverse=false) const; //q^-1 p q or q p q^-1 , applied to all terms
PermutationAlgebra commutator(const NRPerm<T> &q, bool inverse=false) const; //applied to all terms
PermutationAlgebra operator&(const PermutationAlgebra &rhs) const; //each term with each
PermutationAlgebra operator|(const PermutationAlgebra &rhs) const; //each term with each
@@ -237,7 +237,7 @@ public:
void readfrom(const std::string &line);
CyclePerm operator*(const CyclePerm &q) const; //q is rhs and applied first, this applied second
NRPerm<T> operator*(const NRPerm<T> &r) const;
CyclePerm conjugate_by(const CyclePerm &q) const; //q^-1 p q
CyclePerm conjugated_by(const CyclePerm &q) const; //q^-1 p q
PERM_RANK_TYPE order() const; //lcm of cycle lengths
bool operator==(const CyclePerm &rhs) const {return NRPerm<T>(*this) == NRPerm<T>(rhs);}; //cycle representation is not unique, cannot inherit operator== from NRVec
void simplify(bool keep1=false); //remove cycles of size 0 or 1

2
t.cc
View File

@@ -4254,7 +4254,7 @@ for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int
Tensor<double> zzz(s);
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)
{
zzz.lhs(k,l,m,n) = 1/3.*(zz(k,l,m,n)-zz(l,k,m,n)+zz(m,k,l,n));
zzz.lhs(k,l,m,n) = (zz(k,l,m,n)-zz(l,k,m,n)+zz(m,k,l,n));
}
cout <<"Error = "<<(z-zzz).norm()<<endl;

87
tensor.cc
View File

@@ -85,7 +85,11 @@ LA_largeindex subindex(int *sign, const INDEXGROUP &g, const NRVec<LA_index> &I)
#ifdef DEBUG
if(I.size()<=0) laerror("empty index group in subindex");
if(g.number!=I.size()) laerror("mismatch in the number of indices in a group");
for(int i=0; i<I.size(); ++i) if(I[i]<g.offset || I[i] >= g.offset+g.range) laerror("index out of range in tensor subindex");
for(int i=0; i<I.size(); ++i) if(I[i]<g.offset || I[i] >= g.offset+g.range)
{
std::cout<<"TENSOR INDEX PROBLEM in group " <<g<<" with index "<<I<<std::endl;
laerror("index out of range in tensor subindex");
}
#endif
switch(I.size()) //a few special cases for efficiency
@@ -232,7 +236,7 @@ for(int i=0; i<I.size(); ++i)
{
if(I[i][j] <shape[i].offset || I[i][j] >= shape[i].offset+shape[i].range)
{
std::cerr<<"error in index group no. "<<i<<" index no. "<<j<<std::endl;
std::cout<<"TENSOR INDEX PROBLEM group no. "<<i<<" index no. "<<j<<" should be between "<<shape[i].offset<<" and "<<shape[i].offset+shape[i].range-1<<std::endl;
laerror("tensor index out of range");
}
}
@@ -1883,6 +1887,8 @@ if(is_named()) r.names=names.permuted(basicperm,true);
return r;
}
template<typename T>
void Tensor<T>::canonicalize_shape()
{
@@ -1999,7 +2005,7 @@ foundname:
finished:
names0= NRVec<INDEXNAME> (names);
std::cout <<"names parsed "<<names0;
//std::cout <<"names parsed "<<names0;
groups0.resize(groups.size());
int i=0;
@@ -2007,7 +2013,7 @@ for(typename std::list<std::list<int> >::const_iterator ii=groups.begin(); ii!=g
{
groups0[i++] = NRVec_from1<int>(*ii);
}
std::cout<<"groups parsed "<<groups0;
//std::cout<<"groups parsed "<<groups0;
free(txt);
}
@@ -2033,10 +2039,14 @@ for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs1.names.size(); ++j)
}
INDEXLIST il1(nc),il2(nc);
bitvector is_c_index1(rhs1.names.size());
bitvector is_c_index2(rhs2.names.size());
int ii=0;
for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs1.names.size(); ++j)
for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs2.names.size(); ++j)
if(rhs1.names[i]==rhs2.names[j])
{
is_c_index1.set(i);
is_c_index2.set(j);
il1[ii] = indexposition(i,rhs1.shape);
il2[ii] = indexposition(j,rhs2.shape);
++ii;
@@ -2045,8 +2055,11 @@ for(int i=0; i<rhs1.names.size(); ++i) for(int j=0; j<rhs1.names.size(); ++j)
//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");
//make a dry-run to get only the list of names, concatenate names of both tensors unless they are the contraction ones
std::list<INDEXNAME> listnames;
for(int j=0; j<rhs2.names.size(); ++j) if(!is_c_index2[j]) listnames.push_back(rhs2.names[j]);
for(int j=0; j<rhs1.names.size(); ++j) if(!is_c_index1[j]) listnames.push_back(rhs1.names[j]);
NRVec<INDEXNAME> tmpnames(listnames);
//generate the antisymmetrizer, adding also indices not involved as a constant subpermutation
//
@@ -2056,14 +2069,14 @@ NRVec<NRVec_from1<int> > antigroups;
parse_antisymmetrizer(antisymmetrizer,antigroups,antinames);
//check the names make sense and fill in the possibly missing ones as separate group
if(antinames.size()>tmp.names.size()) laerror("too many indices in the antisymmetrizet");
bitvector isexplicit(tmp.names.size());
if(antinames.size()>tmpnames.size()) laerror("too many indices in the antisymmetrizet");
bitvector isexplicit(tmpnames.size());
isexplicit.clear();
for(int i=0; i<antinames.size(); ++i)
{
for(int j=0; j<i; ++j) if(antinames[i]==antinames[j]) laerror("repeated names in the antisymmetrizer");
for(int j=0; j<tmp.names.size(); ++j)
if(antinames[i]==tmp.names[j])
for(int j=0; j<tmpnames.size(); ++j)
if(antinames[i]==tmpnames[j])
{
isexplicit.set(j);
goto namefound;
@@ -2074,15 +2087,15 @@ for(int i=0; i<antinames.size(); ++i)
if(isexplicit.population()!=antinames.size()) laerror("internal error in add_permuted_contractions");
//fill in additional names
if(antinames.size()<tmp.names.size())
if(antinames.size()<tmpnames.size())
{
int lastgroup=antigroups.size()-1;
int lastclass;
if(antigroups.size()==0) lastclass=0;
else lastclass=antigroups[antigroups.size()-1][antigroups[antigroups.size()-1].size()];
int lastname=antinames.size()-1;
antinames.resize(tmp.names.size(),true);
antigroups.resize(antigroups.size()+tmp.names.size()-antinames.size(),true);
antigroups.resize(antigroups.size()+tmpnames.size()-antinames.size(),true);
antinames.resize(tmpnames.size(),true);
for(int j=0; j<names.size(); ++j)
if(!isexplicit[j])
{
@@ -2090,37 +2103,49 @@ if(antinames.size()<tmp.names.size())
++lastname;
++lastgroup;
antigroups[lastgroup].resize(1);
antigroups[lastgroup][0]=lastclass;
antinames[lastname] = tmp.names[j];
antigroups[lastgroup][1]=lastclass;
antinames[lastname] = tmpnames[j];
}
}
std::cout <<"final antisymmmetrizer names and groups"<<antinames<<antigroups;
//std::cout <<"LHS names = "<<names<<std::endl;
//std::cout <<"TMP names = "<<tmpnames<<std::endl;
//prepare the antisymmetrizer
PermutationAlgebra<int,int> pa = general_antisymmetrizer(antigroups,-2,true);
std::cout <<"initial antisymmetrizer = "<<pa;
//find permutation between antisym and TMP index order
NRPerm<int> antiperm=find_indexperm(antinames,tmp.names);
NRPerm<int> antiperm=find_indexperm(antinames,tmpnames);
std::cout<<"permutation from rhs to antisymmetrizer = "<<antiperm<<std::endl;
//@@@conjugate the PA by antiperm or its inverse
//conjugate the PA by antiperm or its inverse
pa= pa.conjugated_by(antiperm,true);//@@@ or false?
//@@@recast the PA
//find permutation which will bring indices of tmp to the order as in *this: this->names[i] == tmpnames[p[i]]
NRPerm<int> basicperm=find_indexperm(names,tmpnames);
std::cout <<"permutation from rhs to lhs = "<<basicperm<<std::endl;
//@@@permutationalgebra::is_identity() a pokd ano, neaplikovat ale vratit tmp tenzor resp. s nim udelat axpy
//@@@PermutationAlgebra<int,T> pb = basicperm*pa; //for apply_permutation_algebra(tmp,pb,false,(T)1,beta);
PermutationAlgebra<int,T> pb = basicperm.inverse()*pa;
//@@@ OR PermutationAlgebra<int,T> pb = basicperm.inverse()*pa; apply_permutation_algebra(tmp,pb,true,(T)1,beta);
//std::cout <<"LHS names = "<<names<<std::endl;
//std::cout <<"TMP names = "<<tmp.names<<std::endl;
//save some work if the PA is trivial
if(pb.is_identity())
{
std::cout <<"simplified version\n";
addcontractions(rhs1,il1,rhs2,il2,alpha,beta,false,conjugate1,conjugate2);
return;
}
//find permutation which will bring indices of tmp to the order as in *this: this->names[i] == tmp.names[p[i]]
NRPerm<int> basicperm=find_indexperm(names,tmp.names);
std::cout <<"Basic permutation = "<<basicperm<<std::endl;
std::cout <<"full version\n";
//create an intermediate contracted tensir and the permute it
Tensor<T> tmp=rhs1.contractions(il1,rhs2,il2,alpha,conjugate1,conjugate2);
if(rank()!=tmp.rank()) laerror("rank mismatch in add_permuted_contractions");
if(tmp.names!=tmpnames) laerror("internal error in add_permuted_contractions");
//@@@probably only onw of those will work
PermutationAlgebra<int,T> pb = basicperm*pa;
apply_permutation_algebra(tmp,pb,false,(T)1/(T)pb.size(),beta);
//equivalently possible (for trivial pa)
//PermutationAlgebra<int,T> pb = basicperm.inverse()*pa;
//apply_permutation_algebra(tmp,pb,true,(T)1/(T)pb.size(),beta);
//@@@apply_permutation_algebra(tmp,pb,false,(T)1,beta);
apply_permutation_algebra(tmp,pb,true,(T)1,beta);
}

20
tensor.h
View File

@@ -267,23 +267,29 @@ public:
NRVec<INDEX> findindexlist(const NRVec<INDEXNAME> &names) const;
void renameindex(const INDEXNAME namfrom, const INDEXNAME nameto) {int i=findflatindex(namfrom); names[i]=nameto;};
inline Tensor& operator+=(const Tensor &rhs)
Tensor& operator+=(const Tensor &rhs)
{
#ifdef DEBUG
if(shape!=rhs.shape) laerror("incompatible tensors for operation");
#endif
if(is_named() && rhs.is_named() && names!=rhs.names) laerror("incompatible names for operation");
data+=rhs.data;
return *this;
}
inline Tensor& operator-=(const Tensor &rhs)
Tensor& operator-=(const Tensor &rhs)
{
#ifdef DEBUG
if(shape!=rhs.shape) laerror("incompatible tensors for operation");
#endif
if(is_named() && rhs.is_named() && names!=rhs.names) laerror("incompatible names for operation");
data-=rhs.data;
return *this;
}
Tensor& axpy(const T alpha, const Tensor &rhs)
{
if(shape!=rhs.shape) laerror("incompatible tensors for operation");
if(is_named() && rhs.is_named() && names!=rhs.names) laerror("incompatible names for operation");
data.axpy(alpha,rhs.data);
return *this;
}
inline Tensor operator+(const Tensor &rhs) const {Tensor r(*this); r+=rhs; return r;};
inline Tensor operator-(const Tensor &rhs) const {Tensor r(*this); r-=rhs; return r;};
@@ -613,7 +619,7 @@ NRMat<T> mat(t.data,range*(range-1)/2,range*(range-1)/2);
f=fourindex_dense<antisymtwoelectronrealdirac,T,I>(range,NRSMat<T>(mat)); //symmetrize mat
}
//@@@formal permutation of names inside a sym/antisy group (with possible sign change)
template <typename T>