From 0ff55b66bb694c6b0d646b7030f95f59ab22f534 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Wed, 10 Apr 2024 18:28:50 +0200 Subject: [PATCH] working on tensor : stream I/O --- t.cc | 6 +- tensor.cc | 179 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ tensor.h | 58 ++++++++++++++---- 3 files changed, 231 insertions(+), 12 deletions(-) diff --git a/t.cc b/t.cc index a4c9d55..e99d72e 100644 --- a/t.cc +++ b/t.cc @@ -3206,9 +3206,10 @@ INDEXGROUP g; g.number=3; g.symmetry= -1; g.offset=1; -g.range=10; +g.range=5; Tensor epsilon(g); +epsilon.clear(); cout < +Tensor::Tensor(const NRVec &x) +: data(x) +{ +myrank=1; +shape.resize(1); +shape[0].number=1; +shape[0].symmetry=0; +#ifndef LA_TENSOR_ZERO_OFFSET +shape[0].offset=0; +#endif +shape[0].range=x.size(); +calcsize(); +} + +template +Tensor::Tensor(const NRMat &x) +: data(&x(0,0),x.nrows()*x.ncols()) +{ +myrank=2; +if(x.nrows()==x.ncols()) + { + shape.resize(1); + shape[0].number=2; + shape[0].symmetry=0; +#ifndef LA_TENSOR_ZERO_OFFSET + shape[0].offset=0; +#endif + shape[0].range=x.nrows(); + } +else + { + shape.resize(2); + shape[0].number=1; shape[1].number=1; + shape[0].symmetry=0; shape[1].symmetry=0; +#ifndef LA_TENSOR_ZERO_OFFSET + shape[0].offset=0; shape[1].offset=0; +#endif + shape[0].range=x.ncols(); + shape[1].range=x.nrows(); + } +calcsize(); +} + + +template +Tensor::Tensor(const NRSMat &x) +: data(NRVec(x)) +{ +myrank=2; +shape.resize(1); +shape[0].number=2; +shape[0].symmetry=1; +#ifndef LA_TENSOR_ZERO_OFFSET + shape[0].offset=0; +#endif +shape[0].range=x.nrows(); +calcsize(); +} + + +template +void loopingroups(Tensor &t, int ngroup, int igroup, T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, T *)) +{ +LA_index istart,iend; +switch(t.shape[ngroup].symmetry) + { + case 0: + istart= t.shape[ngroup].offset; + iend= t.shape[ngroup].offset+t.shape[ngroup].range-1; + break; + case 1: + istart= t.shape[ngroup].offset; + if(igroup==t.shape[ngroup].number-1) iend= t.shape[ngroup].offset+t.shape[ngroup].range-1; + else iend = I[ngroup][igroup+1]; + break; + case -1: + istart= t.shape[ngroup].offset + igroup; + if(igroup==t.shape[ngroup].number-1) iend= t.shape[ngroup].offset+t.shape[ngroup].range-1; + else iend = I[ngroup][igroup+1]-1; + break; + } + +for(LA_index i = istart; i<=iend; ++i) + { + I[ngroup][igroup]=i; + if(ngroup==0 && igroup==0) + { + int sign; + //std::cout <<"TEST "< +void Tensor::loopover(void (*callback)(const SUPERINDEX &, T *)) +{ +SUPERINDEX I(shape.size()); +for(int i=0; i +static void outputcallback(const SUPERINDEX &I, T *v) +{ +//print indices flat +for(int i=0; i>(std::istream &s, INDEXGROUP &x) +{ +s>>x.number>>x.symmetry; +#ifndef LA_TENSOR_ZERO_OFFSET +s>>x.offset; +#endif +s>>x.range; +return s; +} + + + +template +std::ostream & operator<<(std::ostream &s, const Tensor &x) +{ +s< *>(&x)->loopover(&outputcallback); +return s; +} + +template +std::istream & operator>>(std::istream &s, Tensor &x) +{ +s>>x.shape; +x.data.resize(x.calcsize()); x.calcrank(); +FLATINDEX I(x.rank()); +for(LA_largeindex i=0; i>I[j]; + T val; s>>val; + x.lhs(I) = val; + } +return s; +} + template class Tensor; template class Tensor >; +template std::ostream & operator<<(std::ostream &s, const Tensor &x); +template std::ostream & operator<<(std::ostream &s, const Tensor > &x); +template std::istream & operator>>(std::istream &s, Tensor &x); +template std::istream & operator>>(std::istream &s, Tensor > &x); }//namespace diff --git a/tensor.h b/tensor.h index 5ebd25d..884b0a7 100644 --- a/tensor.h +++ b/tensor.h @@ -19,6 +19,8 @@ //a simple tensor class with arbitrary symmetry of index subgroups //stored in an efficient way +//each index group has a specific symmetry (nosym,sym,antisym) +//additional symmetry between index groups (like in 2-electron integrals) is not supported directly, you would need to nest the class to Tensor > //presently only a rudimentary implementation //presently limited to 2G data size due to NRVec - maybe use a typedef LA_index //to uint64_t in the future in vector and matrix classes @@ -30,6 +32,8 @@ #include #include #include "vec.h" +#include "mat.h" +#include "smat.h" #include "miscfunc.h" @@ -70,6 +74,12 @@ LA_index range; //indices span this range inline bool operator!=(const indexgroup &rhs) const {return !((*this)==rhs);}; } INDEXGROUP; + +std::ostream & operator<<(std::ostream &s, const INDEXGROUP &x); +std::istream & operator>>(std::istream &s, INDEXGROUP &x); + + + template<> class LA_traits { public: @@ -91,11 +101,13 @@ typedef NRVec > SUPERINDEX; //all indices in the INDEXGROUP stru template class Tensor { - int myrank; +public: NRVec shape; + NRVec data; +private: + int myrank; NRVec groupsizes; //group sizes of symmetry index groups (a function of shape but precomputed for efficiency) NRVec cumsizes; //cumulative sizes of symmetry index groups (a function of shape but precomputed for efficiency) - NRVec data; public: LA_largeindex index(int *sign, const SUPERINDEX &I) const; //map the tensor indices to the position in data @@ -105,10 +117,13 @@ public: //constructors Tensor() : myrank(0) {}; - Tensor(const NRVec &s) : shape(s), data((int)calcsize()), myrank(calcrank()) {}; //general tensor - Tensor(const indexgroup &g) {shape.resize(1); shape[0]=g; data.resize(calcsize()); myrank=calcrank();}; //tensor with a single index group + Tensor(const NRVec &s) : shape(s), data((int)calcsize()) {calcrank();}; //general tensor + Tensor(const indexgroup &g) {shape.resize(1); shape[0]=g; data.resize(calcsize()); calcrank();}; //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) {}; 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) {}; + explicit Tensor(const NRVec &x); + explicit Tensor(const NRMat &x); + explicit Tensor(const NRSMat &x); void clear() {data.clear();}; int rank() const {return myrank;}; @@ -131,8 +146,23 @@ public: inline Tensor operator/(const T &a) const {Tensor r(*this); r /=a; return r;}; - inline Tensor& operator+=(const Tensor &rhs) {if(shape!=rhs.shape) laerror("incompatible tensors for operation"); data+=rhs.data; return *this;} - inline Tensor& operator-=(const Tensor &rhs) {if(shape!=rhs.shape) laerror("incompatible tensors for operation"); data-=rhs.data; return *this;} + inline Tensor& operator+=(const Tensor &rhs) + { +#ifdef DEBUG + if(shape!=rhs.shape) laerror("incompatible tensors for operation"); +#endif + data+=rhs.data; + return *this; + } + + inline Tensor& operator-=(const Tensor &rhs) + { +#ifdef DEBUG + if(shape!=rhs.shape) laerror("incompatible tensors for operation"); +#endif + data-=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;}; @@ -143,19 +173,25 @@ public: inline void randomize(const typename LA_traits::normtype &x) {data.randomize(x);}; + void loopover(void (*callback)(const SUPERINDEX &, T *)); + //@@@TODO - unwinding to full size in a specified index //@@@contraction by a whole index group or by individual single index //@@@TODO - contractions - basic and efficient? first contraction in a single index; between a given group+index in group at each tensor - //@@@symmetrize a group, antisymmetrize a group, expand a (anti)symmetric grtoup - obecne symmetry change krome +1 na -1 vse mozne //@@@outer product and product with a contraction - //@@@@permuteindexgroups - //@@@@@@explicit constructors from vec mat smat and dense fourindex - //@@@@@@ dvojite rekurzivni loopover s callbackem - nebo iterator s funkci next??? - //@@@@@@ stream i/o na zaklade tohoto + //@@@@symmetrize a group, antisymmetrize a group, expand a (anti)symmetric grtoup - obecne symmetry change krome +1 na -1 vse mozne + //@@@@@@permute index groups }; +template +std::ostream & operator<<(std::ostream &s, const Tensor &x); + +template +std::istream & operator>>(std::istream &s, Tensor &x); + +