diff --git a/t.cc b/t.cc index f01b2e0..0a866c5 100644 --- a/t.cc +++ b/t.cc @@ -4025,7 +4025,7 @@ cout <<"Error = "<<(xx-x).norm()< xx = xf.merge_indices(il,sym); cout <<"xx = "<number; } @@ -52,7 +52,7 @@ cumsizes.resize(shape.size()); LA_largeindex s=1; for(int i=0; i *>(&shape))[i]; + const INDEXGROUP *sh = &(* const_cast *>(&shape))[i]; if(sh->number==0) laerror("empty index group"); if(sh->range==0) return 0; cumsizes[i]=s; @@ -396,7 +396,7 @@ template void loopingroups(Tensor &t, int ngroup, int igroup, T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, T *)) { LA_index istart,iend; -const indexgroup *sh = &(* const_cast *>(&t.shape))[ngroup]; +const INDEXGROUP *sh = &(* const_cast *>(&t.shape))[ngroup]; switch(sh->symmetry) { case 0: @@ -431,7 +431,7 @@ for(LA_index i = istart; i<=iend; ++i) if(newigroup<0) { --newngroup; - const indexgroup *sh2 = &(* const_cast *>(&t.shape))[newngroup]; + const INDEXGROUP *sh2 = &(* const_cast *>(&t.shape))[newngroup]; newigroup=sh2->number-1; } loopingroups(t,newngroup,newigroup,p,I,callback); @@ -446,13 +446,13 @@ void Tensor::loopover(void (*callback)(const SUPERINDEX &, T *)) SUPERINDEX I(shape.size()); for(int i=0; i *>(&shape))[i]; + const INDEXGROUP *sh = &(* const_cast *>(&shape))[i]; I[i].resize(sh->number); I[i] = sh->offset; } T *pp=&data[0]; int ss=shape.size()-1; -const indexgroup *sh = &(* const_cast *>(&shape))[ss]; +const INDEXGROUP *sh = &(* const_cast *>(&shape))[ss]; loopingroups(*this,ss,sh->number-1,&pp,I,callback); } @@ -461,7 +461,7 @@ template void constloopingroups(const Tensor &t, int ngroup, int igroup, const T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, const T *)) { LA_index istart,iend; -const indexgroup *sh = &t.shape[ngroup]; +const INDEXGROUP *sh = &t.shape[ngroup]; switch(sh->symmetry) { case 0: @@ -496,7 +496,7 @@ for(LA_index i = istart; i<=iend; ++i) if(newigroup<0) { --newngroup; - const indexgroup *sh2 = &(* const_cast *>(&t.shape))[newngroup]; + const INDEXGROUP *sh2 = &(* const_cast *>(&t.shape))[newngroup]; newigroup=sh2->number-1; } constloopingroups(t,newngroup,newigroup,p,I,callback); @@ -511,13 +511,13 @@ void Tensor::constloopover(void (*callback)(const SUPERINDEX &, const T *)) c SUPERINDEX I(shape.size()); for(int i=0; inumber); I[i] = sh->offset; } const T *pp=&data[0]; int ss=shape.size()-1; -const indexgroup *sh = &shape[ss]; +const INDEXGROUP *sh = &shape[ss]; constloopingroups(*this,ss,sh->number-1,&pp,I,callback); } @@ -664,7 +664,7 @@ help_t->data[target] = *v; //permutation of individual indices from permutation of index groups -NRPerm group2flat_perm(const NRVec &shape, const NRPerm &p) +NRPerm group2flat_perm(const NRVec &shape, const NRPerm &p) { int rank=0; for(int i=0; i Tensor Tensor::permute_index_groups(const NRPerm &p) const { //std::cout <<"permute_index_groups permutation = "< newshape=shape.permuted(p,true); +NRVec newshape=shape.permuted(p,true); //std::cout <<"permute_index_groups newshape = "< r(newshape); if(is_named()) @@ -718,10 +718,10 @@ for(int i=0; i &shape) +int flatposition(const INDEX &i, const NRVec &shape) {return flatposition(i.group,i.index,shape);} -int flatposition(int group, int index, const NRVec &shape) +int flatposition(int group, int index, const NRVec &shape) { int ii=0; for(int g=0; g &shape) +INDEX indexposition(int flatindex, const NRVec &shape) { INDEX I={0,0}; if(flatindex<0) laerror("illegal index in indexposition"); @@ -792,7 +792,7 @@ if(group==0 && index==0 && shape[0].symmetry==0) //formal split of 1 index witho 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); +NRVec newshape(shape.size()+1); newshape[0].number=1; newshape[0].symmetry=0; newshape[0].range=shape[group].range; @@ -914,7 +914,7 @@ if(group<0) newsize=rank(); else newsize=shape.size()+shape[group].number-1; //build new shape -NRVec newshape(newsize); +NRVec newshape(newsize); int gg=0; for(int g=0; g oldshape(shape); +NRVec oldshape(shape); oldshape.copyonwrite(); -NRVec newshape(shape.size()+nonsolo); +NRVec newshape(shape.size()+nonsolo); //first the unwound indices as solo groups for(int i=0; i 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); +NRVec newshape(u.shape.size()+rhsu.shape.size()-2); int ii=0; for(int i=1; i u = conjugate1? (rhs1.unwind_index(group,index)).conjugate() : r const Tensor rhsu = rhs.unwind_index(rhsgroup,rhsindex); -NRVec newshape(u.shape.size()+rhsu.shape.size()-2); +NRVec newshape(u.shape.size()+rhsu.shape.size()-2); int ii=0; for(int i=1; i(nn,mm,kk,&data[0],&u.data[0], &rhsu.data[0],alpha,beta,conjugate); template void Tensor::addcontractions(const Tensor &rhs1, const INDEXLIST &il1, const Tensor &rhs2, const INDEXLIST &il2, T alpha, T beta, bool doresize, bool conjugate1, bool conjugate2) { -if(il1.size()==0) laerror("empty contraction - outer product not implemented"); if(il1.size()!=il2.size()) laerror("mismatch in index lists in addcontractions"); +if(il1.size()==0) //add just outer product + { + if(beta!=(T)0) *this *= beta; else this->clear(); + *this += (conjugate1?rhs1.conjugate():rhs1) * (conjugate2?rhs2.conjugate():rhs2) * alpha; + return; + } + for(int i=0; i=rhs1.shape.size()) laerror("wrong group1 number in contractions"); @@ -1250,7 +1256,7 @@ const Tensor u = conjugate1? (rhs1.unwind_indices(il1)).conjugateme() : rhs1. const Tensor rhsu = rhs2.unwind_indices(il2); -NRVec newshape(u.shape.size()+rhsu.shape.size()-2*il1.size()); +NRVec newshape(u.shape.size()+rhsu.shape.size()-2*il1.size()); int ii=0; for(int i=il1.size(); i(nn,mm,kk,&data[0],&u.data[0], &rhsu.data[0],alpha,beta,conjugate2) } +int help_tdim; +const SUPERINDEX *help_dummyindex; + +template +static void innercontraction_callback(const SUPERINDEX &I, T *r) +{ +//compute address like in operarator() but we need true address +//this could be done more efficiently if needed, avoiding explicit SUPERINDEX J +SUPERINDEX J = help_dummyindex->concat(I); +int sign; +LA_largeindex ii=help_tt->index(&sign,J); +if(sign<0) laerror("internal error in innercontraction"); +const T *matrix00 = &help_tt->data[ii]; +*r = 0; +for(int i=0; i +Tensor Tensor::innercontraction(const INDEXLIST &il1, const INDEXLIST &il2) const +{ +if(il1.size()!=il2.size()) laerror("mismatch in index lists in innercontraction"); +if(il1.size()==0) return *this; +for(int i=0; i=shape.size()) laerror("wrong group1 number in innercontraction"); + if(il2[i].group<0||il2[i].group>=shape.size()) laerror("wrong group2 number in innercontraction"); + if(il1[i].index<0||il1[i].index>=shape[il1[i].group].number) laerror("wrong index1 number in innerconntraction"); + if(il2[i].index<0||il2[i].index>=shape[il2[i].group].number) laerror("wrong index2 number in innerconntraction"); + 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"); +#ifdef LA_TENSOR_INDEXPOSITION + if(shape[il1[i].group].upperindex ^ shape[il2[i].group].upperindex == false) laerror("can contact only upper with lower index"); +#endif + for(int j=0; j(s); + } + +//result is tensor +NRVec newshape = (u.shape).subvector(il.size(),(u.shape).size()-1); +//std::cout <<"new shape = "< = &u; +help_tdim = tracedim; +SUPERINDEX dummyindex(il.size()); +for(int i=0; i @@ -1387,12 +1465,12 @@ loopover(permutationalgebra_callback2); template void Tensor::split_index_group(int group) { -const indexgroup *sh = &(* const_cast *>(&shape))[0]; +const INDEXGROUP *sh = &(* const_cast *>(&shape))[0]; if(group<0||group >= shape.size()) laerror("illegal index group number"); if(sh[group].number==1) return; //nothing to split if(sh[group].symmetry!=0) laerror("only non-symmetric index group can be splitted, use flatten instead"); -NRVec newshape(shape.size()+sh[group].number-1); +NRVec newshape(shape.size()+sh[group].number-1); int gg=0; for(int g=0; g= 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); +NRVec newshape(shape.size()+1); int gg=0; for(int g=0; g newshape(shape.size()-(groupto-groupfrom+1)+1); +NRVec newshape(shape.size()-(groupto-groupfrom+1)+1); for(int g=0; g<=groupfrom; ++g) newshape[g]=shape[g]; newshape[groupfrom].number=newnumber; for(int g=groupfrom+1; g um; - NRVec ushape; + NRVec ushape; { Tensor uu=unwind_index(I); ushape=uu.shape; //ushape.copyonwrite(); should not be needed @@ -1682,7 +1760,7 @@ if(samegroup && isordered && il.size()==shape[il[0].group].number) return unwind //calculate new shape and flat index permutation -NRVec workshape(shape); +NRVec workshape(shape); workshape.copyonwrite(); NRPerm basicperm(rank()); @@ -1698,7 +1776,7 @@ for(int i=0; i0) ++newshapesize; //this group survived index removal -NRVec newshape(newshapesize); +NRVec newshape(newshapesize); newshape[0].number=il.size(); newshape[0].symmetry=sym; newshape[0].offset=shape[il[0].group].offset; @@ -1753,7 +1831,7 @@ return r; template void Tensor::canonicalize_shape() { -const indexgroup *sh = &(* const_cast *>(&shape))[0]; +const INDEXGROUP *sh = &(* const_cast *>(&shape))[0]; for(int i=0; i { }; -typedef class indexgroup { +typedef class INDEXGROUP { public: int number; //number of indices int symmetry; //-1 0 or 1, later 2 for hermitian and -2 for antihermitian? - would need change in operator() and Signedpointer @@ -120,12 +117,12 @@ LA_index range; //indices span this range bool upperindex; #endif - bool operator==(const indexgroup &rhs) const {return number==rhs.number && symmetry==rhs.symmetry && offset==rhs.offset && range==rhs.range + bool operator==(const INDEXGROUP &rhs) const {return number==rhs.number && symmetry==rhs.symmetry && offset==rhs.offset && range==rhs.range #ifdef LA_TENSOR_INDEXPOSITION && upperindex == rhs.upperindex #endif ;}; - inline bool operator!=(const indexgroup &rhs) const {return !((*this)==rhs);}; + inline bool operator!=(const INDEXGROUP &rhs) const {return !((*this)==rhs);}; } INDEXGROUP; @@ -139,17 +136,17 @@ std::istream & operator>>(std::istream &s, INDEXNAME &x); template<> -class LA_traits { +class LA_traits { public: static bool is_plaindata() {return true;}; - static void copyonwrite(indexgroup& x) {}; + static void copyonwrite(INDEXGROUP& x) {}; typedef INDEXGROUP normtype; - static inline int gencmp(const indexgroup *a, const indexgroup *b, int n) {return memcmp(a,b,n*sizeof(indexgroup));}; - static inline void put(int fd, const indexgroup &x, bool dimensions=1) {if(sizeof(indexgroup)!=write(fd,&x,sizeof(indexgroup))) laerror("write error 1 in indexgroup put"); } - static inline void multiput(int nn, int fd, const indexgroup *x, bool dimensions=1) {if(nn*sizeof(indexgroup)!=write(fd,x,nn*sizeof(indexgroup))) laerror("write error 1 in indexgroup multiiput"); } - static inline void get(int fd, indexgroup &x, bool dimensions=1) {if(sizeof(indexgroup)!=read(fd,&x,sizeof(indexgroup))) laerror("read error 1 in indexgroup get");} - static inline void multiget(int nn, int fd, indexgroup *x, bool dimensions=1) {if(nn*sizeof(indexgroup)!=read(fd,x,nn*sizeof(indexgroup))) laerror("read error 1 in indexgroup get");} - static inline void clear(indexgroup *dest, size_t n) {memset(dest,0,n*sizeof(indexgroup));} + static inline int gencmp(const INDEXGROUP *a, const INDEXGROUP *b, int n) {return memcmp(a,b,n*sizeof(INDEXGROUP));}; + static inline void put(int fd, const INDEXGROUP &x, bool dimensions=1) {if(sizeof(INDEXGROUP)!=write(fd,&x,sizeof(INDEXGROUP))) laerror("write error 1 in INDEXGROUP put"); } + static inline void multiput(int nn, int fd, const INDEXGROUP *x, bool dimensions=1) {if(nn*sizeof(INDEXGROUP)!=write(fd,x,nn*sizeof(INDEXGROUP))) laerror("write error 1 in INDEXGROUP multiiput"); } + static inline void get(int fd, INDEXGROUP &x, bool dimensions=1) {if(sizeof(INDEXGROUP)!=read(fd,&x,sizeof(INDEXGROUP))) laerror("read error 1 in INDEXGROUP get");} + static inline void multiget(int nn, int fd, INDEXGROUP *x, bool dimensions=1) {if(nn*sizeof(INDEXGROUP)!=read(fd,x,nn*sizeof(INDEXGROUP))) laerror("read error 1 in INDEXGROUP get");} + static inline void clear(INDEXGROUP *dest, size_t n) {memset(dest,0,n*sizeof(INDEXGROUP));} }; @@ -168,16 +165,36 @@ std::ostream & operator<<(std::ostream &s, const INDEX &x); std::istream & operator>>(std::istream &s, INDEX &x); -int flatposition(int group, int index, const NRVec &shape); -int flatposition(const INDEX &i, const NRVec &shape); //position of that index in FLATINDEX -INDEX indexposition(int flatindex, const NRVec &shape); //inverse to flatposition +template<> +class LA_traits { + public: + static bool is_plaindata() {return true;}; + static void copyonwrite(INDEX& x) {}; + typedef INDEX normtype; + static inline int gencmp(const INDEX *a, const INDEX *b, int n) {return memcmp(a,b,n*sizeof(INDEX));}; + static inline void put(int fd, const INDEX &x, bool dimensions=1) {if(sizeof(INDEX)!=write(fd,&x,sizeof(INDEX))) laerror("write error 1 in INDEX put"); } + static inline void multiput(int nn, int fd, const INDEX *x, bool dimensions=1) {if(nn*sizeof(INDEX)!=write(fd,x,nn*sizeof(INDEX))) laerror("write error 1 in INDEX multiiput"); } + static inline void get(int fd, INDEX &x, bool dimensions=1) {if(sizeof(INDEX)!=read(fd,&x,sizeof(INDEX))) laerror("read error 1 in INDEX get");} + static inline void multiget(int nn, int fd, INDEX *x, bool dimensions=1) {if(nn*sizeof(INDEX)!=read(fd,x,nn*sizeof(INDEX))) laerror("read error 1 in INDEX get");} + static inline void clear(INDEX *dest, size_t n) {memset(dest,0,n*sizeof(INDEX));} +}; + + + + + + + +int flatposition(int group, int index, const NRVec &shape); +int flatposition(const INDEX &i, const NRVec &shape); //position of that index in FLATINDEX +INDEX indexposition(int flatindex, const NRVec &shape); //inverse to flatposition FLATINDEX superindex2flat(const SUPERINDEX &I); template class Tensor { public: - NRVec shape; + NRVec shape; NRVec data; int myrank; NRVec groupsizes; //group sizes of symmetry index groups (a function of shape but precomputed for efficiency) @@ -193,15 +210,15 @@ public: //constructors Tensor() : myrank(-1) {}; explicit Tensor(const T &x) : myrank(0), data(1) {data[0]=x;}; //scalar - Tensor(const NRVec &s) : shape(s) { data.resize(calcsize()); calcrank(); canonicalize_shape();}; //general tensor - Tensor(const NRVec &s, const NRVec &newnames) : shape(s), names(newnames) { data.resize(calcsize()); calcrank(); canonicalize_shape(); if(names.size()!=myrank && names.size()!=0) laerror("bad number of index names");}; //general tensor - Tensor(const NRVec &s, const NRVec &mydata) : shape(s) { LA_largeindex dsize=calcsize(); calcrank(); canonicalize_shape(); if(mydata.size()!=dsize) laerror("inconsistent data size with shape"); data=mydata;} - Tensor(const NRVec &s, const NRVec &mydata, const NRVec &newnames) : shape(s), names(newnames) { LA_largeindex dsize=calcsize(); calcrank(); canonicalize_shape(); if(mydata.size()!=dsize) laerror("inconsistent data size with shape"); data=mydata; if(names.size()!=myrank && names.size()!=0) laerror("bad number of index names");} - Tensor(const indexgroup &g) {shape.resize(1); shape[0]=g; data.resize(calcsize()); calcrank(); canonicalize_shape();}; //tensor with a single index group - Tensor(const indexgroup &g, const NRVec &newnames) : names(newnames) {shape.resize(1); shape[0]=g; data.resize(calcsize()); calcrank(); canonicalize_shape(); if(names.size()!=myrank && names.size()!=0) laerror("bad number of index names");}; //tensor with a single index group + Tensor(const NRVec &s) : shape(s) { data.resize(calcsize()); calcrank(); canonicalize_shape();}; //general tensor + Tensor(const NRVec &s, const NRVec &newnames) : shape(s), names(newnames) { data.resize(calcsize()); calcrank(); canonicalize_shape(); if(names.size()!=myrank && names.size()!=0) laerror("bad number of index names");}; //general tensor + Tensor(const NRVec &s, const NRVec &mydata) : shape(s) { LA_largeindex dsize=calcsize(); calcrank(); canonicalize_shape(); if(mydata.size()!=dsize) laerror("inconsistent data size with shape"); data=mydata;} + Tensor(const NRVec &s, const NRVec &mydata, const NRVec &newnames) : shape(s), names(newnames) { LA_largeindex dsize=calcsize(); calcrank(); canonicalize_shape(); if(mydata.size()!=dsize) laerror("inconsistent data size with shape"); data=mydata; if(names.size()!=myrank && names.size()!=0) laerror("bad number of index names");} + Tensor(const INDEXGROUP &g) {shape.resize(1); shape[0]=g; data.resize(calcsize()); calcrank(); canonicalize_shape();}; //tensor with a single index group + Tensor(const INDEXGROUP &g, const NRVec &newnames) : names(newnames) {shape.resize(1); shape[0]=g; data.resize(calcsize()); calcrank(); canonicalize_shape(); if(names.size()!=myrank && names.size()!=0) laerror("bad number of index names");}; //tensor with a single index group Tensor(const Tensor &rhs): myrank(rhs.myrank), shape(rhs.shape), groupsizes(rhs.groupsizes), cumsizes(rhs.cumsizes), data(rhs.data), names(rhs.names) {}; - Tensor(int xrank, const NRVec &xshape, const NRVec &xgroupsizes, const NRVec xcumsizes, const NRVec &xdata) : myrank(xrank), shape(xshape), groupsizes(xgroupsizes), cumsizes(xcumsizes), data(xdata) {}; - Tensor(int xrank, const NRVec &xshape, const NRVec &xgroupsizes, const NRVec xcumsizes, const NRVec &xdata, const NRVec &xnames) : myrank(xrank), shape(xshape), groupsizes(xgroupsizes), cumsizes(xcumsizes), data(xdata), names(xnames) {}; + Tensor(int xrank, const NRVec &xshape, const NRVec &xgroupsizes, const NRVec xcumsizes, const NRVec &xdata) : myrank(xrank), shape(xshape), groupsizes(xgroupsizes), cumsizes(xcumsizes), data(xdata) {}; + Tensor(int xrank, const NRVec &xshape, const NRVec &xgroupsizes, const NRVec xcumsizes, const NRVec &xdata, const NRVec &xnames) : myrank(xrank), shape(xshape), groupsizes(xgroupsizes), cumsizes(xcumsizes), data(xdata), names(xnames) {}; //conversions from/to matrix and vector explicit Tensor(const NRVec &x); @@ -221,7 +238,7 @@ public: LA_largeindex calcsize(); //set redundant data and return total size LA_largeindex size() const {return data.size();}; void copyonwrite() {shape.copyonwrite(); groupsizes.copyonwrite(); cumsizes.copyonwrite(); data.copyonwrite(); names.copyonwrite();}; - void resize(const NRVec &s) {shape=s; data.resize(calcsize()); calcrank(); names.clear();}; + void resize(const NRVec &s) {shape=s; data.resize(calcsize()); calcrank(); names.clear();}; void deallocate() {data.resize(0); shape.resize(0); groupsizes.resize(0); cumsizes.resize(0); names.resize(0);}; inline Signedpointer lhs(const SUPERINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer(&data[i],sign);}; inline T operator()(const SUPERINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); if(sign==0) return 0; return sign>0 ?data[i] : -data[i];}; @@ -295,7 +312,6 @@ public: 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 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 Tensor contraction(int group, int index, const Tensor &rhs, int rhsgroup, int rhsindex, T alpha=1, bool conjugate1=false, bool conjugate=false) const {Tensor 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 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);}; @@ -308,6 +324,11 @@ 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 r; r.addgroupcontraction(*this,group,rhs,rhsgroup,alpha,0,true, conjugate1, conjugate); return r; }; + 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 &nl1, const NRVec &nl2) const {return innercontraction(findindexlist(nl1),findindexlist(nl2));}; + Tensor innercontraction(const INDEXNAME &n1, const INDEXNAME &n2) const {return innercontraction(findindex(n1),findindex(n2));}; + 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)) @@ -345,7 +366,7 @@ void tensor2fourindex(const Tensor &t, fourindex_dense &f); template Tensor fourindex2tensor(const fourindex_dense &f) { -NRVec shape(4); +NRVec shape(4); int n=f.nbas(); for(int i=0; i<4; ++i) { @@ -378,7 +399,7 @@ f=fourindex_dense(range,NRMat(t.data,range*range,range*range) template Tensor fourindex2tensor(const fourindex_dense &f) { -NRVec shape(2); +NRVec shape(2); int n=f.nbas(); for(int i=0; i<2; ++i) { @@ -412,7 +433,7 @@ f=fourindex_dense(NRMat(t.data,range*(range+1) template Tensor fourindex2tensor(const fourindex_dense &f) { -NRVec shape(2); +NRVec shape(2); int n=f.nbas(); for(int i=0; i<2; ++i) { @@ -448,7 +469,7 @@ f=fourindex_dense(NRSMat(mat)); //symmetrize mat template Tensor fourindex2tensor(const fourindex_dense &f) { -NRVec shape(4); +NRVec shape(4); for(int i=0; i<4; ++i) { shape[i].number=1; @@ -487,7 +508,7 @@ f=fourindex_dense(noca,nocb,nvra,nvrb,mat); template Tensor fourindex2tensor(const fourindex_dense &f) { -NRVec shape(2); +NRVec shape(2); for(int i=0; i<2; ++i) { shape[i].number=2; @@ -522,7 +543,7 @@ f=fourindex_dense(nocc,nvrt,mat); template Tensor fourindex2tensor(const fourindex_dense &f) { -NRVec shape(4); +NRVec shape(4); int n=f.nbas(); for(int i=0; i<4; ++i) { @@ -557,7 +578,7 @@ f=fourindex_dense(range,NRSMat(mat)); template Tensor fourindex2tensor(const fourindex_dense &f) { -NRVec shape(2); +NRVec shape(2); int n=f.nbas(); for(int i=0; i<2; ++i) { diff --git a/vec.cc b/vec.cc index c4434cb..3eeec6c 100644 --- a/vec.cc +++ b/vec.cc @@ -863,103 +863,6 @@ return -1; -/***************************************************************************//** - * extract block subvector - * @param[in] from starting position - * @param[in] to final position - * @return extracted block subvector - ******************************************************************************/ -template -const NRVec NRVec::subvector(const int from, const int to) const -{ -#ifdef DEBUG - if(from<0 || from>=nn|| to<0 || to>=nn || from>to){ - laerror("invalid subvector specification"); - } -#endif -//trivial case of whole vector for efficiency - if(from==0 && to == nn-1) return *this; - - const int n = to - from + 1; - NRVec r(n, getlocation()); - if(!LA_traits::is_plaindata()) laerror("only implemented for plain data"); - -#ifdef CUDALA - if(location == cpu){ -#endif - memcpy(r.v, v+from, n*sizeof(T)); -#ifdef CUDALA - }else{ - if(sizeof(T)%sizeof(float) != 0) laerror("cpu memcpy alignment problem"); - cublasScopy(n*sizeof(T)/sizeof(float), (const float *)(v+from), 1, (float*)r.v, 1); - TEST_CUBLAS("cublasScopy"); - } -#endif - return r; -} - - -template -const NRVec NRVec::subvector(const NRVec &selection) const -{ - NOT_GPU(*this); - const int n = selection.size(); - NRVec r(n); - - for(int i=0; i=nn) laerror("bad row index in subvector"); - r[i] = (*this)[ii]; - } -return r; - -} - -/***************************************************************************//** - * places given vector as subvector at given position - * @param[in] from coordinate - * @param[in] rhs input vector - ******************************************************************************/ -template -void NRVec::storesubvector(const int from, const NRVec &rhs) -{ - const int to = from + rhs.size() - 1; -#ifdef DEBUG - if(from<0 || from>=nn || to>=nn) laerror("bad indices in storesubvector"); -#endif - SAME_LOC(*this, rhs); - if(!LA_traits::is_plaindata()) laerror("only implemented for plain data"); - - #ifdef CUDALA - if(location == cpu){ - #endif - memcpy(v+from, rhs.v, rhs.size()*sizeof(T)); - - #ifdef CUDALA - }else{ - if(sizeof(T)%sizeof(float) != 0) laerror("cpu memcpy alignment problem"); - cublasScopy(rhs.size()*sizeof(T)/sizeof(float), (const float *) (rhs.v), 1, (float *)(v + from), 1); - } - #endif -} - - -template -void NRVec::storesubvector(const NRVec &selection, const NRVec &rhs) -{ - NOT_GPU(*this); - const int n = selection.size(); - if(n!=rhs.size()) laerror("size mismatch in storesubvector"); - - for(int i=0; i=nn) laerror("bad index in storesubvector"); - (*this)[ii] = rhs[i]; - } -} - /***************************************************************************//** * conjugate this general vector * @return reference to the (unmodified) matrix diff --git a/vec.h b/vec.h index 6b55b6c..8c4a6e3 100644 --- a/vec.h +++ b/vec.h @@ -2115,6 +2115,104 @@ for(int i=0; i +const NRVec NRVec::subvector(const int from, const int to) const +{ +#ifdef DEBUG + if(from<0 || from>=nn|| to<0 || to>=nn || from>to){ + laerror("invalid subvector specification"); + } +#endif +//trivial case of whole vector for efficiency + if(from==0 && to == nn-1) return *this; + + const int n = to - from + 1; + NRVec r(n, getlocation()); + if(!LA_traits::is_plaindata()) laerror("only implemented for plain data"); + +#ifdef CUDALA + if(location == cpu){ +#endif + memcpy(r.v, v+from, n*sizeof(T)); +#ifdef CUDALA + }else{ + if(sizeof(T)%sizeof(float) != 0) laerror("cpu memcpy alignment problem"); + cublasScopy(n*sizeof(T)/sizeof(float), (const float *)(v+from), 1, (float*)r.v, 1); + TEST_CUBLAS("cublasScopy"); + } +#endif + return r; +} + + +template +const NRVec NRVec::subvector(const NRVec &selection) const +{ + NOT_GPU(*this); + const int n = selection.size(); + NRVec r(n); + + for(int i=0; i=nn) laerror("bad row index in subvector"); + r[i] = (*this)[ii]; + } +return r; + +} + +/***************************************************************************//** + * places given vector as subvector at given position + * @param[in] from coordinate + * @param[in] rhs input vector + ******************************************************************************/ +template +void NRVec::storesubvector(const int from, const NRVec &rhs) +{ + const int to = from + rhs.size() - 1; +#ifdef DEBUG + if(from<0 || from>=nn || to>=nn) laerror("bad indices in storesubvector"); +#endif + SAME_LOC(*this, rhs); + if(!LA_traits::is_plaindata()) laerror("only implemented for plain data"); + + #ifdef CUDALA + if(location == cpu){ + #endif + memcpy(v+from, rhs.v, rhs.size()*sizeof(T)); + + #ifdef CUDALA + }else{ + if(sizeof(T)%sizeof(float) != 0) laerror("cpu memcpy alignment problem"); + cublasScopy(rhs.size()*sizeof(T)/sizeof(float), (const float *) (rhs.v), 1, (float *)(v + from), 1); + } + #endif +} + + +template +void NRVec::storesubvector(const NRVec &selection, const NRVec &rhs) +{ + NOT_GPU(*this); + const int n = selection.size(); + if(n!=rhs.size()) laerror("size mismatch in storesubvector"); + + for(int i=0; i=nn) laerror("bad index in storesubvector"); + (*this)[ii] = rhs[i]; + } +} + }//namespace