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

bugfixes in tensor

Browse Source
This commit is contained in:
jiri committed 2026-10-28 00:04:20 +01:00
1 parent 79a0cf9247
commit cbf054fd7a
2 files changed
+49 -29
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

tensor.cc
+46 -26
View File
@@ -327,7 +327,7 @@ shape.put(fd,true);
groupsizes.put(fd,true); groupsizes.put(fd,true);
cumsizes.put(fd,true); cumsizes.put(fd,true);
data.put(fd,true); data.put(fd,true);
if(with_names && names.size()!=0) names.put(fd,true); if(with_names) names.put(fd,true); //also if empty (writes its length 0), get() always reads the names if with_names
} }
template<typename T> template<typename T>
@@ -335,7 +335,7 @@ void Tensor<T>::get(int fd, bool with_names)
{ {
shape.get(fd,true); shape.get(fd,true);
myrank=calcrank(); //is not stored but recomputed myrank=calcrank(); //is not stored but recomputed
groupsizes.put(fd,true); groupsizes.get(fd,true);
cumsizes.get(fd,true); cumsizes.get(fd,true);
data.get(fd,true); data.get(fd,true);
if(with_names) names.get(fd,true); if(with_names) names.get(fd,true);
@@ -403,7 +403,8 @@ Tensor<T>::Tensor(const NRSMat<T> &x)
myrank=2; myrank=2;
shape.resize(1); shape.resize(1);
shape[0].number=2; shape[0].number=2;
shape[0].symmetry=1; shape[0].symmetry= LA_traits<T>::is_complex() ? 2 : 1; //complex NRSMat is hermitian, same packed layout
//TODO: if NRSMat gets the planned hermitian flag, use 2 only for a hermitian complex x, 1 otherwise
#ifndef LA_TENSOR_ZERO_OFFSET #ifndef LA_TENSOR_ZERO_OFFSET
shape[0].offset=0; shape[0].offset=0;
#endif #endif
@@ -450,7 +451,8 @@ for(LA_index i = istart; i<=iend; ++i)
{ {
int sign; int sign;
//std::cout <<"TEST "<<t.index(&sign,I)<<" "; //std::cout <<"TEST "<<t.index(&sign,I)<<" ";
(*callback)(I,(*p)++); if(skipzeros) (*callback)(I,&t.data[t.index(&sign,I)]); //pointer from the index, the skipped blocks are not counted by p
else (*callback)(I,(*p)++);
} }
else else
{ {
@@ -478,6 +480,7 @@ for(int i=0; i<I.size(); ++i)
I[i].resize(sh->number); I[i].resize(sh->number);
I[i] = sh->offset; I[i] = sh->offset;
} }
data.copyonwrite(); //the callback gets writable pointers
T *pp=&data[0]; T *pp=&data[0];
int ss=shape.size()-1; int ss=shape.size()-1;
const INDEXGROUP *sh = &(* const_cast<const NRVec<INDEXGROUP> *>(&shape))[ss]; const INDEXGROUP *sh = &(* const_cast<const NRVec<INDEXGROUP> *>(&shape))[ss];
@@ -519,7 +522,8 @@ for(LA_index i = istart; i<=iend; ++i)
{ {
int sign; int sign;
//std::cout <<"TEST "<<t.index(&sign,I)<<" "; //std::cout <<"TEST "<<t.index(&sign,I)<<" ";
(*callback)(I,(*p)++); if(skipzeros) (*callback)(I,&t.data[t.index(&sign,I)]); //pointer from the index, the skipped blocks are not counted by p
else (*callback)(I,(*p)++);
} }
else else
{ {
@@ -674,6 +678,7 @@ template<typename T>
void Tensor<T>::grouploopover(void (*callback)(const GROUPINDEX &, T *)) void Tensor<T>::grouploopover(void (*callback)(const GROUPINDEX &, T *))
{ {
GROUPINDEX I(shape.size()); GROUPINDEX I(shape.size());
data.copyonwrite(); //the callback gets writable pointers
T *pp= &data[0]; T *pp= &data[0];
loopovergroups(*this,shape.size()-1,&pp,I,callback); loopovergroups(*this,shape.size()-1,&pp,I,callback);
} }
@@ -965,6 +970,7 @@ if(group>=0) //single group
else else
{ {
Tensor<T> r(*this); Tensor<T> r(*this);
r.shape.copyonwrite(); //data stay shared, the change is formal
r.shape[group].symmetry=0; r.shape[group].symmetry=0;
return r; return r;
} }
@@ -979,7 +985,7 @@ if(group>=0) //single group
if(group<0 && !is_compressed()) if(group<0 && !is_compressed())
{ {
Tensor<T> r(*this); Tensor<T> r(*this);
for(int g=0; g<shape.size(); ++g) for(int g=shape.size()-1; g>=0; --g) //downward: splitting renumbers only the following groups
{ {
if(shape[g].number>1) r.split_index_group(g); if(shape[g].number>1) r.split_index_group(g);
} }
@@ -1165,8 +1171,8 @@ static void auxmatmult(int nn, int mm, int kk, T *r, const T *a, const T *b, T
{ {
for(int i=0; i<nn; ++i) for(int j=0; j<mm; ++j) for(int i=0; i<nn; ++i) for(int j=0; j<mm; ++j)
{ {
if(beta==0) r[i*mm+j]=0; else r[i*mm+j] *= beta; if(beta==(T)0) r[i*mm+j]=0; else r[i*mm+j] *= beta;
for(int k=0; k<kk; ++k) r[i*mm+j] += alpha * a[i*kk+k] * conjugate? LA_traits<T>::conjugate(b[j*kk+k]) : b[j*kk+k]; for(int k=0; k<kk; ++k) r[i*mm+j] += alpha * a[i*kk+k] * (conjugate? LA_traits<T>::conjugate(b[j*kk+k]) : b[j*kk+k]);
} }
} }
@@ -1203,7 +1209,7 @@ if(rhs1.shape[group].range != rhs.shape[rhsgroup].range) laerror("incompatible i
if(rhs1.shape[group].symmetry != rhs.shape[rhsgroup].symmetry) laerror("incompatible index symmetry in addgroupcontraction"); if(rhs1.shape[group].symmetry != rhs.shape[rhsgroup].symmetry) laerror("incompatible index symmetry in addgroupcontraction");
if(rhs1.shape[group].symmetry !=0 && rhs1.shape[group].symmetry != -1) laerror("addgroupcontraction only implemented for nonsymmetric and antisymmetric index groups"); if(rhs1.shape[group].symmetry !=0 && rhs1.shape[group].symmetry != -1) laerror("addgroupcontraction only implemented for nonsymmetric and antisymmetric index groups");
#ifdef LA_TENSOR_INDEXPOSITION #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 contact only upper with lower index");
#endif #endif
const Tensor<T> u = conjugate1? (rhs1.unwind_index_group(group)).conjugate() : rhs1.unwind_index_group(group); const Tensor<T> u = conjugate1? (rhs1.unwind_index_group(group)).conjugate() : rhs1.unwind_index_group(group);
@@ -1271,7 +1277,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].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].range != rhs.shape[rhsgroup].range) laerror("incompatible index range in contraction");
#ifdef LA_TENSOR_INDEXPOSITION #ifdef LA_TENSOR_INDEXPOSITION
if(rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex == false) laerror("can contract only upper with lower index"); if((rhs1.shape[group].upperindex ^ rhs.shape[rhsgroup].upperindex) == false) laerror("can contract only upper with lower index");
#endif #endif
const Tensor<T> u = conjugate1? (rhs1.unwind_index(group,index)).conjugate() : rhs1.unwind_index(group,index); const Tensor<T> u = conjugate1? (rhs1.unwind_index(group,index)).conjugate() : rhs1.unwind_index(group,index);
@@ -1331,8 +1337,19 @@ void Tensor<T>::addcontractions(const Tensor &rhs1, const INDEXLIST &il1, const
if(il1.size()!=il2.size()) laerror("mismatch in index lists in addcontractions"); if(il1.size()!=il2.size()) laerror("mismatch in index lists in addcontractions");
if(il1.size()==0) //add just outer product if(il1.size()==0) //add just outer product
{ {
if(beta!=(T)0) *this *= beta; else this->clear(); Tensor<T> p = (conjugate1?rhs1.conjugate():rhs1) * (conjugate2?rhs2.conjugate():rhs2);
*this += (conjugate1?rhs1.conjugate():rhs1) * (conjugate2?rhs2.conjugate():rhs2) * alpha; if(alpha!=(T)1) p *= alpha; //in place, no extra copy
if(doresize)
{
if(beta!= (T)0) laerror("resize in addcontractions requires beta=0");
*this = p;
}
else
{
if(shape!=p.shape) laerror("tensor shape mismatch in addcontraction");
if(is_named() && names!=p.names) laerror("remaining index names mismatch in addcontractions");
if(beta==(T)0) data=p.data; else {if(beta!=(T)1) *this *= beta; *this += p;}
}
return; return;
} }
@@ -1345,7 +1362,7 @@ for(int i=0; i<il1.size(); ++i)
if(rhs1.shape[il1[i].group].offset != rhs2.shape[il2[i].group].offset) laerror("incompatible index offset in contractions"); if(rhs1.shape[il1[i].group].offset != rhs2.shape[il2[i].group].offset) laerror("incompatible index offset in contractions");
if(rhs1.shape[il1[i].group].range != rhs2.shape[il2[i].group].range) laerror("incompatible index range in contractions"); if(rhs1.shape[il1[i].group].range != rhs2.shape[il2[i].group].range) laerror("incompatible index range in contractions");
#ifdef LA_TENSOR_INDEXPOSITION #ifdef LA_TENSOR_INDEXPOSITION
if(rhs1.shape[il1[i].group].upperindex ^ rhs2.shape[il2[i].group].upperindex == false) laerror("can contact only upper with lower index"); if((rhs1.shape[il1[i].group].upperindex ^ rhs2.shape[il2[i].group].upperindex) == false) laerror("can contact only upper with lower index");
#endif #endif
for(int j=0; j<i; ++j) if(il1[i]==il1[j]||il2[i]==il2[j]) laerror("repeated index in the list"); for(int j=0; j<i; ++j) if(il1[i]==il1[j]||il2[i]==il2[j]) laerror("repeated index in the list");
} }
@@ -1436,7 +1453,7 @@ for(int i=0; i<il1.size(); ++i)
if(shape[il1[i].group].offset != shape[il2[i].group].offset) laerror("incompatible index offset in innercontraction"); if(shape[il1[i].group].offset != shape[il2[i].group].offset) laerror("incompatible index offset in innercontraction");
if(shape[il1[i].group].range != shape[il2[i].group].range) laerror("incompatible index range in innercontraction"); if(shape[il1[i].group].range != shape[il2[i].group].range) laerror("incompatible index range in innercontraction");
#ifdef LA_TENSOR_INDEXPOSITION #ifdef LA_TENSOR_INDEXPOSITION
if(shape[il1[i].group].upperindex ^ shape[il2[i].group].upperindex == false) laerror("can contact only upper with lower index"); if((shape[il1[i].group].upperindex ^ shape[il2[i].group].upperindex) == false) laerror("can contact only upper with lower index");
#endif #endif
for(int j=0; j<i; ++j) if(il1[i]==il1[j]||il2[i]==il2[j]) laerror("repeated index in the innercontraction list"); for(int j=0; j<i; ++j) if(il1[i]==il1[j]||il2[i]==il2[j]) laerror("repeated index in the innercontraction list");
for(int j=0; j<il1.size(); ++j) if(il1[i]==il2[j]) laerror("repeated index between the innercontraction lists"); for(int j=0; j<il1.size(); ++j) if(il1[i]==il2[j]) laerror("repeated index between the innercontraction lists");
@@ -1701,13 +1718,13 @@ selected.clear();
for(int g=0; g<groups.size(); ++g) for(int g=0; g<groups.size(); ++g)
{ {
if(groups[g]<0||groups[g]>=shape.size()) laerror("illegal group number in merge_index_groups"); if(groups[g]<0||groups[g]>=shape.size()) laerror("illegal group number in merge_index_groups");
if(selected[g]) laerror("repeated group number in merge_index_groups"); if(selected[groups[g]]) laerror("repeated group number in merge_index_groups");
selected.set(g); selected.set(groups[g]);
p[gg++] = 1+groups[g]; p[1+gg++] = 1+groups[g];
} }
for(int g=0; g<shape.size(); ++g) for(int g=0; g<shape.size(); ++g)
if(!selected[g]) if(!selected[g])
p[gg++] = 1+g; p[1+gg++] = 1+g;
if(gg!=shape.size() || !p.is_valid()) laerror("internal error in merge_index_groups"); if(gg!=shape.size() || !p.is_valid()) laerror("internal error in merge_index_groups");
Tensor<T> r = permute_index_groups(p); //takes care of names permutation too Tensor<T> r = permute_index_groups(p); //takes care of names permutation too
@@ -2006,7 +2023,7 @@ for(int i=0; i<il.size(); ++i)
if(il[i].index<0||il[i].index>=shape[il[i].group].number) laerror("wrong index number in merge_indices"); if(il[i].index<0||il[i].index>=shape[il[i].group].number) laerror("wrong index number in merge_indices");
for(int j=0; j<i; ++j) if(il[i]==il[j]) laerror("repeated index in the list"); for(int j=0; j<i; ++j) if(il[i]==il[j]) laerror("repeated index in the list");
#ifdef LA_TENSOR_INDEXPOSITION #ifdef LA_TENSOR_INDEXPOSITION
if(shape[il[0].group].upperindex != shape[il[i].group].upperindex == false) laerror("can merge only within lower or upper separately"); if(shape[il[0].group].upperindex != shape[il[i].group].upperindex) laerror("can merge only within lower or upper separately");
#endif #endif
if(shape[il[0].group].range != shape[il[i].group].range) if(shape[il[0].group].range != shape[il[i].group].range)
{ {
@@ -2231,14 +2248,14 @@ if(!rhs1.is_uniquely_named()||!rhs2.is_uniquely_named()|| !is_uniquely_named())
//find contraction indices //find contraction indices
int nc=0; int nc=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]) if(rhs1.names[i]==rhs2.names[j])
{ {
//std::cout << "found contraction "<<nc<<" th. index = "<<rhs1.names[i]<<std::endl; //std::cout << "found contraction "<<nc<<" th. index = "<<rhs1.names[i]<<std::endl;
#ifdef LA_TENSOR_INDEXPOSITION #ifdef LA_TENSOR_INDEXPOSITION
int group1=indexposition(i,rhs1.shape); INDEX I1=indexposition(i,rhs1.shape);
int group2=indexposition(j,rhs2.shape); INDEX I2=indexposition(j,rhs2.shape);
if(rhs1.shape[group1].upperindex ^ rhs2.shape[group2].upperindex == false) laerror("can contract only upper with lower index"); if((rhs1.shape[I1.group].upperindex ^ rhs2.shape[I2.group].upperindex) == false) laerror("can contract only upper with lower index");
#endif #endif
++nc; ++nc;
} }
@@ -2386,6 +2403,7 @@ if(is_named())
int ii=0; int ii=0;
for(int i=0; i<g; ++i) ii+= shape[i].number; for(int i=0; i<g; ++i) ii+= shape[i].number;
NRVec<INDEXNAME> tmp=(names.subvector(ii,ii+p.size()-1)).permuted(p,inverse); NRVec<INDEXNAME> tmp=(names.subvector(ii,ii+p.size()-1)).permuted(p,inverse);
names.copyonwrite();
for(int i=0; i<p.size(); ++i) names[ii+i]=tmp[i]; for(int i=0; i<p.size(); ++i) names[ii+i]=tmp[i];
} }
} }
@@ -2519,7 +2537,7 @@ template<typename T>
void Tensor<T>::numpy_write(const char *name, const char *descr) const void Tensor<T>::numpy_write(const char *name, const char *descr) const
{ {
if(!is_flat()) laerror("numpy_write only for flat tensors"); if(!is_flat()) laerror("numpy_write only for flat tensors");
int fd=open(name,O_CREAT|O_LARGEFILE|O_RDWR,0777); int fd=open(name,O_CREAT|O_TRUNC|O_LARGEFILE|O_WRONLY,0666);
if(fd<0) laerror("cannot open in numpy_write"); if(fd<0) laerror("cannot open in numpy_write");
char magic[]="\x93NUMPY\x01\x00"; char magic[]="\x93NUMPY\x01\x00";
if(8!=write(fd,magic,8)) laerror("cannot write 1 in numpy_write"); if(8!=write(fd,magic,8)) laerror("cannot write 1 in numpy_write");
@@ -2571,16 +2589,18 @@ char readmagic[8];
if(8!=read(fd,readmagic,8)) laerror("cannot read 1 in numpy_read"); if(8!=read(fd,readmagic,8)) laerror("cannot read 1 in numpy_read");
if(memcmp(magic,readmagic,8)) laerror("bad magic or version != 1.0 in numpy_read"); if(memcmp(magic,readmagic,8)) laerror("bad magic or version != 1.0 in numpy_read");
int16_t header_len; uint16_t header_len; //unsigned in the npy format 1.0
if(2!=read(fd,&header_len,2)) laerror("cannot read 2 in numpy_read"); if(2!=read(fd,&header_len,2)) laerror("cannot read 2 in numpy_read");
char header[header_len]; char header[header_len+1];
if(header_len!=read(fd,header,header_len)) laerror("cannot read 3 in numpy_read"); if(header_len!=read(fd,header,header_len)) laerror("cannot read 3 in numpy_read");
header[header_len]=0; //for the string functions below
char *p=strstr(header,"fortran_order"); char *p=strstr(header,"fortran_order");
if(!p) laerror("cannot find fortran_order in numpy_read"); if(!p) laerror("cannot find fortran_order in numpy_read");
bool fortranorder=false; bool fortranorder=false;
p=strstr(p,": "); p=strstr(p,": ");
if(!p) laerror("cannot parse fortran_order in numpy_read");
if(strncmp(p+2,"False",5)) fortranorder=true; if(strncmp(p+2,"False",5)) fortranorder=true;
//std::cout <<"Fortran order "<<fortranorder<<std::endl; //std::cout <<"Fortran order "<<fortranorder<<std::endl;
tensor.h
+3 -3
View File
@@ -343,7 +343,7 @@ public:
int findflatindex(const INDEXNAME nam) const; int findflatindex(const INDEXNAME nam) const;
INDEX findindex(const INDEXNAME nam) const; INDEX findindex(const INDEXNAME nam) const;
NRVec<INDEX> findindexlist(const NRVec<INDEXNAME> &names) const; NRVec<INDEX> findindexlist(const NRVec<INDEXNAME> &names) const;
void renameindex(const INDEXNAME namfrom, const INDEXNAME nameto) {int i=findflatindex(namfrom); names[i]=nameto;}; void renameindex(const INDEXNAME namfrom, const INDEXNAME nameto) {int i=findflatindex(namfrom); names.copyonwrite(); names[i]=nameto;};
Tensor& operator+=(const Tensor &rhs) Tensor& operator+=(const Tensor &rhs)
{ {
@@ -400,8 +400,8 @@ public:
Tensor unwind_indices(const NRVec<INDEXNAME> &names) const {return unwind_indices(findindexlist(names));}; Tensor unwind_indices(const NRVec<INDEXNAME> &names) const {return unwind_indices(findindexlist(names));};
void addcontraction(const Tensor &rhs1, int group, int index, const Tensor &rhs2, int rhsgroup, int rhsindex, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false); //rhs1 will have more significant non-contracted indices in the result than rhs2 void addcontraction(const Tensor &rhs1, int group, int index, const Tensor &rhs2, int rhsgroup, int rhsindex, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false); //rhs1 will have more significant non-contracted indices in the result than rhs2
inline void addcontraction(const Tensor &rhs1, const INDEX &I1, const Tensor &rhs2, const INDEX &I2, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false) {addcontraction(rhs1, I1.group, I1.index, rhs2, I2.group, I2.index, alpha, beta, doresize, conjugate, conjugate);}; inline void addcontraction(const Tensor &rhs1, const INDEX &I1, const Tensor &rhs2, const INDEX &I2, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false) {addcontraction(rhs1, I1.group, I1.index, rhs2, I2.group, I2.index, alpha, beta, doresize, conjugate1, conjugate);};
inline void addcontraction(const Tensor &rhs1, const Tensor &rhs2, const INDEXNAME &iname, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false) {addcontraction(rhs1, rhs1.findindex(iname), rhs2, rhs2.findindex(iname), alpha, beta, doresize, conjugate, conjugate);}; inline void addcontraction(const Tensor &rhs1, const Tensor &rhs2, const INDEXNAME &iname, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false) {addcontraction(rhs1, rhs1.findindex(iname), rhs2, rhs2.findindex(iname), alpha, beta, doresize, conjugate1, conjugate);};
inline Tensor contraction(int group, int index, const Tensor &rhs, int rhsgroup, int rhsindex, T alpha=1, bool conjugate1=false, bool conjugate=false) const {Tensor<T> r; r.addcontraction(*this,group,index,rhs,rhsgroup,rhsindex,alpha,0,true, conjugate1, conjugate); return r; }; inline Tensor contraction(int group, int index, const Tensor &rhs, int rhsgroup, int rhsindex, T alpha=1, bool conjugate1=false, bool conjugate=false) const {Tensor<T> r; r.addcontraction(*this,group,index,rhs,rhsgroup,rhsindex,alpha,0,true, conjugate1, conjugate); return r; };
inline Tensor contraction(const INDEX &I, const Tensor &rhs, const INDEX &RHSI, T alpha=1, bool conjugate1=false, bool conjugate=false) const {return contraction(I.group,I.index, rhs, RHSI.group, RHSI.index,alpha, conjugate1, conjugate);}; inline Tensor contraction(const INDEX &I, const Tensor &rhs, const INDEX &RHSI, T alpha=1, bool conjugate1=false, bool conjugate=false) const {return contraction(I.group,I.index, rhs, RHSI.group, RHSI.index,alpha, conjugate1, conjugate);};
inline Tensor contraction(const Tensor &rhs, const INDEXNAME &iname, T alpha=1, bool conjugate1=false, bool conjugate=false) const {return contraction(findindex(iname),rhs,rhs.findindex(iname),alpha, conjugate1, conjugate);}; inline Tensor contraction(const Tensor &rhs, const INDEXNAME &iname, T alpha=1, bool conjugate1=false, bool conjugate=false) const {return contraction(findindex(iname),rhs,rhs.findindex(iname),alpha, conjugate1, conjugate);};