diff --git a/tensor.cc b/tensor.cc index 1951334..9fe4f49 100644 --- a/tensor.cc +++ b/tensor.cc @@ -16,25 +16,27 @@ along with this program. If not, see . */ +#include #include "tensor.h" #include "laerror.h" -#include +#include "qsort.h" +#include "miscfunc.h" namespace LA { template -LA_largeindex Tensor::index(const SUPERINDEX &I) +LA_largeindex Tensor::index(int *sign, const SUPERINDEX &I) { //check index structure and ranges -#ifndef DEBUG +#ifdef DEBUG if(I.size()!=shape.size()) laerror("mismatch in the number of tensor index groups"); -for(int i=0; i= shape[i].offset+shape[i].size()) + if(I[i][j] = shape[i].offset+shape[i].range) { std::cerr<<"error in index group no. "< +class Signedpointer +{ +T *ptr; +int sgn; +public: + Signedpointer(const T *p, int s) : ptr(p),sgn(s) {}; +//@@@@@@operations on singedpointer as LHS of the non-const tensor.operator() expressions +}; + + typedef int LA_index; typedef int LA_largeindex; @@ -42,7 +54,7 @@ typedef class indexgroup { int number; //number of indices int symmetry; //-1 0 or 1 LA_index offset; //indices start at -LA_index size; //indices span this range +LA_index range; //indices span this range } INDEXGROUP; typedef NRVec FLATINDEX; //all indices but in a single vector @@ -51,26 +63,31 @@ typedef NRVec > SUPERINDEX; //all indices in the INDEXGROUP stru template class Tensor { +//essential data NRVec shape; NRVec data; +//redundant data to facilitate efficient indexing + NRVec cumsizes; //cumulative sizes of symmetry index groups + private: - LA_largeindex index(const SUPERINDEX &I); //map the tensor indices to the position in data - LA_largeindex index(const FLATINDEX &I); //map the tensor indices to the position in data - LA_largeindex vindex(int i1,va_list args); //map list of indices to the position in data @@@must call va_end + LA_largeindex index(int *sign, const SUPERINDEX &I); //map the tensor indices to the position in data + LA_largeindex index(int *sign, const FLATINDEX &I); //map the tensor indices to the position in data + LA_largeindex vindex(int *sign, int i1, va_list args); //map list of indices to the position in data @@@must call va_end public: Tensor() {}; - Tensor(const NRVec &s) : shape(s), data((int)size()) {data.clear();}; - int rank() const; //is computed from shape - LA_largeindex size() const; //expensive, is computed from shape + Tensor(const NRVec &s) : shape(s), data((int)getsize()) {data.clear();}; + int getrank() const; //is computed from shape + LA_largeindex getsize(); //set redundant data and return total size void copyonwrite() {shape.copyonwrite(); data.copyonwrite();}; - inline T& operator()(const SUPERINDEX &I) {return data[index(I)];}; - inline const T& operator()(const SUPERINDEX &I) const {return data[index(I)];}; - inline T& operator()(const FLATINDEX &I) {return data[index(I)];}; - inline const T& operator()(const FLATINDEX &I) const {return data[index(I)];}; - inline T& operator()(int i1...) {va_list args; va_start(args,i1); return data[vindex(i1,args)];}; - inline const T& operator()(int i1...) const {va_list args; va_start(args,i1); return data[vindex(i1,args)];}; + inline Signedpointer operator()(const SUPERINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer(&data[i],sign);}; + inline const 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];}; + inline Signedpointer operator()(const FLATINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer(&data[i],sign);}; + inline const T& operator()(const FLATINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); if(sign==0) return 0; return sign>0 ?data[i] : -data[i];}; + inline Signedpointer operator()(int i1...) {va_list args; int sign; LA_largeindex i; va_start(args,i1); i= vindex(&sign, i1,args); return Signedpointer(&data[i],sign); }; + inline const T& operator()(int i1...) const {va_list args; ; int sign; LA_largeindex i; va_start(args,i1); i= vindex(&sign, i1,args); if(sign==0) return 0; return sign>0 ?data[i] : -data[i];}; + //NOTE: for sign==0 data[i] can be undefined pointer, avoid dereferencing it //@@@TODO - unwinding to full size in a specified index //@@@TODO - contractions - basic and efficient @@ -81,7 +98,7 @@ public: template -int Tensor:: rank() const +int Tensor:: getrank() const { int r=0; for(int i=0; i -LA_largeindex Tensor::size() const +LA_largeindex Tensor::getsize() { +cumsizes.resize(shape.size()); LA_largeindex s=1; for(int i=0; i