diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/normal_moments.h b/mex/sources/korderpert/Dyn_pp/tl/cc/normal_moments.h index f179f78a349ef31b2c146fc6f2652bdb8f4559e7..8a2cf91946c9b140fdeb5453cb16db3b2f0cf33a 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/normal_moments.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/normal_moments.h @@ -9,10 +9,10 @@ class UNormalMoments:public TensorContainer<URSingleTensor> { public: -UNormalMoments(int maxdim,const TwoDMatrix&v); + UNormalMoments(int maxdim,const TwoDMatrix&v); private: -void generateMoments(int maxdim,const TwoDMatrix&v); -static bool selectEquiv(const Equivalence&e); + void generateMoments(int maxdim,const TwoDMatrix&v); + static bool selectEquiv(const Equivalence&e); }; /*:2*/ @@ -21,7 +21,7 @@ static bool selectEquiv(const Equivalence&e); class FNormalMoments:public TensorContainer<FRSingleTensor> { public: -FNormalMoments(const UNormalMoments&moms); + FNormalMoments(const UNormalMoments&moms); }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/permutation.h b/mex/sources/korderpert/Dyn_pp/tl/cc/permutation.h index 0ce3345bee364b646739eeb4f610fc6ef97b0463..3388aa2a2e7e1b9aa79ada975cf33ad14a56f27a 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/permutation.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/permutation.h @@ -12,40 +12,40 @@ class Permutation{ protected: -IntSequence permap; + IntSequence permap; public: -Permutation(int len) -:permap(len){for(int i= 0;i<len;i++)permap[i]= i;} -Permutation(const Equivalence&e) -:permap(e.getN()){e.trace(permap);} -Permutation(const Equivalence&e,const Permutation&per) -:permap(e.getN()){e.trace(permap,per);} -Permutation(const IntSequence&s) -:permap(s.size()){computeSortingMap(s);}; -Permutation(const Permutation&p) -:permap(p.permap){} -Permutation(const Permutation&p1,const Permutation&p2) -:permap(p2.permap){p1.apply(permap);} -Permutation(const Permutation&p,int i) -:permap(p.size(),p.permap,i){} -const Permutation&operator= (const Permutation&p) -{permap= p.permap;return*this;} -bool operator==(const Permutation&p) -{return permap==p.permap;} -int size()const -{return permap.size();} -void print()const -{permap.print();} -void apply(const IntSequence&src,IntSequence&tar)const; -void apply(IntSequence&tar)const; -void inverse(); -int tailIdentity()const; -const IntSequence&getMap()const -{return permap;} -IntSequence&getMap() -{return permap;} + Permutation(int len) + :permap(len){for(int i= 0;i<len;i++)permap[i]= i;} + Permutation(const Equivalence&e) + :permap(e.getN()){e.trace(permap);} + Permutation(const Equivalence&e,const Permutation&per) + :permap(e.getN()){e.trace(permap,per);} + Permutation(const IntSequence&s) + :permap(s.size()){computeSortingMap(s);}; + Permutation(const Permutation&p) + :permap(p.permap){} + Permutation(const Permutation&p1,const Permutation&p2) + :permap(p2.permap){p1.apply(permap);} + Permutation(const Permutation&p,int i) + :permap(p.size(),p.permap,i){} + const Permutation&operator= (const Permutation&p) + {permap= p.permap;return*this;} + bool operator==(const Permutation&p) + {return permap==p.permap;} + int size()const + {return permap.size();} + void print()const + {permap.print();} + void apply(const IntSequence&src,IntSequence&tar)const; + void apply(IntSequence&tar)const; + void inverse(); + int tailIdentity()const; + const IntSequence&getMap()const + {return permap;} + IntSequence&getMap() + {return permap;} protected: -void computeSortingMap(const IntSequence&s); + void computeSortingMap(const IntSequence&s); }; @@ -54,18 +54,18 @@ void computeSortingMap(const IntSequence&s); /*3:*/ class PermutationSet{ -int order; -int size; -const Permutation**const pers; + int order; + int size; + const Permutation**const pers; public: -PermutationSet(); -PermutationSet(const PermutationSet&ps,int n); -~PermutationSet(); -int getNum()const -{return size;} -const Permutation&get(int i)const -{return*(pers[i]);} -vector<const Permutation*> getPreserving(const IntSequence&s)const; + PermutationSet(); + PermutationSet(const PermutationSet&ps,int n); + ~PermutationSet(); + int getNum()const + {return size;} + const Permutation&get(int i)const + {return*(pers[i]);} + vector<const Permutation*> getPreserving(const IntSequence&s)const; }; @@ -74,12 +74,12 @@ vector<const Permutation*> getPreserving(const IntSequence&s)const; /*4:*/ class PermutationBundle{ -vector<PermutationSet*> bundle; + vector<PermutationSet*> bundle; public: -PermutationBundle(int nmax); -~PermutationBundle(); -const PermutationSet&get(int n)const; -void generateUpTo(int nmax); + PermutationBundle(int nmax); + ~PermutationBundle(); + const PermutationSet&get(int n)const; + void generateUpTo(int nmax); }; /*:4*/ diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/ps_tensor.h b/mex/sources/korderpert/Dyn_pp/tl/cc/ps_tensor.h index b2df33b392453cbf27e69b60ed6908e44b2d5576..081d6e19cbbf32e4c34bf04c72d16082150751d4 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/ps_tensor.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/ps_tensor.h @@ -15,8 +15,8 @@ class SortIntSequence:public IntSequence{ public: -SortIntSequence(const IntSequence&s) -:IntSequence(s){sort();} + SortIntSequence(const IntSequence&s) + :IntSequence(s){sort();} }; @@ -26,31 +26,31 @@ SortIntSequence(const IntSequence&s) class PerTensorDimens:public TensorDimens{ protected: -Permutation per; + Permutation per; public: -PerTensorDimens(const Symmetry&s,const IntSequence&nvars, -const Equivalence&e) -:TensorDimens(s,nvars),per(e) -{per.apply(nvmax);} -PerTensorDimens(const TensorDimens&td,const Equivalence&e) -:TensorDimens(td),per(e) -{per.apply(nvmax);} -PerTensorDimens(const TensorDimens&td,const Permutation&p) -:TensorDimens(td),per(p) -{per.apply(nvmax);} -PerTensorDimens(const IntSequence&ss,const IntSequence&coor) -:TensorDimens(ss,SortIntSequence(coor)),per(coor) -{per.apply(nvmax);} -PerTensorDimens(const PerTensorDimens&td) -:TensorDimens(td),per(td.per){} -const PerTensorDimens&operator= (const PerTensorDimens&td) -{TensorDimens::operator= (td);per= td.per;return*this;} -bool operator==(const PerTensorDimens&td) -{return TensorDimens::operator==(td)&&per==td.per;} -int tailIdentity()const -{return per.tailIdentity();} -const Permutation&getPer()const -{return per;} + PerTensorDimens(const Symmetry&s,const IntSequence&nvars, + const Equivalence&e) + :TensorDimens(s,nvars),per(e) + {per.apply(nvmax);} + PerTensorDimens(const TensorDimens&td,const Equivalence&e) + :TensorDimens(td),per(e) + {per.apply(nvmax);} + PerTensorDimens(const TensorDimens&td,const Permutation&p) + :TensorDimens(td),per(p) + {per.apply(nvmax);} + PerTensorDimens(const IntSequence&ss,const IntSequence&coor) + :TensorDimens(ss,SortIntSequence(coor)),per(coor) + {per.apply(nvmax);} + PerTensorDimens(const PerTensorDimens&td) + :TensorDimens(td),per(td.per){} + const PerTensorDimens&operator= (const PerTensorDimens&td) + {TensorDimens::operator= (td);per= td.per;return*this;} + bool operator==(const PerTensorDimens&td) + {return TensorDimens::operator==(td)&&per==td.per;} + int tailIdentity()const + {return per.tailIdentity();} + const Permutation&getPer()const + {return per;} }; /*:3*/ @@ -58,54 +58,54 @@ const Permutation&getPer()const /*4:*/ class UPSTensor:public UTensor{ -const PerTensorDimens tdims; + const PerTensorDimens tdims; public: -/*5:*/ - -UPSTensor(const TensorDimens&td,const Equivalence&e, -const ConstTwoDMatrix&a,const KronProdAll&kp) -:UTensor(along_col,PerTensorDimens(td,e).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,e) -{kp.mult(a,*this);} -UPSTensor(const TensorDimens&td,const Equivalence&e, -const ConstTwoDMatrix&a,const KronProdAllOptim&kp) -:UTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,kp.getPer())) -{kp.mult(a,*this);} -UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, -const ConstTwoDMatrix&a,const KronProdAll&kp) -:UTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,p)) -{kp.mult(a,*this);} -UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, -const ConstTwoDMatrix&a,const KronProdAllOptim&kp) -:UTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,Permutation(p,kp.getPer()))) -{kp.mult(a,*this);} - -/*:5*/ -; -UPSTensor(const FSSparseTensor&t,const IntSequence&ss, -const IntSequence&coor,const PerTensorDimens&ptd); -UPSTensor(const UPSTensor&ut) -:UTensor(ut),tdims(ut.tdims){} - -void increment(IntSequence&v)const; -void decrement(IntSequence&v)const; -FTensor&fold()const; - -int getOffset(const IntSequence&v)const; -void addTo(FGSTensor&out)const; -void addTo(UGSTensor&out)const; - -enum fill_method{first,second}; -static fill_method decideFillMethod(const FSSparseTensor&t); + /*5:*/ + + UPSTensor(const TensorDimens&td,const Equivalence&e, + const ConstTwoDMatrix&a,const KronProdAll&kp) + :UTensor(along_col,PerTensorDimens(td,e).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,e) + {kp.mult(a,*this);} + UPSTensor(const TensorDimens&td,const Equivalence&e, + const ConstTwoDMatrix&a,const KronProdAllOptim&kp) + :UTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,kp.getPer())) + {kp.mult(a,*this);} + UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, + const ConstTwoDMatrix&a,const KronProdAll&kp) + :UTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,p)) + {kp.mult(a,*this);} + UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, + const ConstTwoDMatrix&a,const KronProdAllOptim&kp) + :UTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,Permutation(p,kp.getPer()))) + {kp.mult(a,*this);} + + /*:5*/ + ; + UPSTensor(const FSSparseTensor&t,const IntSequence&ss, + const IntSequence&coor,const PerTensorDimens&ptd); + UPSTensor(const UPSTensor&ut) + :UTensor(ut),tdims(ut.tdims){} + + void increment(IntSequence&v)const; + void decrement(IntSequence&v)const; + FTensor&fold()const; + + int getOffset(const IntSequence&v)const; + void addTo(FGSTensor&out)const; + void addTo(UGSTensor&out)const; + + enum fill_method{first,second}; + static fill_method decideFillMethod(const FSSparseTensor&t); private: -int tailIdentitySize()const; -void fillFromSparseOne(const FSSparseTensor&t,const IntSequence&ss, -const IntSequence&coor); -void fillFromSparseTwo(const FSSparseTensor&t,const IntSequence&ss, -const IntSequence&coor); + int tailIdentitySize()const; + void fillFromSparseOne(const FSSparseTensor&t,const IntSequence&ss, + const IntSequence&coor); + void fillFromSparseTwo(const FSSparseTensor&t,const IntSequence&ss, + const IntSequence&coor); }; /*:4*/ @@ -113,30 +113,30 @@ const IntSequence&coor); /*6:*/ class PerTensorDimens2:public PerTensorDimens{ -InducedSymmetries syms; -IntSequence ds; + InducedSymmetries syms; + IntSequence ds; public: -PerTensorDimens2(const TensorDimens&td,const Equivalence&e, -const Permutation&p) -:PerTensorDimens(td,Permutation(e,p)), -syms(e,p,td.getSym()), -ds(syms.size()) -{setDimensionSizes();} -PerTensorDimens2(const TensorDimens&td,const Equivalence&e) -:PerTensorDimens(td,e), -syms(e,td.getSym()), -ds(syms.size()) -{setDimensionSizes();} -int numSyms()const -{return(int)syms.size();} -const Symmetry&getSym(int i)const -{return syms[i];} -int calcMaxOffset()const -{return ds.mult();} -int calcOffset(const IntSequence&coor)const; -void print()const; + PerTensorDimens2(const TensorDimens&td,const Equivalence&e, + const Permutation&p) + :PerTensorDimens(td,Permutation(e,p)), + syms(e,p,td.getSym()), + ds(syms.size()) + {setDimensionSizes();} + PerTensorDimens2(const TensorDimens&td,const Equivalence&e) + :PerTensorDimens(td,e), + syms(e,td.getSym()), + ds(syms.size()) + {setDimensionSizes();} + int numSyms()const + {return(int)syms.size();} + const Symmetry&getSym(int i)const + {return syms[i];} + int calcMaxOffset()const + {return ds.mult();} + int calcOffset(const IntSequence&coor)const; + void print()const; protected: -void setDimensionSizes(); + void setDimensionSizes(); }; /*:6*/ @@ -146,46 +146,46 @@ void setDimensionSizes(); template<typename _Ttype> class StackProduct; class FPSTensor:public FTensor{ -const PerTensorDimens2 tdims; + const PerTensorDimens2 tdims; public: -/*8:*/ - -FPSTensor(const TensorDimens&td,const Equivalence&e, -const ConstTwoDMatrix&a,const KronProdAll&kp) -:FTensor(along_col,PerTensorDimens(td,e).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,e) -{kp.mult(a,*this);} -FPSTensor(const TensorDimens&td,const Equivalence&e, -const ConstTwoDMatrix&a,const KronProdAllOptim&kp) -:FTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,kp.getPer()) -{kp.mult(a,*this);} -FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, -const ConstTwoDMatrix&a,const KronProdAll&kp) -:FTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,p) -{kp.mult(a,*this);} -FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, -const ConstTwoDMatrix&a,const KronProdAllOptim&kp) -:FTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(), -a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,Permutation(p,kp.getPer())) -{kp.mult(a,*this);} - -FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, -const GSSparseTensor&t,const KronProdAll&kp); - -FPSTensor(const FPSTensor&ft) -:FTensor(ft),tdims(ft.tdims){} - -/*:8*/ -; - -void increment(IntSequence&v)const; -void decrement(IntSequence&v)const; -UTensor&unfold()const; - -int getOffset(const IntSequence&v)const; -void addTo(FGSTensor&out)const; + /*8:*/ + + FPSTensor(const TensorDimens&td,const Equivalence&e, + const ConstTwoDMatrix&a,const KronProdAll&kp) + :FTensor(along_col,PerTensorDimens(td,e).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,e) + {kp.mult(a,*this);} + FPSTensor(const TensorDimens&td,const Equivalence&e, + const ConstTwoDMatrix&a,const KronProdAllOptim&kp) + :FTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,kp.getPer()) + {kp.mult(a,*this);} + FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, + const ConstTwoDMatrix&a,const KronProdAll&kp) + :FTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,p) + {kp.mult(a,*this);} + FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, + const ConstTwoDMatrix&a,const KronProdAllOptim&kp) + :FTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(), + a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,Permutation(p,kp.getPer())) + {kp.mult(a,*this);} + + FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p, + const GSSparseTensor&t,const KronProdAll&kp); + + FPSTensor(const FPSTensor&ft) + :FTensor(ft),tdims(ft.tdims){} + + /*:8*/ + ; + + void increment(IntSequence&v)const; + void decrement(IntSequence&v)const; + UTensor&unfold()const; + + int getOffset(const IntSequence&v)const; + void addTo(FGSTensor&out)const; }; /*:7*/ diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod.h b/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod.h index f3a51927366079478f194b09ffad0960d57c265f..795e30c115b0c6700e7183bbfc1ac8f20927067f 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod.h @@ -16,10 +16,10 @@ using namespace std; class USubTensor:public URTensor{ public: -USubTensor(const TensorDimens&bdims,const TensorDimens&hdims, -const FGSContainer&cont,const vector<IntSequence> &lst); -void addKronColumn(int i,const vector<const FGSTensor*> &ts, -const IntSequence&pindex); + USubTensor(const TensorDimens&bdims,const TensorDimens&hdims, + const FGSContainer&cont,const vector<IntSequence> &lst); + void addKronColumn(int i,const vector<const FGSTensor*> &ts, + const IntSequence&pindex); }; /*:2*/ diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod2.h b/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod2.h index ee0b80657973c4d9a15dcbd57672a216576ff6ee..266056931784191a188d5f78a04f84502efce920 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod2.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/pyramid_prod2.h @@ -15,20 +15,20 @@ class IrregTensor; class IrregTensorHeader{ -friend class IrregTensor; -int nv; -IntSequence unit_flag; -Vector**const cols; -IntSequence end_seq; + friend class IrregTensor; + int nv; + IntSequence unit_flag; + Vector**const cols; + IntSequence end_seq; public: -IrregTensorHeader(const StackProduct<FGSTensor> &sp,const IntSequence&c); -~IrregTensorHeader(); -int dimen()const -{return unit_flag.size();} -void increment(IntSequence&v)const; -int calcMaxOffset()const; + IrregTensorHeader(const StackProduct<FGSTensor> &sp,const IntSequence&c); + ~IrregTensorHeader(); + int dimen()const + {return unit_flag.size();} + void increment(IntSequence&v)const; + int calcMaxOffset()const; private: -IrregTensorHeader(const IrregTensorHeader&); + IrregTensorHeader(const IrregTensorHeader&); }; @@ -37,16 +37,16 @@ IrregTensorHeader(const IrregTensorHeader&); /*3:*/ class IrregTensor:public Tensor{ -const IrregTensorHeader&header; + const IrregTensorHeader&header; public: -IrregTensor(const IrregTensorHeader&h); -void addTo(FRSingleTensor&out)const; -void increment(IntSequence&v)const -{header.increment(v);} -void decrement(IntSequence&v)const -{TL_RAISE("Not implemented error in IrregTensor::decrement");} -int getOffset(const IntSequence&v)const -{TL_RAISE("Not implemented error in IrregTensor::getOffset");return 0;} + IrregTensor(const IrregTensorHeader&h); + void addTo(FRSingleTensor&out)const; + void increment(IntSequence&v)const + {header.increment(v);} + void decrement(IntSequence&v)const + {TL_RAISE("Not implemented error in IrregTensor::decrement");} + int getOffset(const IntSequence&v)const + {TL_RAISE("Not implemented error in IrregTensor::getOffset");return 0;} }; /*:3*/ diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/rfs_tensor.h b/mex/sources/korderpert/Dyn_pp/tl/cc/rfs_tensor.h index 10fbf247917176750187caf345849ddc63c0b109..00b3072bbd41b7c52928adf2618bb6505860a286 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/rfs_tensor.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/rfs_tensor.h @@ -11,30 +11,30 @@ class FRTensor; class URTensor:public UTensor{ -int nv; + int nv; public: -/*3:*/ - -URTensor(int c,int nvar,int d) -:UTensor(along_row,IntSequence(d,nvar), -UFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){} -URTensor(const URTensor&ut) -:UTensor(ut),nv(ut.nv){} -URTensor(const FRTensor&ft); - -/*:3*/ -; -virtual~URTensor(){} - -void increment(IntSequence&v)const; -void decrement(IntSequence&v)const; -FTensor&fold()const; - -int getOffset(const IntSequence&v)const; -int nvar()const -{return nv;} -Symmetry getSym()const -{return Symmetry(dimen());} + /*3:*/ + + URTensor(int c,int nvar,int d) + :UTensor(along_row,IntSequence(d,nvar), + UFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){} + URTensor(const URTensor&ut) + :UTensor(ut),nv(ut.nv){} + URTensor(const FRTensor&ft); + + /*:3*/ + ; + virtual~URTensor(){} + + void increment(IntSequence&v)const; + void decrement(IntSequence&v)const; + FTensor&fold()const; + + int getOffset(const IntSequence&v)const; + int nvar()const + {return nv;} + Symmetry getSym()const + {return Symmetry(dimen());} }; /*:2*/ @@ -42,31 +42,31 @@ Symmetry getSym()const /*4:*/ class FRTensor:public FTensor{ -int nv; + int nv; public: -/*5:*/ - -FRTensor(int c,int nvar,int d) -:FTensor(along_row,IntSequence(d,nvar), -FFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){} -FRTensor(const FRTensor&ft) -:FTensor(ft),nv(ft.nv){} -FRTensor(const URTensor&ut); - -/*:5*/ -; -virtual~FRTensor(){} - -void increment(IntSequence&v)const; -void decrement(IntSequence&v)const; -UTensor&unfold()const; - -int nvar()const -{return nv;} -int getOffset(const IntSequence&v)const -{return FTensor::getOffset(v,nv);} -Symmetry getSym()const -{return Symmetry(dimen());} + /*5:*/ + + FRTensor(int c,int nvar,int d) + :FTensor(along_row,IntSequence(d,nvar), + FFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){} + FRTensor(const FRTensor&ft) + :FTensor(ft),nv(ft.nv){} + FRTensor(const URTensor&ut); + + /*:5*/ + ; + virtual~FRTensor(){} + + void increment(IntSequence&v)const; + void decrement(IntSequence&v)const; + UTensor&unfold()const; + + int nvar()const + {return nv;} + int getOffset(const IntSequence&v)const + {return FTensor::getOffset(v,nv);} + Symmetry getSym()const + {return Symmetry(dimen());} }; /*:4*/ @@ -75,14 +75,14 @@ Symmetry getSym()const class URSingleTensor:public URTensor{ public: -URSingleTensor(int nvar,int d) -:URTensor(1,nvar,d){} -URSingleTensor(const vector<ConstVector> &cols); -URSingleTensor(const ConstVector&v,int d); -URSingleTensor(const URSingleTensor&ut) -:URTensor(ut){} -virtual~URSingleTensor(){} -FTensor&fold()const; + URSingleTensor(int nvar,int d) + :URTensor(1,nvar,d){} + URSingleTensor(const vector<ConstVector> &cols); + URSingleTensor(const ConstVector&v,int d); + URSingleTensor(const URSingleTensor&ut) + :URTensor(ut){} + virtual~URSingleTensor(){} + FTensor&fold()const; }; /*:6*/ @@ -91,12 +91,12 @@ FTensor&fold()const; class FRSingleTensor:public FRTensor{ public: -FRSingleTensor(int nvar,int d) -:FRTensor(1,nvar,d){} -FRSingleTensor(const URSingleTensor&ut); -FRSingleTensor(const FRSingleTensor&ft) -:FRTensor(ft){} -virtual~FRSingleTensor(){} + FRSingleTensor(int nvar,int d) + :FRTensor(1,nvar,d){} + FRSingleTensor(const URSingleTensor&ut); + FRSingleTensor(const FRSingleTensor&ft) + :FRTensor(ft){} + virtual~FRSingleTensor(){} }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/sparse_tensor.h b/mex/sources/korderpert/Dyn_pp/tl/cc/sparse_tensor.h index 66976d7ce26a27b9c3c828f81ee3cf9fb2915021..110ee3950bfaebbb65d361a4868ee772bc804893 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/sparse_tensor.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/sparse_tensor.h @@ -15,8 +15,8 @@ using namespace std; /*2:*/ struct ltseq{ -bool operator()(const IntSequence&s1,const IntSequence&s2)const -{return s1<s2;} + bool operator()(const IntSequence&s1,const IntSequence&s2)const + {return s1<s2;} }; /*:2*/ @@ -25,46 +25,46 @@ bool operator()(const IntSequence&s1,const IntSequence&s2)const class SparseTensor{ public: -typedef pair<int,double> Item; -typedef multimap<IntSequence,Item,ltseq> Map; -typedef Map::const_iterator const_iterator; + typedef pair<int,double> Item; + typedef multimap<IntSequence,Item,ltseq> Map; + typedef Map::const_iterator const_iterator; protected: -typedef Map::iterator iterator; - -Map m; -const int dim; -const int nr; -const int nc; -int first_nz_row; -int last_nz_row; + typedef Map::iterator iterator; + + Map m; + const int dim; + const int nr; + const int nc; + int first_nz_row; + int last_nz_row; public: -SparseTensor(int d,int nnr,int nnc) -:dim(d),nr(nnr),nc(nnc),first_nz_row(nr),last_nz_row(-1){} -SparseTensor(const SparseTensor&t) -:m(t.m),dim(t.dim),nr(t.nr),nc(t.nc){} -virtual~SparseTensor(){} -void insert(const IntSequence&s,int r,double c); -const Map&getMap()const -{return m;} -int dimen()const -{return dim;} -int nrows()const -{return nr;} -int ncols()const -{return nc;} -double getFillFactor()const -{return((double)m.size())/(nrows()*ncols());} -double getFoldIndexFillFactor()const; -double getUnfoldIndexFillFactor()const; -int getNumNonZero()const -{return m.size();} -int getFirstNonZeroRow()const -{return first_nz_row;} -int getLastNonZeroRow()const -{return last_nz_row;} -virtual const Symmetry&getSym()const= 0; -void print()const; -bool isFinite()const; + SparseTensor(int d,int nnr,int nnc) + :dim(d),nr(nnr),nc(nnc),first_nz_row(nr),last_nz_row(-1){} + SparseTensor(const SparseTensor&t) + :m(t.m),dim(t.dim),nr(t.nr),nc(t.nc){} + virtual~SparseTensor(){} + void insert(const IntSequence&s,int r,double c); + const Map&getMap()const + {return m;} + int dimen()const + {return dim;} + int nrows()const + {return nr;} + int ncols()const + {return nc;} + double getFillFactor()const + {return((double)m.size())/(nrows()*ncols());} + double getFoldIndexFillFactor()const; + double getUnfoldIndexFillFactor()const; + int getNumNonZero()const + {return m.size();} + int getFirstNonZeroRow()const + {return first_nz_row;} + int getLastNonZeroRow()const + {return last_nz_row;} + virtual const Symmetry&getSym()const= 0; + void print()const; + bool isFinite()const; } /*:3*/ @@ -73,20 +73,20 @@ bool isFinite()const; class FSSparseTensor:public SparseTensor{ public: -typedef SparseTensor::const_iterator const_iterator; + typedef SparseTensor::const_iterator const_iterator; private: -const int nv; -const Symmetry sym; + const int nv; + const Symmetry sym; public: -FSSparseTensor(int d,int nvar,int r); -FSSparseTensor(const FSSparseTensor&t); -void insert(const IntSequence&s,int r,double c); -void multColumnAndAdd(const Tensor&t,Vector&v)const; -const Symmetry&getSym()const -{return sym;} -int nvar()const -{return nv;} -void print()const; + FSSparseTensor(int d,int nvar,int r); + FSSparseTensor(const FSSparseTensor&t); + void insert(const IntSequence&s,int r,double c); + void multColumnAndAdd(const Tensor&t,Vector&v)const; + const Symmetry&getSym()const + {return sym;} + int nvar()const + {return nv;} + void print()const; }; @@ -96,21 +96,21 @@ void print()const; class GSSparseTensor:public SparseTensor{ public: -typedef SparseTensor::const_iterator const_iterator; + typedef SparseTensor::const_iterator const_iterator; private: -const TensorDimens tdims; + const TensorDimens tdims; public: -GSSparseTensor(const FSSparseTensor&t,const IntSequence&ss, -const IntSequence&coor,const TensorDimens&td); -GSSparseTensor(const GSSparseTensor&t) -:SparseTensor(t),tdims(t.tdims){} -void insert(const IntSequence&s,int r,double c); -const Symmetry&getSym()const -{return tdims.getSym();} -const TensorDimens&getDims()const -{return tdims;} -void print()const; - + GSSparseTensor(const FSSparseTensor&t,const IntSequence&ss, + const IntSequence&coor,const TensorDimens&td); + GSSparseTensor(const GSSparseTensor&t) + :SparseTensor(t),tdims(t.tdims){} + void insert(const IntSequence&s,int r,double c); + const Symmetry&getSym()const + {return tdims.getSym();} + const TensorDimens&getDims()const + {return tdims;} + void print()const; + }; /*:5*/ diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/stack_container.h b/mex/sources/korderpert/Dyn_pp/tl/cc/stack_container.h index 55fb960394b65340d7ee082f47127ff23cbf00d8..e9bc5ad67a817bf4cf49752af018b9aaa1357b17 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/stack_container.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/stack_container.h @@ -16,32 +16,32 @@ template<class _Ttype> class StackContainerInterface{ public: -typedef TensorContainer<_Ttype> _Ctype; -typedef enum{matrix,unit,zero}itype; + typedef TensorContainer<_Ttype> _Ctype; + typedef enum{matrix,unit,zero}itype; protected: -const EquivalenceBundle&ebundle; + const EquivalenceBundle&ebundle; public: -StackContainerInterface() -:ebundle(*(tls.ebundle)){} -virtual~StackContainerInterface(){} -virtual const IntSequence&getStackSizes()const= 0; -virtual IntSequence&getStackSizes()= 0; -virtual const IntSequence&getStackOffsets()const= 0; -virtual IntSequence&getStackOffsets()= 0; -virtual int numConts()const= 0; -virtual const _Ctype*getCont(int i)const= 0; -virtual itype getType(int i,const Symmetry&s)const= 0; -virtual int numStacks()const= 0; -virtual bool isZero(int i,const Symmetry&s)const= 0; -virtual const _Ttype*getMatrix(int i,const Symmetry&s)const= 0; -virtual int getLengthOfMatrixStacks(const Symmetry&s)const= 0; -virtual int getUnitPos(const Symmetry&s)const= 0; -virtual Vector*createPackedColumn(const Symmetry&s, -const IntSequence&coor, -int&iu)const= 0; -int getAllSize()const -{return getStackOffsets()[numStacks()-1] -+getStackSizes()[numStacks()-1];} + StackContainerInterface() + :ebundle(*(tls.ebundle)){} + virtual~StackContainerInterface(){} + virtual const IntSequence&getStackSizes()const= 0; + virtual IntSequence&getStackSizes()= 0; + virtual const IntSequence&getStackOffsets()const= 0; + virtual IntSequence&getStackOffsets()= 0; + virtual int numConts()const= 0; + virtual const _Ctype*getCont(int i)const= 0; + virtual itype getType(int i,const Symmetry&s)const= 0; + virtual int numStacks()const= 0; + virtual bool isZero(int i,const Symmetry&s)const= 0; + virtual const _Ttype*getMatrix(int i,const Symmetry&s)const= 0; + virtual int getLengthOfMatrixStacks(const Symmetry&s)const= 0; + virtual int getUnitPos(const Symmetry&s)const= 0; + virtual Vector*createPackedColumn(const Symmetry&s, + const IntSequence&coor, + int&iu)const= 0; + int getAllSize()const + {return getStackOffsets()[numStacks()-1] + +getStackSizes()[numStacks()-1];} }; /*:2*/ @@ -51,131 +51,131 @@ int getAllSize()const template<class _Ttype> class StackContainer:virtual public StackContainerInterface<_Ttype> { public: -typedef StackContainerInterface<_Ttype> _Stype; -typedef typename StackContainerInterface<_Ttype> ::_Ctype _Ctype; -typedef typename StackContainerInterface<_Ttype> ::itype itype; + typedef StackContainerInterface<_Ttype> _Stype; + typedef typename StackContainerInterface<_Ttype> ::_Ctype _Ctype; + typedef typename StackContainerInterface<_Ttype> ::itype itype; protected: -int num_conts; -IntSequence stack_sizes; -IntSequence stack_offsets; -const _Ctype**const conts; + int num_conts; + IntSequence stack_sizes; + IntSequence stack_offsets; + const _Ctype**const conts; public: -StackContainer(int ns,int nc) -:num_conts(nc),stack_sizes(ns,0),stack_offsets(ns,0), -conts(new const _Ctype*[nc]){} -virtual~StackContainer(){delete[]conts;} -const IntSequence&getStackSizes()const -{return stack_sizes;} -IntSequence&getStackSizes() -{return stack_sizes;} -const IntSequence&getStackOffsets()const -{return stack_offsets;} -IntSequence&getStackOffsets() -{return stack_offsets;} -int numConts()const -{return num_conts;} -const _Ctype*getCont(int i)const -{return conts[i];} -virtual itype getType(int i,const Symmetry&s)const= 0; -int numStacks()const -{return stack_sizes.size();} -/*4:*/ - -bool isZero(int i,const Symmetry&s)const -{ -TL_RAISE_IF(i<0||i>=numStacks(), -"Wrong index to stack in StackContainer::isZero."); -return(getType(i,s)==_Stype::zero|| -(getType(i,s)==_Stype::matrix&&!conts[i]->check(s))); -} - -/*:4*/ -; -/*5:*/ - -const _Ttype*getMatrix(int i,const Symmetry&s)const -{ -TL_RAISE_IF(isZero(i,s)||getType(i,s)==_Stype::unit, -"Matrix is not returned in StackContainer::getMatrix"); -return conts[i]->get(s); -} - -/*:5*/ -; -/*6:*/ - -int getLengthOfMatrixStacks(const Symmetry&s)const -{ -int res= 0; -int i= 0; -while(i<numStacks()&&getType(i,s)==_Stype::matrix) -res+= stack_sizes[i++]; -return res; -} - - -/*:6*/ -; -/*7:*/ - -int getUnitPos(const Symmetry&s)const -{ -if(s.dimen()!=1) -return-1; -int i= numStacks()-1; -while(i>=0&&getType(i,s)!=_Stype::unit) -i--; -return i; -} - - -/*:7*/ -; -/*8:*/ - -Vector*createPackedColumn(const Symmetry&s, -const IntSequence&coor,int&iu)const -{ -TL_RAISE_IF(s.dimen()!=coor.size(), -"Incompatible coordinates for symmetry in StackContainer::createPackedColumn"); - -int len= getLengthOfMatrixStacks(s); -iu= -1; -int i= 0; -if(-1!=(i= getUnitPos(s))){ -iu= stack_offsets[i]+coor[0]; -len++; -} - -Vector*res= new Vector(len); -i= 0; -while(i<numStacks()&&getType(i,s)==_Stype::matrix){ -const _Ttype*t= getMatrix(i,s); -Tensor::index ind(t,coor); -Vector subres(*res,stack_offsets[i],stack_sizes[i]); -subres= ConstVector(ConstGeneralMatrix(*t),*ind); -i++; -} -if(iu!=-1) -(*res)[len-1]= 1; - -return res; -} - -/*:8*/ -; + StackContainer(int ns,int nc) + :num_conts(nc),stack_sizes(ns,0),stack_offsets(ns,0), + conts(new const _Ctype*[nc]){} + virtual~StackContainer(){delete[]conts;} + const IntSequence&getStackSizes()const + {return stack_sizes;} + IntSequence&getStackSizes() + {return stack_sizes;} + const IntSequence&getStackOffsets()const + {return stack_offsets;} + IntSequence&getStackOffsets() + {return stack_offsets;} + int numConts()const + {return num_conts;} + const _Ctype*getCont(int i)const + {return conts[i];} + virtual itype getType(int i,const Symmetry&s)const= 0; + int numStacks()const + {return stack_sizes.size();} + /*4:*/ + + bool isZero(int i,const Symmetry&s)const + { + TL_RAISE_IF(i<0||i>=numStacks(), + "Wrong index to stack in StackContainer::isZero."); + return(getType(i,s)==_Stype::zero|| + (getType(i,s)==_Stype::matrix&&!conts[i]->check(s))); + } + + /*:4*/ + ; + /*5:*/ + + const _Ttype*getMatrix(int i,const Symmetry&s)const + { + TL_RAISE_IF(isZero(i,s)||getType(i,s)==_Stype::unit, + "Matrix is not returned in StackContainer::getMatrix"); + return conts[i]->get(s); + } + + /*:5*/ + ; + /*6:*/ + + int getLengthOfMatrixStacks(const Symmetry&s)const + { + int res= 0; + int i= 0; + while(i<numStacks()&&getType(i,s)==_Stype::matrix) + res+= stack_sizes[i++]; + return res; + } + + + /*:6*/ + ; + /*7:*/ + + int getUnitPos(const Symmetry&s)const + { + if(s.dimen()!=1) + return-1; + int i= numStacks()-1; + while(i>=0&&getType(i,s)!=_Stype::unit) + i--; + return i; + } + + + /*:7*/ + ; + /*8:*/ + + Vector*createPackedColumn(const Symmetry&s, + const IntSequence&coor,int&iu)const + { + TL_RAISE_IF(s.dimen()!=coor.size(), + "Incompatible coordinates for symmetry in StackContainer::createPackedColumn"); + + int len= getLengthOfMatrixStacks(s); + iu= -1; + int i= 0; + if(-1!=(i= getUnitPos(s))){ + iu= stack_offsets[i]+coor[0]; + len++; + } + + Vector*res= new Vector(len); + i= 0; + while(i<numStacks()&&getType(i,s)==_Stype::matrix){ + const _Ttype*t= getMatrix(i,s); + Tensor::index ind(t,coor); + Vector subres(*res,stack_offsets[i],stack_sizes[i]); + subres= ConstVector(ConstGeneralMatrix(*t),*ind); + i++; + } + if(iu!=-1) + (*res)[len-1]= 1; + + return res; + } + + /*:8*/ + ; protected: -/*9:*/ - -void calculateOffsets() -{ -stack_offsets[0]= 0; -for(int i= 1;i<stack_offsets.size();i++) -stack_offsets[i]= stack_offsets[i-1]+stack_sizes[i-1]; -} - -/*:9*/ -; + /*9:*/ + + void calculateOffsets() + { + stack_offsets[0]= 0; + for(int i= 1;i<stack_offsets.size();i++) + stack_offsets[i]= stack_offsets[i-1]+stack_sizes[i-1]; + } + + /*:9*/ + ; }; /*:3*/ @@ -187,26 +187,26 @@ class WorkerFoldMAASparse1; class WorkerFoldMAASparse2; class WorkerFoldMAASparse4; class FoldedStackContainer:virtual public StackContainerInterface<FGSTensor> { -friend class WorkerFoldMAADense; -friend class WorkerFoldMAASparse1; -friend class WorkerFoldMAASparse2; -friend class WorkerFoldMAASparse4; + friend class WorkerFoldMAADense; + friend class WorkerFoldMAASparse1; + friend class WorkerFoldMAASparse2; + friend class WorkerFoldMAASparse4; public: -static double fill_threshold; -void multAndAdd(int dim,const TensorContainer<FSSparseTensor> &c, -FGSTensor&out)const -{if(c.check(Symmetry(dim)))multAndAdd(*(c.get(Symmetry(dim))),out);} -void multAndAdd(const FSSparseTensor&t,FGSTensor&out)const; -void multAndAdd(int dim,const FGSContainer&c,FGSTensor&out)const; + static double fill_threshold; + void multAndAdd(int dim,const TensorContainer<FSSparseTensor> &c, + FGSTensor&out)const + {if(c.check(Symmetry(dim)))multAndAdd(*(c.get(Symmetry(dim))),out);} + void multAndAdd(const FSSparseTensor&t,FGSTensor&out)const; + void multAndAdd(int dim,const FGSContainer&c,FGSTensor&out)const; protected: -void multAndAddSparse1(const FSSparseTensor&t,FGSTensor&out)const; -void multAndAddSparse2(const FSSparseTensor&t,FGSTensor&out)const; -void multAndAddSparse3(const FSSparseTensor&t,FGSTensor&out)const; -void multAndAddSparse4(const FSSparseTensor&t,FGSTensor&out)const; -void multAndAddStacks(const IntSequence&fi,const FGSTensor&g, -FGSTensor&out,const void*ad)const; -void multAndAddStacks(const IntSequence&fi,const GSSparseTensor&g, -FGSTensor&out,const void*ad)const; + void multAndAddSparse1(const FSSparseTensor&t,FGSTensor&out)const; + void multAndAddSparse2(const FSSparseTensor&t,FGSTensor&out)const; + void multAndAddSparse3(const FSSparseTensor&t,FGSTensor&out)const; + void multAndAddSparse4(const FSSparseTensor&t,FGSTensor&out)const; + void multAndAddStacks(const IntSequence&fi,const FGSTensor&g, + FGSTensor&out,const void*ad)const; + void multAndAddStacks(const IntSequence&fi,const GSSparseTensor&g, + FGSTensor&out,const void*ad)const; }; @@ -218,21 +218,21 @@ class WorkerUnfoldMAADense; class WorkerUnfoldMAASparse1; class WorkerUnfoldMAASparse2; class UnfoldedStackContainer:virtual public StackContainerInterface<UGSTensor> { -friend class WorkerUnfoldMAADense; -friend class WorkerUnfoldMAASparse1; -friend class WorkerUnfoldMAASparse2; + friend class WorkerUnfoldMAADense; + friend class WorkerUnfoldMAASparse1; + friend class WorkerUnfoldMAASparse2; public: -static double fill_threshold; -void multAndAdd(int dim,const TensorContainer<FSSparseTensor> &c, -UGSTensor&out)const -{if(c.check(Symmetry(dim)))multAndAdd(*(c.get(Symmetry(dim))),out);} -void multAndAdd(const FSSparseTensor&t,UGSTensor&out)const; -void multAndAdd(int dim,const UGSContainer&c,UGSTensor&out)const; + static double fill_threshold; + void multAndAdd(int dim,const TensorContainer<FSSparseTensor> &c, + UGSTensor&out)const + {if(c.check(Symmetry(dim)))multAndAdd(*(c.get(Symmetry(dim))),out);} + void multAndAdd(const FSSparseTensor&t,UGSTensor&out)const; + void multAndAdd(int dim,const UGSContainer&c,UGSTensor&out)const; protected: -void multAndAddSparse1(const FSSparseTensor&t,UGSTensor&out)const; -void multAndAddSparse2(const FSSparseTensor&t,UGSTensor&out)const; -void multAndAddStacks(const IntSequence&fi,const UGSTensor&g, -UGSTensor&out,const void*ad)const; + void multAndAddSparse1(const FSSparseTensor&t,UGSTensor&out)const; + void multAndAddSparse2(const FSSparseTensor&t,UGSTensor&out)const; + void multAndAddStacks(const IntSequence&fi,const UGSTensor&g, + UGSTensor&out,const void*ad)const; }; /*:11*/ @@ -242,49 +242,49 @@ UGSTensor&out,const void*ad)const; template<class _Ttype> class ZContainer:public StackContainer<_Ttype> { public: -typedef StackContainer<_Ttype> _Tparent; -typedef StackContainerInterface<_Ttype> _Stype; -typedef typename _Tparent::_Ctype _Ctype; -typedef typename _Tparent::itype itype; -ZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, -int ny,int nu) -:_Tparent(4,2) -{ -_Tparent::stack_sizes[0]= ngss;_Tparent::stack_sizes[1]= ng; -_Tparent::stack_sizes[2]= ny;_Tparent::stack_sizes[3]= nu; -_Tparent::conts[0]= gss; -_Tparent::conts[1]= g; -_Tparent::calculateOffsets(); -} - -/*13:*/ - -itype getType(int i,const Symmetry&s)const -{ -if(i==0) -return _Stype::matrix; -if(i==1) -if(s[2]> 0) -return _Stype::zero; -else -return _Stype::matrix; -if(i==2) -if(s==Symmetry(1,0,0,0)) -return _Stype::unit; -else -return _Stype::zero; -if(i==3) -if(s==Symmetry(0,1,0,0)) -return _Stype::unit; -else -return _Stype::zero; - -TL_RAISE("Wrong stack index in ZContainer::getType"); -return _Stype::zero; -} - -/*:13*/ -; + typedef StackContainer<_Ttype> _Tparent; + typedef StackContainerInterface<_Ttype> _Stype; + typedef typename _Tparent::_Ctype _Ctype; + typedef typename _Tparent::itype itype; + ZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, + int ny,int nu) + :_Tparent(4,2) + { + _Tparent::stack_sizes[0]= ngss;_Tparent::stack_sizes[1]= ng; + _Tparent::stack_sizes[2]= ny;_Tparent::stack_sizes[3]= nu; + _Tparent::conts[0]= gss; + _Tparent::conts[1]= g; + _Tparent::calculateOffsets(); + } + + /*13:*/ + + itype getType(int i,const Symmetry&s)const + { + if(i==0) + return _Stype::matrix; + if(i==1) + if(s[2]> 0) + return _Stype::zero; + else + return _Stype::matrix; + if(i==2) + if(s==Symmetry(1,0,0,0)) + return _Stype::unit; + else + return _Stype::zero; + if(i==3) + if(s==Symmetry(0,1,0,0)) + return _Stype::unit; + else + return _Stype::zero; + + TL_RAISE("Wrong stack index in ZContainer::getType"); + return _Stype::zero; + } + + /*:13*/ + ; }; /*:12*/ @@ -294,10 +294,10 @@ return _Stype::zero; class FoldedZContainer:public ZContainer<FGSTensor> , public FoldedStackContainer{ public: -typedef TensorContainer<FGSTensor> _Ctype; -FoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, -int ny,int nu) -:ZContainer<FGSTensor> (gss,ngss,g,ng,ny,nu){} + typedef TensorContainer<FGSTensor> _Ctype; + FoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, + int ny,int nu) + :ZContainer<FGSTensor> (gss,ngss,g,ng,ny,nu){} }; /*:14*/ @@ -307,10 +307,10 @@ int ny,int nu) class UnfoldedZContainer:public ZContainer<UGSTensor> , public UnfoldedStackContainer{ public: -typedef TensorContainer<UGSTensor> _Ctype; -UnfoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, -int ny,int nu) -:ZContainer<UGSTensor> (gss,ngss,g,ng,ny,nu){} + typedef TensorContainer<UGSTensor> _Ctype; + UnfoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, + int ny,int nu) + :ZContainer<UGSTensor> (gss,ngss,g,ng,ny,nu){} }; /*:15*/ @@ -320,48 +320,48 @@ int ny,int nu) template<class _Ttype> class GContainer:public StackContainer<_Ttype> { public: -typedef StackContainer<_Ttype> _Tparent; -typedef StackContainerInterface<_Ttype> _Stype; -typedef typename StackContainer<_Ttype> ::_Ctype _Ctype; -typedef typename StackContainer<_Ttype> ::itype itype; -GContainer(const _Ctype*gs,int ngs,int nu) -:StackContainer<_Ttype> (4,1) -{ -_Tparent::stack_sizes[0]= ngs;_Tparent::stack_sizes[1]= nu; -_Tparent::stack_sizes[2]= nu;_Tparent::stack_sizes[3]= 1; -_Tparent::conts[0]= gs; -_Tparent::calculateOffsets(); -} - -/*17:*/ - -itype getType(int i,const Symmetry&s)const -{ -if(i==0) -if(s[2]> 0||s==Symmetry(0,0,0,1)) -return _Stype::zero; -else -return _Stype::matrix; -if(i==1) -if(s==Symmetry(0,0,1,0)) -return _Stype::unit; -else -return _Stype::zero; -if(i==2) -return _Stype::zero; -if(i==3) -if(s==Symmetry(0,0,0,1)) -return _Stype::unit; -else -return _Stype::zero; - -TL_RAISE("Wrong stack index in GContainer::getType"); -return _Stype::zero; -} - - -/*:17*/ -; + typedef StackContainer<_Ttype> _Tparent; + typedef StackContainerInterface<_Ttype> _Stype; + typedef typename StackContainer<_Ttype> ::_Ctype _Ctype; + typedef typename StackContainer<_Ttype> ::itype itype; + GContainer(const _Ctype*gs,int ngs,int nu) + :StackContainer<_Ttype> (4,1) + { + _Tparent::stack_sizes[0]= ngs;_Tparent::stack_sizes[1]= nu; + _Tparent::stack_sizes[2]= nu;_Tparent::stack_sizes[3]= 1; + _Tparent::conts[0]= gs; + _Tparent::calculateOffsets(); + } + + /*17:*/ + + itype getType(int i,const Symmetry&s)const + { + if(i==0) + if(s[2]> 0||s==Symmetry(0,0,0,1)) + return _Stype::zero; + else + return _Stype::matrix; + if(i==1) + if(s==Symmetry(0,0,1,0)) + return _Stype::unit; + else + return _Stype::zero; + if(i==2) + return _Stype::zero; + if(i==3) + if(s==Symmetry(0,0,0,1)) + return _Stype::unit; + else + return _Stype::zero; + + TL_RAISE("Wrong stack index in GContainer::getType"); + return _Stype::zero; + } + + + /*:17*/ + ; }; /*:16*/ @@ -371,9 +371,9 @@ return _Stype::zero; class FoldedGContainer:public GContainer<FGSTensor> , public FoldedStackContainer{ public: -typedef TensorContainer<FGSTensor> _Ctype; -FoldedGContainer(const _Ctype*gs,int ngs,int nu) -:GContainer<FGSTensor> (gs,ngs,nu){} + typedef TensorContainer<FGSTensor> _Ctype; + FoldedGContainer(const _Ctype*gs,int ngs,int nu) + :GContainer<FGSTensor> (gs,ngs,nu){} }; /*:18*/ @@ -383,9 +383,9 @@ FoldedGContainer(const _Ctype*gs,int ngs,int nu) class UnfoldedGContainer:public GContainer<UGSTensor> , public UnfoldedStackContainer{ public: -typedef TensorContainer<UGSTensor> _Ctype; -UnfoldedGContainer(const _Ctype*gs,int ngs,int nu) -:GContainer<UGSTensor> (gs,ngs,nu){} + typedef TensorContainer<UGSTensor> _Ctype; + UnfoldedGContainer(const _Ctype*gs,int ngs,int nu) + :GContainer<UGSTensor> (gs,ngs,nu){} }; @@ -396,112 +396,112 @@ UnfoldedGContainer(const _Ctype*gs,int ngs,int nu) template<class _Ttype> class StackProduct{ public: -typedef StackContainerInterface<_Ttype> _Stype; -typedef typename _Stype::_Ctype _Ctype; -typedef typename _Stype::itype itype; + typedef StackContainerInterface<_Ttype> _Stype; + typedef typename _Stype::_Ctype _Ctype; + typedef typename _Stype::itype itype; protected: -const _Stype&stack_cont; -InducedSymmetries syms; -Permutation per; + const _Stype&stack_cont; + InducedSymmetries syms; + Permutation per; public: -StackProduct(const _Stype&sc,const Equivalence&e, -const Symmetry&os) -:stack_cont(sc),syms(e,os),per(e){} -StackProduct(const _Stype&sc,const Equivalence&e, -const Permutation&p,const Symmetry&os) -:stack_cont(sc),syms(e,p,os),per(e,p){} -int dimen()const -{return syms.size();} -int getAllSize()const -{return stack_cont.getAllSize();} -const Symmetry&getProdSym(int ip)const -{return syms[ip];} -/*21:*/ - -bool isZero(const IntSequence&istacks)const -{ -TL_RAISE_IF(istacks.size()!=dimen(), -"Wrong istacks coordinates for StackProduct::isZero"); - -bool res= false; -int i= 0; -while(i<dimen()&&!(res= stack_cont.isZero(istacks[i],syms[i]))) -i++; -return res; -} - -/*:21*/ -; -/*22:*/ - -itype getType(int is,int ip)const -{ -TL_RAISE_IF(is<0||is>=stack_cont.numStacks(), -"Wrong index to stack in StackProduct::getType"); -TL_RAISE_IF(ip<0||ip>=dimen(), -"Wrong index to stack container in StackProduct::getType"); -return stack_cont.getType(is,syms[ip]); -} - -/*:22*/ -; -/*23:*/ - -const _Ttype*getMatrix(int is,int ip)const -{ -return stack_cont.getMatrix(is,syms[ip]); -} - -/*:23*/ -; -/*24:*/ - -void createPackedColumns(const IntSequence&coor, -Vector**vs,IntSequence&iu)const -{ -TL_RAISE_IF(iu.size()!=dimen(), -"Wrong storage length for unit flags in StackProduct::createPackedColumn"); -TL_RAISE_IF(coor.size()!=per.size(), -"Wrong size of index coor in StackProduct::createPackedColumn"); -IntSequence perindex(coor.size()); -per.apply(coor,perindex); -int off= 0; -for(int i= 0;i<dimen();i++){ -IntSequence percoor(perindex,off,syms[i].dimen()+off); -vs[i]= stack_cont.createPackedColumn(syms[i],percoor,iu[i]); -off+= syms[i].dimen(); -} -} - -/*:24*/ -; -/*25:*/ - -int getSize(int is)const -{ -return stack_cont.getStackSizes()[is]; -} - - -/*:25*/ -; -/*26:*/ - -int numMatrices(const IntSequence&istacks)const -{ -TL_RAISE_IF(istacks.size()!=dimen(), -"Wrong size of stack coordinates in StackContainer::numMatrices"); -int ret= 0; -int ip= 0; -while(ip<dimen()&&getType(istacks[ip],ip)==_Stype::matrix){ -ret++; -ip++; -} -return ret; -} - -/*:26*/ -; + StackProduct(const _Stype&sc,const Equivalence&e, + const Symmetry&os) + :stack_cont(sc),syms(e,os),per(e){} + StackProduct(const _Stype&sc,const Equivalence&e, + const Permutation&p,const Symmetry&os) + :stack_cont(sc),syms(e,p,os),per(e,p){} + int dimen()const + {return syms.size();} + int getAllSize()const + {return stack_cont.getAllSize();} + const Symmetry&getProdSym(int ip)const + {return syms[ip];} + /*21:*/ + + bool isZero(const IntSequence&istacks)const + { + TL_RAISE_IF(istacks.size()!=dimen(), + "Wrong istacks coordinates for StackProduct::isZero"); + + bool res= false; + int i= 0; + while(i<dimen()&&!(res= stack_cont.isZero(istacks[i],syms[i]))) + i++; + return res; + } + + /*:21*/ + ; + /*22:*/ + + itype getType(int is,int ip)const + { + TL_RAISE_IF(is<0||is>=stack_cont.numStacks(), + "Wrong index to stack in StackProduct::getType"); + TL_RAISE_IF(ip<0||ip>=dimen(), + "Wrong index to stack container in StackProduct::getType"); + return stack_cont.getType(is,syms[ip]); + } + + /*:22*/ + ; + /*23:*/ + + const _Ttype*getMatrix(int is,int ip)const + { + return stack_cont.getMatrix(is,syms[ip]); + } + + /*:23*/ + ; + /*24:*/ + + void createPackedColumns(const IntSequence&coor, + Vector**vs,IntSequence&iu)const + { + TL_RAISE_IF(iu.size()!=dimen(), + "Wrong storage length for unit flags in StackProduct::createPackedColumn"); + TL_RAISE_IF(coor.size()!=per.size(), + "Wrong size of index coor in StackProduct::createPackedColumn"); + IntSequence perindex(coor.size()); + per.apply(coor,perindex); + int off= 0; + for(int i= 0;i<dimen();i++){ + IntSequence percoor(perindex,off,syms[i].dimen()+off); + vs[i]= stack_cont.createPackedColumn(syms[i],percoor,iu[i]); + off+= syms[i].dimen(); + } + } + + /*:24*/ + ; + /*25:*/ + + int getSize(int is)const + { + return stack_cont.getStackSizes()[is]; + } + + + /*:25*/ + ; + /*26:*/ + + int numMatrices(const IntSequence&istacks)const + { + TL_RAISE_IF(istacks.size()!=dimen(), + "Wrong size of stack coordinates in StackContainer::numMatrices"); + int ret= 0; + int ip= 0; + while(ip<dimen()&&getType(istacks[ip],ip)==_Stype::matrix){ + ret++; + ip++; + } + return ret; + } + + /*:26*/ + ; }; /*:20*/ @@ -511,33 +511,33 @@ return ret; template<class _Ttype> class KronProdStack:public KronProdAllOptim{ public: -typedef StackProduct<_Ttype> _Ptype; -typedef StackContainerInterface<_Ttype> _Stype; -/*28:*/ - -KronProdStack(const _Ptype&sp,const IntSequence&istack) -:KronProdAllOptim(sp.dimen()) -{ -TL_RAISE_IF(sp.dimen()!=istack.size(), -"Wrong stack product dimension for KronProdStack constructor"); - -for(int i= 0;i<sp.dimen();i++){ -TL_RAISE_IF(sp.getType(istack[i],i)==_Stype::zero, -"Attempt to construct KronProdStack from zero matrix"); -if(sp.getType(istack[i],i)==_Stype::unit) -setUnit(i,sp.getSize(istack[i])); -if(sp.getType(istack[i],i)==_Stype::matrix){ -const TwoDMatrix*m= sp.getMatrix(istack[i],i); -TL_RAISE_IF(m->nrows()!=sp.getSize(istack[i]), -"Wrong size of returned matrix in KronProdStack constructor"); -setMat(i,*m); -} -} -} - - -/*:28*/ -; + typedef StackProduct<_Ttype> _Ptype; + typedef StackContainerInterface<_Ttype> _Stype; + /*28:*/ + + KronProdStack(const _Ptype&sp,const IntSequence&istack) + :KronProdAllOptim(sp.dimen()) + { + TL_RAISE_IF(sp.dimen()!=istack.size(), + "Wrong stack product dimension for KronProdStack constructor"); + + for(int i= 0;i<sp.dimen();i++){ + TL_RAISE_IF(sp.getType(istack[i],i)==_Stype::zero, + "Attempt to construct KronProdStack from zero matrix"); + if(sp.getType(istack[i],i)==_Stype::unit) + setUnit(i,sp.getSize(istack[i])); + if(sp.getType(istack[i],i)==_Stype::matrix){ + const TwoDMatrix*m= sp.getMatrix(istack[i],i); + TL_RAISE_IF(m->nrows()!=sp.getSize(istack[i]), + "Wrong size of returned matrix in KronProdStack constructor"); + setMat(i,*m); + } + } + } + + + /*:28*/ + ; }; /*:27*/ @@ -545,16 +545,16 @@ setMat(i,*m); /*29:*/ class WorkerFoldMAADense:public THREAD{ -const FoldedStackContainer&cont; -Symmetry sym; -const FGSContainer&dense_cont; -FGSTensor&out; + const FoldedStackContainer&cont; + Symmetry sym; + const FGSContainer&dense_cont; + FGSTensor&out; public: -WorkerFoldMAADense(const FoldedStackContainer&container, -const Symmetry&s, -const FGSContainer&dcontainer, -FGSTensor&outten); -void operator()(); + WorkerFoldMAADense(const FoldedStackContainer&container, + const Symmetry&s, + const FGSContainer&dcontainer, + FGSTensor&outten); + void operator()(); }; /*:29*/ @@ -562,16 +562,16 @@ void operator()(); /*30:*/ class WorkerFoldMAASparse1:public THREAD{ -const FoldedStackContainer&cont; -const FSSparseTensor&t; -FGSTensor&out; -IntSequence coor; -const EquivalenceBundle&ebundle; + const FoldedStackContainer&cont; + const FSSparseTensor&t; + FGSTensor&out; + IntSequence coor; + const EquivalenceBundle&ebundle; public: -WorkerFoldMAASparse1(const FoldedStackContainer&container, -const FSSparseTensor&ten, -FGSTensor&outten,const IntSequence&c); -void operator()(); + WorkerFoldMAASparse1(const FoldedStackContainer&container, + const FSSparseTensor&ten, + FGSTensor&outten,const IntSequence&c); + void operator()(); }; /*:30*/ @@ -579,15 +579,15 @@ void operator()(); /*31:*/ class WorkerFoldMAASparse2:public THREAD{ -const FoldedStackContainer&cont; -const FSSparseTensor&t; -FGSTensor&out; -IntSequence coor; + const FoldedStackContainer&cont; + const FSSparseTensor&t; + FGSTensor&out; + IntSequence coor; public: -WorkerFoldMAASparse2(const FoldedStackContainer&container, -const FSSparseTensor&ten, -FGSTensor&outten,const IntSequence&c); -void operator()(); + WorkerFoldMAASparse2(const FoldedStackContainer&container, + const FSSparseTensor&ten, + FGSTensor&outten,const IntSequence&c); + void operator()(); }; /*:31*/ @@ -595,15 +595,15 @@ void operator()(); /*32:*/ class WorkerFoldMAASparse4:public THREAD{ -const FoldedStackContainer&cont; -const FSSparseTensor&t; -FGSTensor&out; -IntSequence coor; + const FoldedStackContainer&cont; + const FSSparseTensor&t; + FGSTensor&out; + IntSequence coor; public: -WorkerFoldMAASparse4(const FoldedStackContainer&container, -const FSSparseTensor&ten, -FGSTensor&outten,const IntSequence&c); -void operator()(); + WorkerFoldMAASparse4(const FoldedStackContainer&container, + const FSSparseTensor&ten, + FGSTensor&outten,const IntSequence&c); + void operator()(); }; /*:32*/ @@ -611,16 +611,16 @@ void operator()(); /*33:*/ class WorkerUnfoldMAADense:public THREAD{ -const UnfoldedStackContainer&cont; -Symmetry sym; -const UGSContainer&dense_cont; -UGSTensor&out; + const UnfoldedStackContainer&cont; + Symmetry sym; + const UGSContainer&dense_cont; + UGSTensor&out; public: -WorkerUnfoldMAADense(const UnfoldedStackContainer&container, -const Symmetry&s, -const UGSContainer&dcontainer, -UGSTensor&outten); -void operator()(); + WorkerUnfoldMAADense(const UnfoldedStackContainer&container, + const Symmetry&s, + const UGSContainer&dcontainer, + UGSTensor&outten); + void operator()(); }; /*:33*/ @@ -628,16 +628,16 @@ void operator()(); /*34:*/ class WorkerUnfoldMAASparse1:public THREAD{ -const UnfoldedStackContainer&cont; -const FSSparseTensor&t; -UGSTensor&out; -IntSequence coor; -const EquivalenceBundle&ebundle; + const UnfoldedStackContainer&cont; + const FSSparseTensor&t; + UGSTensor&out; + IntSequence coor; + const EquivalenceBundle&ebundle; public: -WorkerUnfoldMAASparse1(const UnfoldedStackContainer&container, -const FSSparseTensor&ten, -UGSTensor&outten,const IntSequence&c); -void operator()(); + WorkerUnfoldMAASparse1(const UnfoldedStackContainer&container, + const FSSparseTensor&ten, + UGSTensor&outten,const IntSequence&c); + void operator()(); }; /*:34*/ @@ -645,15 +645,15 @@ void operator()(); /*35:*/ class WorkerUnfoldMAASparse2:public THREAD{ -const UnfoldedStackContainer&cont; -const FSSparseTensor&t; -UGSTensor&out; -IntSequence coor; + const UnfoldedStackContainer&cont; + const FSSparseTensor&t; + UGSTensor&out; + IntSequence coor; public: -WorkerUnfoldMAASparse2(const UnfoldedStackContainer&container, -const FSSparseTensor&ten, -UGSTensor&outten,const IntSequence&c); -void operator()(); + WorkerUnfoldMAASparse2(const UnfoldedStackContainer&container, + const FSSparseTensor&ten, + UGSTensor&outten,const IntSequence&c); + void operator()(); }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/sthread.h b/mex/sources/korderpert/Dyn_pp/tl/cc/sthread.h index 2e82be645455193e360c9cb42e3dae94ce4ae0ae..27964291e121c95ece5ae485727b2c925db23c3e 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/sthread.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/sthread.h @@ -12,464 +12,464 @@ #include <map> namespace sthread{ -using namespace std; - -class Empty{}; -/*2:*/ - -template<bool condition,class Then,class Else> -struct IF{ -typedef Then RET; -}; - -template<class Then,class Else> -struct IF<false,Then,Else> { -typedef Else RET; -}; - - - -/*:2*/ -; -enum{posix,empty}; -template<int> class thread_traits; -template<int> class detach_thread; -/*3:*/ - -template<int thread_impl> -class thread{ -typedef thread_traits<thread_impl> _Ttraits; -typedef typename _Ttraits::_Tthread _Tthread; -_Tthread th; -public: -virtual~thread(){} -_Tthread&getThreadIden() -{return th;} -const _Tthread&getThreadIden()const -{return th;} -virtual void operator()()= 0; -void run() -{_Ttraits::run(this);} -void detach_run() -{_Ttraits::detach_run(this);} -void exit() -{_Ttraits::exit();} -}; - -/*:3*/ -; -/*4:*/ - -template<int thread_impl> -class thread_group{ -typedef thread_traits<thread_impl> _Ttraits; -typedef thread<thread_impl> _Ctype; -list<_Ctype*> tlist; -typedef typename list<_Ctype*> ::iterator iterator; -public: -static int max_parallel_threads; -void insert(_Ctype*c) -{tlist.push_back(c);} -/*5:*/ - -~thread_group() -{ -while(!tlist.empty()){ -delete tlist.front(); -tlist.pop_front(); -} -} - -/*:5*/ -; -/*7:*/ - -void run() -{ -int rem= tlist.size(); -iterator pfirst= tlist.begin(); -while(rem> 2*max_parallel_threads){ -pfirst= run_portion(pfirst,max_parallel_threads); -rem-= max_parallel_threads; -} -if(rem> max_parallel_threads){ -pfirst= run_portion(pfirst,rem/2); -rem-= rem/2; -} -run_portion(pfirst,rem); -} - - - - -/*:7*/ -; -private: -/*6:*/ - -iterator run_portion(iterator start,int n) -{ -int c= 0; -for(iterator i= start;c<n;++i,c++){ -(*i)->run(); -} -iterator ret; -c= 0; -for(ret= start;c<n;++ret,c++){ -_Ttraits::join(*ret); -} -return ret; -} - - -/*:6*/ -; -}; - -/*:4*/ -; -/*8:*/ - -template<int thread_impl> -struct thread_traits{ -typedef typename IF<thread_impl==posix,pthread_t,Empty> ::RET _Tthread; -typedef thread<thread_impl> _Ctype; -typedef detach_thread<thread_impl> _Dtype; -static void run(_Ctype*c); -static void detach_run(_Dtype*c); -static void exit(); -static void join(_Ctype*c); -}; - -/*:8*/ -; -/*9:*/ - -struct ltmmkey; -typedef pair<const void*,const char*> mmkey; - -template<int thread_impl> -struct mutex_traits{ -typedef typename IF<thread_impl==posix,pthread_mutex_t,Empty> ::RET _Tmutex; -typedef map<mmkey,pair<_Tmutex,int> ,ltmmkey> mutex_int_map; -static void init(_Tmutex&m); -static void lock(_Tmutex&m); -static void unlock(_Tmutex&m); -}; - -/*:9*/ -; -/*10:*/ - -struct ltmmkey{ -bool operator()(const mmkey&k1,const mmkey&k2)const -{return k1.first<k2.first|| -(k1.first==k2.first&&strcmp(k1.second,k2.second)<0);} -}; - -template<int thread_impl> -class mutex_map -:public mutex_traits<thread_impl> ::mutex_int_map -{ -typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; -typedef mutex_traits<thread_impl> _Mtraits; -typedef pair<_Tmutex,int> mmval; -typedef map<mmkey,mmval,ltmmkey> _Tparent; -typedef typename _Tparent::iterator iterator; -typedef typename _Tparent::value_type _mvtype; -_Tmutex m; -public: -mutex_map() -{_Mtraits::init(m);} -void insert(const void*c,const char*id,const _Tmutex&m) -{_Tparent::insert(_mvtype(mmkey(c,id),mmval(m,0)));} -bool check(const void*c,const char*id)const -{return _Tparent::find(mmkey(c,id))!=_Tparent::end();} -/*11:*/ - -mmval*get(const void*c,const char*id) -{ -iterator it= _Tparent::find(mmkey(c,id)); -if(it==_Tparent::end()) -return NULL; -return&((*it).second); -} - -/*:11*/ -; -/*12:*/ - -void remove(const void*c,const char*id) -{ -iterator it= _Tparent::find(mmkey(c,id)); -if(it!=_Tparent::end()) -erase(it); -} - -/*:12*/ -; -void lock_map() -{_Mtraits::lock(m);} -void unlock_map() -{_Mtraits::unlock(m);} - -}; - -/*:10*/ -; -/*13:*/ - -template<int thread_impl> -class synchro{ -typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; -typedef mutex_traits<thread_impl> _Mtraits; -public: -typedef mutex_map<thread_impl> mutex_map_t; -private: -const void*caller; -const char*iden; -mutex_map_t&mutmap; -public: -synchro(const void*c,const char*id,mutex_map_t&mmap) -:caller(c),iden(id),mutmap(mmap) -{lock();} -~synchro() -{unlock();} -private: -/*14:*/ - -void lock(){ -mutmap.lock_map(); -if(!mutmap.check(caller,iden)){ -_Tmutex mut; -_Mtraits::init(mut); -mutmap.insert(caller,iden,mut); -} -mutmap.get(caller,iden)->second++; -mutmap.unlock_map(); -_Mtraits::lock(mutmap.get(caller,iden)->first); -} - -/*:14*/ -; -/*15:*/ - -void unlock(){ -mutmap.lock_map(); -if(mutmap.check(caller,iden)){ -_Mtraits::unlock(mutmap.get(caller,iden)->first); -mutmap.get(caller,iden)->second--; -if(mutmap.get(caller,iden)->second==0) -mutmap.remove(caller,iden); -} -mutmap.unlock_map(); -} - -/*:15*/ -; -}; - -/*:13*/ -; -/*16:*/ - -template<int thread_impl> -struct cond_traits{ -typedef typename IF<thread_impl==posix,pthread_cond_t,Empty> ::RET _Tcond; -typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; -static void init(_Tcond&cond); -static void broadcast(_Tcond&cond); -static void wait(_Tcond&cond,_Tmutex&mutex); -static void destroy(_Tcond&cond); -}; - -/*:16*/ -; -/*17:*/ - -template<int thread_impl> -class condition_counter{ -typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; -typedef typename cond_traits<thread_impl> ::_Tcond _Tcond; -int counter; -_Tmutex mut; -_Tcond cond; -bool changed; + using namespace std; + + class Empty{}; + /*2:*/ + + template<bool condition,class Then,class Else> + struct IF{ + typedef Then RET; + }; + + template<class Then,class Else> + struct IF<false,Then,Else> { + typedef Else RET; + }; + + + + /*:2*/ + ; + enum{posix,empty}; + template<int> class thread_traits; + template<int> class detach_thread; + /*3:*/ + + template<int thread_impl> + class thread{ + typedef thread_traits<thread_impl> _Ttraits; + typedef typename _Ttraits::_Tthread _Tthread; + _Tthread th; + public: + virtual~thread(){} + _Tthread&getThreadIden() + {return th;} + const _Tthread&getThreadIden()const + {return th;} + virtual void operator()()= 0; + void run() + {_Ttraits::run(this);} + void detach_run() + {_Ttraits::detach_run(this);} + void exit() + {_Ttraits::exit();} + }; + + /*:3*/ + ; + /*4:*/ + + template<int thread_impl> + class thread_group{ + typedef thread_traits<thread_impl> _Ttraits; + typedef thread<thread_impl> _Ctype; + list<_Ctype*> tlist; + typedef typename list<_Ctype*> ::iterator iterator; + public: + static int max_parallel_threads; + void insert(_Ctype*c) + {tlist.push_back(c);} + /*5:*/ + + ~thread_group() + { + while(!tlist.empty()){ + delete tlist.front(); + tlist.pop_front(); + } + } + + /*:5*/ + ; + /*7:*/ + + void run() + { + int rem= tlist.size(); + iterator pfirst= tlist.begin(); + while(rem> 2*max_parallel_threads){ + pfirst= run_portion(pfirst,max_parallel_threads); + rem-= max_parallel_threads; + } + if(rem> max_parallel_threads){ + pfirst= run_portion(pfirst,rem/2); + rem-= rem/2; + } + run_portion(pfirst,rem); + } + + + + + /*:7*/ + ; + private: + /*6:*/ + + iterator run_portion(iterator start,int n) + { + int c= 0; + for(iterator i= start;c<n;++i,c++){ + (*i)->run(); + } + iterator ret; + c= 0; + for(ret= start;c<n;++ret,c++){ + _Ttraits::join(*ret); + } + return ret; + } + + + /*:6*/ + ; + }; + + /*:4*/ + ; + /*8:*/ + + template<int thread_impl> + struct thread_traits{ + typedef typename IF<thread_impl==posix,pthread_t,Empty> ::RET _Tthread; + typedef thread<thread_impl> _Ctype; + typedef detach_thread<thread_impl> _Dtype; + static void run(_Ctype*c); + static void detach_run(_Dtype*c); + static void exit(); + static void join(_Ctype*c); + }; + + /*:8*/ + ; + /*9:*/ + + struct ltmmkey; + typedef pair<const void*,const char*> mmkey; + + template<int thread_impl> + struct mutex_traits{ + typedef typename IF<thread_impl==posix,pthread_mutex_t,Empty> ::RET _Tmutex; + typedef map<mmkey,pair<_Tmutex,int> ,ltmmkey> mutex_int_map; + static void init(_Tmutex&m); + static void lock(_Tmutex&m); + static void unlock(_Tmutex&m); + }; + + /*:9*/ + ; + /*10:*/ + + struct ltmmkey{ + bool operator()(const mmkey&k1,const mmkey&k2)const + {return k1.first<k2.first|| + (k1.first==k2.first&&strcmp(k1.second,k2.second)<0);} + }; + + template<int thread_impl> + class mutex_map + :public mutex_traits<thread_impl> ::mutex_int_map + { + typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; + typedef mutex_traits<thread_impl> _Mtraits; + typedef pair<_Tmutex,int> mmval; + typedef map<mmkey,mmval,ltmmkey> _Tparent; + typedef typename _Tparent::iterator iterator; + typedef typename _Tparent::value_type _mvtype; + _Tmutex m; + public: + mutex_map() + {_Mtraits::init(m);} + void insert(const void*c,const char*id,const _Tmutex&m) + {_Tparent::insert(_mvtype(mmkey(c,id),mmval(m,0)));} + bool check(const void*c,const char*id)const + {return _Tparent::find(mmkey(c,id))!=_Tparent::end();} + /*11:*/ + + mmval*get(const void*c,const char*id) + { + iterator it= _Tparent::find(mmkey(c,id)); + if(it==_Tparent::end()) + return NULL; + return&((*it).second); + } + + /*:11*/ + ; + /*12:*/ + + void remove(const void*c,const char*id) + { + iterator it= _Tparent::find(mmkey(c,id)); + if(it!=_Tparent::end()) + erase(it); + } + + /*:12*/ + ; + void lock_map() + {_Mtraits::lock(m);} + void unlock_map() + {_Mtraits::unlock(m);} + + }; + + /*:10*/ + ; + /*13:*/ + + template<int thread_impl> + class synchro{ + typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; + typedef mutex_traits<thread_impl> _Mtraits; + public: + typedef mutex_map<thread_impl> mutex_map_t; + private: + const void*caller; + const char*iden; + mutex_map_t&mutmap; + public: + synchro(const void*c,const char*id,mutex_map_t&mmap) + :caller(c),iden(id),mutmap(mmap) + {lock();} + ~synchro() + {unlock();} + private: + /*14:*/ + + void lock(){ + mutmap.lock_map(); + if(!mutmap.check(caller,iden)){ + _Tmutex mut; + _Mtraits::init(mut); + mutmap.insert(caller,iden,mut); + } + mutmap.get(caller,iden)->second++; + mutmap.unlock_map(); + _Mtraits::lock(mutmap.get(caller,iden)->first); + } + + /*:14*/ + ; + /*15:*/ + + void unlock(){ + mutmap.lock_map(); + if(mutmap.check(caller,iden)){ + _Mtraits::unlock(mutmap.get(caller,iden)->first); + mutmap.get(caller,iden)->second--; + if(mutmap.get(caller,iden)->second==0) + mutmap.remove(caller,iden); + } + mutmap.unlock_map(); + } + + /*:15*/ + ; + }; + + /*:13*/ + ; + /*16:*/ + + template<int thread_impl> + struct cond_traits{ + typedef typename IF<thread_impl==posix,pthread_cond_t,Empty> ::RET _Tcond; + typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; + static void init(_Tcond&cond); + static void broadcast(_Tcond&cond); + static void wait(_Tcond&cond,_Tmutex&mutex); + static void destroy(_Tcond&cond); + }; + + /*:16*/ + ; + /*17:*/ + + template<int thread_impl> + class condition_counter{ + typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex; + typedef typename cond_traits<thread_impl> ::_Tcond _Tcond; + int counter; + _Tmutex mut; + _Tcond cond; + bool changed; public: -/*18:*/ - -condition_counter() -:counter(0),changed(true) -{ -mutex_traits<thread_impl> ::init(mut); -cond_traits<thread_impl> ::init(cond); -} - -/*:18*/ -; -/*19:*/ - -~condition_counter() -{ -cond_traits<thread_impl> ::destroy(cond); -} - -/*:19*/ -; -/*20:*/ - -void increase() -{ -mutex_traits<thread_impl> ::lock(mut); -counter++; -changed= true; -cond_traits<thread_impl> ::broadcast(cond); -mutex_traits<thread_impl> ::unlock(mut); -} - -/*:20*/ -; -/*21:*/ - -void decrease() -{ -mutex_traits<thread_impl> ::lock(mut); -counter--; -changed= true; -cond_traits<thread_impl> ::broadcast(cond); -mutex_traits<thread_impl> ::unlock(mut); -} - -/*:21*/ -; -/*22:*/ - -int waitForChange() -{ -mutex_traits<thread_impl> ::lock(mut); -if(!changed){ -cond_traits<thread_impl> ::wait(cond,mut); -} -changed= false; -int res= counter; -mutex_traits<thread_impl> ::unlock(mut); -return res; -} - - -/*:22*/ -; -}; - -/*:17*/ -; -/*23:*/ - -template<int thread_impl> -class detach_thread:public thread<thread_impl> { + /*18:*/ + + condition_counter() + :counter(0),changed(true) + { + mutex_traits<thread_impl> ::init(mut); + cond_traits<thread_impl> ::init(cond); + } + + /*:18*/ + ; + /*19:*/ + + ~condition_counter() + { + cond_traits<thread_impl> ::destroy(cond); + } + + /*:19*/ + ; + /*20:*/ + + void increase() + { + mutex_traits<thread_impl> ::lock(mut); + counter++; + changed= true; + cond_traits<thread_impl> ::broadcast(cond); + mutex_traits<thread_impl> ::unlock(mut); + } + + /*:20*/ + ; + /*21:*/ + + void decrease() + { + mutex_traits<thread_impl> ::lock(mut); + counter--; + changed= true; + cond_traits<thread_impl> ::broadcast(cond); + mutex_traits<thread_impl> ::unlock(mut); + } + + /*:21*/ + ; + /*22:*/ + + int waitForChange() + { + mutex_traits<thread_impl> ::lock(mut); + if(!changed){ + cond_traits<thread_impl> ::wait(cond,mut); + } + changed= false; + int res= counter; + mutex_traits<thread_impl> ::unlock(mut); + return res; + } + + + /*:22*/ + ; + }; + + /*:17*/ + ; + /*23:*/ + + template<int thread_impl> + class detach_thread:public thread<thread_impl> { public: -condition_counter<thread_impl> *counter; -detach_thread():counter(NULL){} -void installCounter(condition_counter<thread_impl> *c) -{counter= c;} -void run() -{thread_traits<thread_impl> ::detach_run(this);} -}; - -/*:23*/ -; -/*24:*/ - -template<int thread_impl> -class detach_thread_group{ -typedef thread_traits<thread_impl> _Ttraits; -typedef cond_traits<thread_impl> _Ctraits; -typedef detach_thread<thread_impl> _Ctype; -list<_Ctype*> tlist; -typedef typename list<_Ctype*> ::iterator iterator; -condition_counter<thread_impl> counter; -public: -static int max_parallel_threads; -/*25:*/ - -void insert(_Ctype*c) -{ -tlist.push_back(c); -c->installCounter(&counter); -} - -/*:25*/ -; -/*26:*/ - -~detach_thread_group() -{ -while(!tlist.empty()){ -delete tlist.front(); -tlist.pop_front(); -} -} - -/*:26*/ -; -/*27:*/ - -void run() -{ -int mpt= max_parallel_threads; -iterator it= tlist.begin(); -while(it!=tlist.end()){ -if(counter.waitForChange()<mpt){ -counter.increase(); -(*it)->run(); -++it; -} -} -while(counter.waitForChange()> 0){} -} - - -/*:27*/ -; -}; - -/*:24*/ -; + condition_counter<thread_impl> *counter; + detach_thread():counter(NULL){} + void installCounter(condition_counter<thread_impl> *c) + {counter= c;} + void run() + {thread_traits<thread_impl> ::detach_run(this);} + }; + + /*:23*/ + ; + /*24:*/ + + template<int thread_impl> + class detach_thread_group{ + typedef thread_traits<thread_impl> _Ttraits; + typedef cond_traits<thread_impl> _Ctraits; + typedef detach_thread<thread_impl> _Ctype; + list<_Ctype*> tlist; + typedef typename list<_Ctype*> ::iterator iterator; + condition_counter<thread_impl> counter; + public: + static int max_parallel_threads; + /*25:*/ + + void insert(_Ctype*c) + { + tlist.push_back(c); + c->installCounter(&counter); + } + + /*:25*/ + ; + /*26:*/ + + ~detach_thread_group() + { + while(!tlist.empty()){ + delete tlist.front(); + tlist.pop_front(); + } + } + + /*:26*/ + ; + /*27:*/ + + void run() + { + int mpt= max_parallel_threads; + iterator it= tlist.begin(); + while(it!=tlist.end()){ + if(counter.waitForChange()<mpt){ + counter.increase(); + (*it)->run(); + ++it; + } + } + while(counter.waitForChange()> 0){} + } + + + /*:27*/ + ; + }; + + /*:24*/ + ; #ifdef POSIX_THREADS -/*28:*/ - -typedef detach_thread<posix> PosixThread; -typedef detach_thread_group<posix> PosixThreadGroup; -typedef synchro<posix> posix_synchro; -class PosixSynchro:public posix_synchro{ -public: -PosixSynchro(const void*c,const char*id); -}; - + /*28:*/ + + typedef detach_thread<posix> PosixThread; + typedef detach_thread_group<posix> PosixThreadGroup; + typedef synchro<posix> posix_synchro; + class PosixSynchro:public posix_synchro{ + public: + PosixSynchro(const void*c,const char*id); + }; + #define THREAD sthread::PosixThread #define THREAD_GROUP sthread::PosixThreadGroup #define SYNCHRO sthread::PosixSynchro - -/*:28*/ -; + + /*:28*/ + ; #else -/*29:*/ - -typedef thread<empty> NoThread; -typedef thread_group<empty> NoThreadGroup; -typedef synchro<empty> no_synchro; -class NoSynchro{ -public: -NoSynchro(const void*c,const char*id){} -~NoSynchro(){} -}; - + /*29:*/ + + typedef thread<empty> NoThread; + typedef thread_group<empty> NoThreadGroup; + typedef synchro<empty> no_synchro; + class NoSynchro{ + public: + NoSynchro(const void*c,const char*id){} + ~NoSynchro(){} + }; + #define THREAD sthread::NoThread #define THREAD_GROUP sthread::NoThreadGroup #define SYNCHRO sthread::NoSynchro - -/*:29*/ -; + + /*:29*/ + ; #endif }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/symmetry.h b/mex/sources/korderpert/Dyn_pp/tl/cc/symmetry.h index 8700b33f4fad41a180f42e89f5272cdb696bdfc2..3a5c8447c9cd271a7c3a1683db914546622ffb3b 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/symmetry.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/symmetry.h @@ -13,45 +13,45 @@ class Symmetry:public IntSequence{ public: -/*3:*/ - -Symmetry(int len,const char*dummy) -:IntSequence(len,0){} -Symmetry(int i1) -:IntSequence(1,i1){} -Symmetry(int i1,int i2) -:IntSequence(2){operator[](0)= i1;operator[](1)= i2;} -Symmetry(int i1,int i2,int i3) -:IntSequence(3) -{ -operator[](0)= i1; -operator[](1)= i2; -operator[](2)= i3; -} -Symmetry(int i1,int i2,int i3,int i4) -:IntSequence(4) -{ -operator[](0)= i1; -operator[](1)= i2; -operator[](2)= i3; -operator[](3)= i4; -} -Symmetry(const Symmetry&s) -:IntSequence(s){} -Symmetry(const Symmetry&s,const OrdSequence&cl) -:IntSequence(s,cl.getData()){} -Symmetry(Symmetry&s,int len) -:IntSequence(s,s.size()-len,s.size()){} -Symmetry(const IntSequence&s); - -/*:3*/ -; -int num()const -{return size();} -int dimen()const -{return sum();} -int findClass(int i)const; -bool isFull()const; + /*3:*/ + + Symmetry(int len,const char*dummy) + :IntSequence(len,0){} + Symmetry(int i1) + :IntSequence(1,i1){} + Symmetry(int i1,int i2) + :IntSequence(2){operator[](0)= i1;operator[](1)= i2;} + Symmetry(int i1,int i2,int i3) + :IntSequence(3) + { + operator[](0)= i1; + operator[](1)= i2; + operator[](2)= i3; + } + Symmetry(int i1,int i2,int i3,int i4) + :IntSequence(4) + { + operator[](0)= i1; + operator[](1)= i2; + operator[](2)= i3; + operator[](3)= i4; + } + Symmetry(const Symmetry&s) + :IntSequence(s){} + Symmetry(const Symmetry&s,const OrdSequence&cl) + :IntSequence(s,cl.getData()){} + Symmetry(Symmetry&s,int len) + :IntSequence(s,s.size()-len,s.size()){} + Symmetry(const IntSequence&s); + + /*:3*/ + ; + int num()const + {return size();} + int dimen()const + {return sum();} + int findClass(int i)const; + bool isFull()const; }; /*:2*/ @@ -59,21 +59,21 @@ bool isFull()const; /*4:*/ class SymmetrySet{ -Symmetry run; -int dim; + Symmetry run; + int dim; public: -SymmetrySet(int d,int length) -:run(length,""),dim(d){} -SymmetrySet(SymmetrySet&s,int d) -:run(s.run,s.size()-1),dim(d){} -int dimen()const -{return dim;} -const Symmetry&sym()const -{return run;} -Symmetry&sym() -{return run;} -int size()const -{return run.size();} + SymmetrySet(int d,int length) + :run(length,""),dim(d){} + SymmetrySet(SymmetrySet&s,int d) + :run(s.run,s.size()-1),dim(d){} + int dimen()const + {return dim;} + const Symmetry&sym()const + {return run;} + Symmetry&sym() + {return run;} + int size()const + {return run.size();} }; /*:4*/ @@ -81,18 +81,18 @@ int size()const /*5:*/ class symiterator{ -SymmetrySet&s; -symiterator*subit; -SymmetrySet*subs; -bool end_flag; + SymmetrySet&s; + symiterator*subit; + SymmetrySet*subs; + bool end_flag; public: -symiterator(SymmetrySet&ss); -~symiterator(); -symiterator&operator++(); -bool isEnd()const -{return end_flag;} -const Symmetry&operator*()const -{return s.sym();} + symiterator(SymmetrySet&ss); + ~symiterator(); + symiterator&operator++(); + bool isEnd()const + {return end_flag;} + const Symmetry&operator*()const + {return s.sym();} }; @@ -102,9 +102,9 @@ const Symmetry&operator*()const class InducedSymmetries:public vector<Symmetry> { public: -InducedSymmetries(const Equivalence&e,const Symmetry&s); -InducedSymmetries(const Equivalence&e,const Permutation&p,const Symmetry&s); -void print()const; + InducedSymmetries(const Equivalence&e,const Symmetry&s); + InducedSymmetries(const Equivalence&e,const Permutation&p,const Symmetry&s); + void print()const; }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/t_container.h b/mex/sources/korderpert/Dyn_pp/tl/cc/t_container.h index 392c27be1167b274acb07c2c9ec24629bda6474e..3ec7f86d9f041dca1b56dd3fb5af040619b7afc6 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/t_container.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/t_container.h @@ -17,8 +17,8 @@ /*2:*/ struct ltsym{ -bool operator()(const Symmetry&s1,const Symmetry&s2)const -{return s1<s2;} + bool operator()(const Symmetry&s1,const Symmetry&s2)const + {return s1<s2;} }; /*:2*/ @@ -27,224 +27,224 @@ bool operator()(const Symmetry&s1,const Symmetry&s2)const template<class _Ttype> class TensorContainer{ protected: -typedef const _Ttype*_const_ptr; -typedef _Ttype*_ptr; -typedef map<Symmetry,_ptr,ltsym> _Map; -typedef typename _Map::value_type _mvtype; + typedef const _Ttype*_const_ptr; + typedef _Ttype*_ptr; + typedef map<Symmetry,_ptr,ltsym> _Map; + typedef typename _Map::value_type _mvtype; public: -typedef typename _Map::iterator iterator; -typedef typename _Map::const_iterator const_iterator; + typedef typename _Map::iterator iterator; + typedef typename _Map::const_iterator const_iterator; private: -int n; -_Map m; + int n; + _Map m; protected: -const EquivalenceBundle&ebundle; + const EquivalenceBundle&ebundle; public: -TensorContainer(int nn) -:n(nn),ebundle(*(tls.ebundle)){} -/*5:*/ - -TensorContainer(const TensorContainer<_Ttype> &c) -:n(c.n),m(),ebundle(c.ebundle) -{ -for(const_iterator it= c.m.begin();it!=c.m.end();++it){ -_Ttype*ten= new _Ttype(*((*it).second)); -insert(ten); -} -} - -/*:5*/ -; -/*6:*/ - -TensorContainer(int first_row,int num,TensorContainer<_Ttype> &c) -:n(c.n),ebundle(*(tls.ebundle)) -{ -for(iterator it= c.m.begin();it!=c.m.end();++it){ -_Ttype*t= new _Ttype(first_row,num,*((*it).second)); -insert(t); -} -} - - -/*:6*/ -; -/*7:*/ - -_const_ptr get(const Symmetry&s)const -{ -TL_RAISE_IF(s.num()!=num(), -"Incompatible symmetry lookup in TensorContainer::get"); -const_iterator it= m.find(s); -if(it==m.end()){ -TL_RAISE("Symmetry not found in TensorContainer::get"); -return NULL; -}else{ -return(*it).second; -} -} - - -_ptr get(const Symmetry&s) -{ -TL_RAISE_IF(s.num()!=num(), -"Incompatible symmetry lookup in TensorContainer::get"); -iterator it= m.find(s); -if(it==m.end()){ -TL_RAISE("Symmetry not found in TensorContainer::get"); -return NULL; -}else{ -return(*it).second; -} -} - -/*:7*/ -; -/*8:*/ - -bool check(const Symmetry&s)const -{ -TL_RAISE_IF(s.num()!=num(), -"Incompatible symmetry lookup in TensorContainer::check"); -const_iterator it= m.find(s); -return it!=m.end(); -} - -/*:8*/ -; -/*9:*/ - -void insert(_ptr t) -{ -TL_RAISE_IF(t->getSym().num()!=num(), -"Incompatible symmetry insertion in TensorContainer::insert"); -TL_RAISE_IF(check(t->getSym()), -"Tensor already in container in TensorContainer::insert"); -m.insert(_mvtype(t->getSym(),t)); -if(!t->isFinite()){ -throw TLException(__FILE__,__LINE__,"NaN or Inf asserted in TensorContainer::insert"); -} -} - -/*:9*/ -; -/*10:*/ - -void remove(const Symmetry&s) -{ -iterator it= m.find(s); -if(it!=m.end()){ -_ptr t= (*it).second; -m.erase(it); -delete t; -} -} - - -/*:10*/ -; -/*11:*/ - -void clear() -{ -while(!m.empty()){ -delete(*(m.begin())).second; -m.erase(m.begin()); -} -} - -/*:11*/ -; -/*15:*/ - -vector<_const_ptr> -fetchTensors(const Symmetry&rsym,const Equivalence&e)const -{ -vector<_const_ptr> res(e.numClasses()); -int i= 0; -for(Equivalence::const_seqit it= e.begin(); -it!=e.end();++it,i++){ -Symmetry s(rsym,*it); -res[i]= get(s); -} -return res; -} - -/*:15*/ -; -/*12:*/ - -int getMaxDim()const -{ -int res= -1; -for(const_iterator run= m.begin();run!=m.end();++run){ -int dim= (*run).first.dimen(); -if(dim> res) -res= dim; -} -return res; -} - - -/*:12*/ -; -/*13:*/ - -void print()const -{ -printf("Tensor container: nvars=%d, tensors=%d\n",n,m.size()); -for(const_iterator it= m.begin();it!=m.end();++it){ -printf("Symmetry: "); -(*it).first.print(); -((*it).second)->print(); -} -} - -/*:13*/ -; -/*14:*/ - -void writeMat4(FILE*fd,const char*prefix)const -{ -for(const_iterator it= begin();it!=end();++it){ -char lname[100]; -sprintf(lname,"%s_g",prefix); -const Symmetry&sym= (*it).first; -for(int i= 0;i<sym.num();i++){ -char tmp[10]; -sprintf(tmp,"_%d",sym[i]); -strcat(lname,tmp); -} -ConstTwoDMatrix m(*((*it).second)); -m.writeMat4(fd,lname); -} -} - - -/*:14*/ -; - -virtual~TensorContainer() -{clear();} - -/*4:*/ - -int num()const -{return n;} -const EquivalenceBundle&getEqBundle()const -{return ebundle;} - -const_iterator begin()const -{return m.begin();} -const_iterator end()const -{return m.end();} -iterator begin() -{return m.begin();} -iterator end() -{return m.end();} - -/*:4*/ -; + TensorContainer(int nn) + :n(nn),ebundle(*(tls.ebundle)){} + /*5:*/ + + TensorContainer(const TensorContainer<_Ttype> &c) + :n(c.n),m(),ebundle(c.ebundle) + { + for(const_iterator it= c.m.begin();it!=c.m.end();++it){ + _Ttype*ten= new _Ttype(*((*it).second)); + insert(ten); + } + } + + /*:5*/ + ; + /*6:*/ + + TensorContainer(int first_row,int num,TensorContainer<_Ttype> &c) + :n(c.n),ebundle(*(tls.ebundle)) + { + for(iterator it= c.m.begin();it!=c.m.end();++it){ + _Ttype*t= new _Ttype(first_row,num,*((*it).second)); + insert(t); + } + } + + + /*:6*/ + ; + /*7:*/ + + _const_ptr get(const Symmetry&s)const + { + TL_RAISE_IF(s.num()!=num(), + "Incompatible symmetry lookup in TensorContainer::get"); + const_iterator it= m.find(s); + if(it==m.end()){ + TL_RAISE("Symmetry not found in TensorContainer::get"); + return NULL; + }else{ + return(*it).second; + } + } + + + _ptr get(const Symmetry&s) + { + TL_RAISE_IF(s.num()!=num(), + "Incompatible symmetry lookup in TensorContainer::get"); + iterator it= m.find(s); + if(it==m.end()){ + TL_RAISE("Symmetry not found in TensorContainer::get"); + return NULL; + }else{ + return(*it).second; + } + } + + /*:7*/ + ; + /*8:*/ + + bool check(const Symmetry&s)const + { + TL_RAISE_IF(s.num()!=num(), + "Incompatible symmetry lookup in TensorContainer::check"); + const_iterator it= m.find(s); + return it!=m.end(); + } + + /*:8*/ + ; + /*9:*/ + + void insert(_ptr t) + { + TL_RAISE_IF(t->getSym().num()!=num(), + "Incompatible symmetry insertion in TensorContainer::insert"); + TL_RAISE_IF(check(t->getSym()), + "Tensor already in container in TensorContainer::insert"); + m.insert(_mvtype(t->getSym(),t)); + if(!t->isFinite()){ + throw TLException(__FILE__,__LINE__,"NaN or Inf asserted in TensorContainer::insert"); + } + } + + /*:9*/ + ; + /*10:*/ + + void remove(const Symmetry&s) + { + iterator it= m.find(s); + if(it!=m.end()){ + _ptr t= (*it).second; + m.erase(it); + delete t; + } + } + + + /*:10*/ + ; + /*11:*/ + + void clear() + { + while(!m.empty()){ + delete(*(m.begin())).second; + m.erase(m.begin()); + } + } + + /*:11*/ + ; + /*15:*/ + + vector<_const_ptr> + fetchTensors(const Symmetry&rsym,const Equivalence&e)const + { + vector<_const_ptr> res(e.numClasses()); + int i= 0; + for(Equivalence::const_seqit it= e.begin(); + it!=e.end();++it,i++){ + Symmetry s(rsym,*it); + res[i]= get(s); + } + return res; + } + + /*:15*/ + ; + /*12:*/ + + int getMaxDim()const + { + int res= -1; + for(const_iterator run= m.begin();run!=m.end();++run){ + int dim= (*run).first.dimen(); + if(dim> res) + res= dim; + } + return res; + } + + + /*:12*/ + ; + /*13:*/ + + void print()const + { + printf("Tensor container: nvars=%d, tensors=%d\n",n,m.size()); + for(const_iterator it= m.begin();it!=m.end();++it){ + printf("Symmetry: "); + (*it).first.print(); + ((*it).second)->print(); + } + } + + /*:13*/ + ; + /*14:*/ + + void writeMat4(FILE*fd,const char*prefix)const + { + for(const_iterator it= begin();it!=end();++it){ + char lname[100]; + sprintf(lname,"%s_g",prefix); + const Symmetry&sym= (*it).first; + for(int i= 0;i<sym.num();i++){ + char tmp[10]; + sprintf(tmp,"_%d",sym[i]); + strcat(lname,tmp); + } + ConstTwoDMatrix m(*((*it).second)); + m.writeMat4(fd,lname); + } + } + + + /*:14*/ + ; + + virtual~TensorContainer() + {clear();} + + /*4:*/ + + int num()const + {return n;} + const EquivalenceBundle&getEqBundle()const + {return ebundle;} + + const_iterator begin()const + {return m.begin();} + const_iterator end()const + {return m.end();} + iterator begin() + {return m.begin();} + iterator end() + {return m.end();} + + /*:4*/ + ; }; /*:3*/ @@ -254,12 +254,12 @@ iterator end() class FGSContainer; class UGSContainer:public TensorContainer<UGSTensor> { public: -UGSContainer(int nn) -:TensorContainer<UGSTensor> (nn){} -UGSContainer(const UGSContainer&uc) -:TensorContainer<UGSTensor> (uc){} -UGSContainer(const FGSContainer&c); -void multAndAdd(const UGSTensor&t,UGSTensor&out)const; + UGSContainer(int nn) + :TensorContainer<UGSTensor> (nn){} + UGSContainer(const UGSContainer&uc) + :TensorContainer<UGSTensor> (uc){} + UGSContainer(const FGSContainer&c); + void multAndAdd(const UGSTensor&t,UGSTensor&out)const; }; @@ -268,20 +268,20 @@ void multAndAdd(const UGSTensor&t,UGSTensor&out)const; /*17:*/ class FGSContainer:public TensorContainer<FGSTensor> { -static const int num_one_time; + static const int num_one_time; public: -FGSContainer(int nn) -:TensorContainer<FGSTensor> (nn){} -FGSContainer(const FGSContainer&fc) -:TensorContainer<FGSTensor> (fc){} -FGSContainer(const UGSContainer&c); -void multAndAdd(const FGSTensor&t,FGSTensor&out)const; -void multAndAdd(const UGSTensor&t,FGSTensor&out)const; + FGSContainer(int nn) + :TensorContainer<FGSTensor> (nn){} + FGSContainer(const FGSContainer&fc) + :TensorContainer<FGSTensor> (fc){} + FGSContainer(const UGSContainer&c); + void multAndAdd(const FGSTensor&t,FGSTensor&out)const; + void multAndAdd(const UGSTensor&t,FGSTensor&out)const; private: -static Tensor::index -getIndices(int num,vector<IntSequence> &out, -const Tensor::index&start, -const Tensor::index&end); + static Tensor::index + getIndices(int num,vector<IntSequence> &out, + const Tensor::index&start, + const Tensor::index&end); }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/t_polynomial.h b/mex/sources/korderpert/Dyn_pp/tl/cc/t_polynomial.h index 535ad7818c68c66c7311b048bcb96d18e9d33514..84e8739df7dd92cb9592eba0fcfb7c9037ee6e43 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/t_polynomial.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/t_polynomial.h @@ -8,16 +8,16 @@ /*2:*/ class PowerProvider{ -Vector origv; -URSingleTensor*ut; -FRSingleTensor*ft; -int nv; + Vector origv; + URSingleTensor*ut; + FRSingleTensor*ft; + int nv; public: -PowerProvider(const ConstVector&v) -:origv(v),ut(NULL),ft(NULL),nv(v.length()){} -~PowerProvider(); -const URSingleTensor&getNext(const URSingleTensor*dummy); -const FRSingleTensor&getNext(const FRSingleTensor*dummy); + PowerProvider(const ConstVector&v) + :origv(v),ut(NULL),ft(NULL),nv(v.length()){} + ~PowerProvider(); + const URSingleTensor&getNext(const URSingleTensor*dummy); + const FRSingleTensor&getNext(const FRSingleTensor*dummy); }; /*:2*/ @@ -26,226 +26,226 @@ const FRSingleTensor&getNext(const FRSingleTensor*dummy); template<class _Ttype,class _TGStype,class _Stype> class TensorPolynomial:public TensorContainer<_Ttype> { -int nr; -int nv; -int maxdim; -typedef TensorContainer<_Ttype> _Tparent; -typedef typename _Tparent::_ptr _ptr; + int nr; + int nv; + int maxdim; + typedef TensorContainer<_Ttype> _Tparent; + typedef typename _Tparent::_ptr _ptr; public: -TensorPolynomial(int rows,int vars) -:TensorContainer<_Ttype> (1), -nr(rows),nv(vars),maxdim(0){} -TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,int k) -:TensorContainer<_Ttype> (tp), -nr(tp.nr),nv(tp.nv),maxdim(0){derivative(k);} -TensorPolynomial(int first_row,int num,TensorPolynomial<_Ttype,_TGStype,_Stype> &tp) -:TensorContainer<_Ttype> (first_row,num,tp), -nr(num),nv(tp.nv),maxdim(tp.maxdim){} -/*4:*/ - -TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,const Vector&xval) -:TensorContainer<_Ttype> (1), -nr(tp.nrows()),nv(tp.nvars()-xval.length()),maxdim(0) -{ -TL_RAISE_IF(nvars()<0, -"Length of xval too big in TensorPolynomial contract constructor"); -IntSequence ss(2);ss[0]= xval.length();ss[1]= nvars(); -IntSequence pp(2);pp[0]= 0;pp[1]= 1; - -/*5:*/ - -PowerProvider pwp(xval); -for(int i= 1;i<=tp.maxdim;i++){ -const _Stype&xpow= pwp.getNext((const _Stype*)NULL); -for(int j= 0;j<=tp.maxdim-i;j++){ -if(tp.check(Symmetry(i+j))){ -/*7:*/ - -_Ttype*ten; -if(_Tparent::check(Symmetry(j))){ -ten= _Tparent::get(Symmetry(j)); -}else{ -ten= new _Ttype(nrows(),nvars(),j); -ten->zeros(); -insert(ten); -} - - -/*:7*/ -; -Symmetry sym(i,j); -IntSequence coor(sym,pp); -_TGStype slice(*(tp.get(Symmetry(i+j))),ss,coor,TensorDimens(sym,ss)); -slice.mult(Tensor::noverk(i+j,j)); -_TGStype tmp(*ten); -slice.contractAndAdd(0,tmp,xpow); -} -} -} - -/*:5*/ -; -/*6:*/ - -for(int j= 0;j<=tp.maxdim;j++){ -if(tp.check(Symmetry(j))){ -/*7:*/ - -_Ttype*ten; -if(_Tparent::check(Symmetry(j))){ -ten= _Tparent::get(Symmetry(j)); -}else{ -ten= new _Ttype(nrows(),nvars(),j); -ten->zeros(); -insert(ten); -} - - -/*:7*/ -; -Symmetry sym(0,j); -IntSequence coor(sym,pp); -_TGStype slice(*(tp.get(Symmetry(j))),ss,coor,TensorDimens(sym,ss)); -ten->add(1.0,slice); -} -} - - -/*:6*/ -; -} - -/*:4*/ -; -TensorPolynomial(const TensorPolynomial&tp) -:TensorContainer<_Ttype> (tp),nr(tp.nr),nv(tp.nv),maxdim(tp.maxdim){} -int nrows()const -{return nr;} -int nvars()const -{return nv;} -/*8:*/ - -void evalTrad(Vector&out,const ConstVector&v)const -{ -if(_Tparent::check(Symmetry(0))) -out= _Tparent::get(Symmetry(0))->getData(); -else -out.zeros(); - -PowerProvider pp(v); -for(int d= 1;d<=maxdim;d++){ -const _Stype&p= pp.getNext((const _Stype*)NULL); -Symmetry cs(d); -if(_Tparent::check(cs)){ -const _Ttype*t= _Tparent::get(cs); -t->multaVec(out,p.getData()); -} -} -} - -/*:8*/ -; -/*9:*/ - -void evalHorner(Vector&out,const ConstVector&v)const -{ -if(_Tparent::check(Symmetry(0))) -out= _Tparent::get(Symmetry(0))->getData(); -else -out.zeros(); - -if(maxdim==0) -return; - -_Ttype*last; -if(maxdim==1) -last= new _Ttype(*(_Tparent::get(Symmetry(1)))); -else -last= new _Ttype(*(_Tparent::get(Symmetry(maxdim))),v); -for(int d= maxdim-1;d>=1;d--){ -Symmetry cs(d); -if(_Tparent::check(cs)){ -const _Ttype*nt= _Tparent::get(cs); -last->add(1.0,ConstTwoDMatrix(*nt)); -} -if(d> 1){ -_Ttype*new_last= new _Ttype(*last,v); -delete last; -last= new_last; -} -} -last->multaVec(out,v); -delete last; -} - -/*:9*/ -; -/*10:*/ - -void insert(_ptr t) -{ -TL_RAISE_IF(t->nrows()!=nr, -"Wrong number of rows in TensorPolynomial::insert"); -TL_RAISE_IF(t->nvar()!=nv, -"Wrong number of variables in TensorPolynomial::insert"); -TensorContainer<_Ttype> ::insert(t); -if(maxdim<t->dimen()) -maxdim= t->dimen(); -} - -/*:10*/ -; -/*11:*/ - -void derivative(int k) -{ -for(int d= 1;d<=maxdim;d++){ -if(_Tparent::check(Symmetry(d))){ -_Ttype*ten= _Tparent::get(Symmetry(d)); -ten->mult((double)max((d-k),0)); -} -} -} - -/*:11*/ -; -/*12:*/ - -_Ttype*evalPartially(int s,const ConstVector&v) -{ -TL_RAISE_IF(v.length()!=nvars(), -"Wrong length of vector for TensorPolynomial::evalPartially"); - -_Ttype*res= new _Ttype(nrows(),nvars(),s); -res->zeros(); - -int sfac= 1; -for(int i= 1;i<=s;i++) -sfac*= i; - -if(_Tparent::check(Symmetry(s))) -res->add(1.0/sfac,*(_Tparent::get(Symmetry(s)))); - -int dfac= sfac*(s+1); -for(int d= s+1;d<=maxdim;d++,dfac*= d){ -if(_Tparent::check(Symmetry(d))){ -const _Ttype<mp= *(_Tparent::get(Symmetry(d))); -_Ttype*last= new _Ttype(ltmp); -for(int j= 0;j<d-s;j++){ -_Ttype*newlast= new _Ttype(*last,v); -delete last; -last= newlast; -} -res->add(1.0/dfac,*last); -delete last; -} -} - -return res; -} - -/*:12*/ -; + TensorPolynomial(int rows,int vars) + :TensorContainer<_Ttype> (1), + nr(rows),nv(vars),maxdim(0){} + TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,int k) + :TensorContainer<_Ttype> (tp), + nr(tp.nr),nv(tp.nv),maxdim(0){derivative(k);} + TensorPolynomial(int first_row,int num,TensorPolynomial<_Ttype,_TGStype,_Stype> &tp) + :TensorContainer<_Ttype> (first_row,num,tp), + nr(num),nv(tp.nv),maxdim(tp.maxdim){} + /*4:*/ + + TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,const Vector&xval) + :TensorContainer<_Ttype> (1), + nr(tp.nrows()),nv(tp.nvars()-xval.length()),maxdim(0) + { + TL_RAISE_IF(nvars()<0, + "Length of xval too big in TensorPolynomial contract constructor"); + IntSequence ss(2);ss[0]= xval.length();ss[1]= nvars(); + IntSequence pp(2);pp[0]= 0;pp[1]= 1; + + /*5:*/ + + PowerProvider pwp(xval); + for(int i= 1;i<=tp.maxdim;i++){ + const _Stype&xpow= pwp.getNext((const _Stype*)NULL); + for(int j= 0;j<=tp.maxdim-i;j++){ + if(tp.check(Symmetry(i+j))){ + /*7:*/ + + _Ttype*ten; + if(_Tparent::check(Symmetry(j))){ + ten= _Tparent::get(Symmetry(j)); + }else{ + ten= new _Ttype(nrows(),nvars(),j); + ten->zeros(); + insert(ten); + } + + + /*:7*/ + ; + Symmetry sym(i,j); + IntSequence coor(sym,pp); + _TGStype slice(*(tp.get(Symmetry(i+j))),ss,coor,TensorDimens(sym,ss)); + slice.mult(Tensor::noverk(i+j,j)); + _TGStype tmp(*ten); + slice.contractAndAdd(0,tmp,xpow); + } + } + } + + /*:5*/ + ; + /*6:*/ + + for(int j= 0;j<=tp.maxdim;j++){ + if(tp.check(Symmetry(j))){ + /*7:*/ + + _Ttype*ten; + if(_Tparent::check(Symmetry(j))){ + ten= _Tparent::get(Symmetry(j)); + }else{ + ten= new _Ttype(nrows(),nvars(),j); + ten->zeros(); + insert(ten); + } + + + /*:7*/ + ; + Symmetry sym(0,j); + IntSequence coor(sym,pp); + _TGStype slice(*(tp.get(Symmetry(j))),ss,coor,TensorDimens(sym,ss)); + ten->add(1.0,slice); + } + } + + + /*:6*/ + ; + } + + /*:4*/ + ; + TensorPolynomial(const TensorPolynomial&tp) + :TensorContainer<_Ttype> (tp),nr(tp.nr),nv(tp.nv),maxdim(tp.maxdim){} + int nrows()const + {return nr;} + int nvars()const + {return nv;} + /*8:*/ + + void evalTrad(Vector&out,const ConstVector&v)const + { + if(_Tparent::check(Symmetry(0))) + out= _Tparent::get(Symmetry(0))->getData(); + else + out.zeros(); + + PowerProvider pp(v); + for(int d= 1;d<=maxdim;d++){ + const _Stype&p= pp.getNext((const _Stype*)NULL); + Symmetry cs(d); + if(_Tparent::check(cs)){ + const _Ttype*t= _Tparent::get(cs); + t->multaVec(out,p.getData()); + } + } + } + + /*:8*/ + ; + /*9:*/ + + void evalHorner(Vector&out,const ConstVector&v)const + { + if(_Tparent::check(Symmetry(0))) + out= _Tparent::get(Symmetry(0))->getData(); + else + out.zeros(); + + if(maxdim==0) + return; + + _Ttype*last; + if(maxdim==1) + last= new _Ttype(*(_Tparent::get(Symmetry(1)))); + else + last= new _Ttype(*(_Tparent::get(Symmetry(maxdim))),v); + for(int d= maxdim-1;d>=1;d--){ + Symmetry cs(d); + if(_Tparent::check(cs)){ + const _Ttype*nt= _Tparent::get(cs); + last->add(1.0,ConstTwoDMatrix(*nt)); + } + if(d> 1){ + _Ttype*new_last= new _Ttype(*last,v); + delete last; + last= new_last; + } + } + last->multaVec(out,v); + delete last; + } + + /*:9*/ + ; + /*10:*/ + + void insert(_ptr t) + { + TL_RAISE_IF(t->nrows()!=nr, + "Wrong number of rows in TensorPolynomial::insert"); + TL_RAISE_IF(t->nvar()!=nv, + "Wrong number of variables in TensorPolynomial::insert"); + TensorContainer<_Ttype> ::insert(t); + if(maxdim<t->dimen()) + maxdim= t->dimen(); + } + + /*:10*/ + ; + /*11:*/ + + void derivative(int k) + { + for(int d= 1;d<=maxdim;d++){ + if(_Tparent::check(Symmetry(d))){ + _Ttype*ten= _Tparent::get(Symmetry(d)); + ten->mult((double)max((d-k),0)); + } + } + } + + /*:11*/ + ; + /*12:*/ + + _Ttype*evalPartially(int s,const ConstVector&v) + { + TL_RAISE_IF(v.length()!=nvars(), + "Wrong length of vector for TensorPolynomial::evalPartially"); + + _Ttype*res= new _Ttype(nrows(),nvars(),s); + res->zeros(); + + int sfac= 1; + for(int i= 1;i<=s;i++) + sfac*= i; + + if(_Tparent::check(Symmetry(s))) + res->add(1.0/sfac,*(_Tparent::get(Symmetry(s)))); + + int dfac= sfac*(s+1); + for(int d= s+1;d<=maxdim;d++,dfac*= d){ + if(_Tparent::check(Symmetry(d))){ + const _Ttype<mp= *(_Tparent::get(Symmetry(d))); + _Ttype*last= new _Ttype(ltmp); + for(int j= 0;j<d-s;j++){ + _Ttype*newlast= new _Ttype(*last,v); + delete last; + last= newlast; + } + res->add(1.0/dfac,*last); + delete last; + } + } + + return res; + } + + /*:12*/ + ; }; @@ -256,15 +256,15 @@ return res; class FTensorPolynomial; class UTensorPolynomial:public TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> { public: -UTensorPolynomial(int rows,int vars) -:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (rows,vars){} -UTensorPolynomial(const UTensorPolynomial&up,int k) -:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (up,k){} -UTensorPolynomial(const FTensorPolynomial&fp); -UTensorPolynomial(const UTensorPolynomial&tp,const Vector&xval) -:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (tp,xval){} -UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp) -:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (first_row,num,tp){} + UTensorPolynomial(int rows,int vars) + :TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (rows,vars){} + UTensorPolynomial(const UTensorPolynomial&up,int k) + :TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (up,k){} + UTensorPolynomial(const FTensorPolynomial&fp); + UTensorPolynomial(const UTensorPolynomial&tp,const Vector&xval) + :TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (tp,xval){} + UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp) + :TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (first_row,num,tp){} }; /*:13*/ @@ -273,15 +273,15 @@ UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp) class FTensorPolynomial:public TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> { public: -FTensorPolynomial(int rows,int vars) -:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (rows,vars){} -FTensorPolynomial(const FTensorPolynomial&fp,int k) -:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fp,k){} -FTensorPolynomial(const UTensorPolynomial&up); -FTensorPolynomial(const FTensorPolynomial&tp,const Vector&xval) -:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (tp,xval){} -FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp) -:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (first_row,num,tp){} + FTensorPolynomial(int rows,int vars) + :TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (rows,vars){} + FTensorPolynomial(const FTensorPolynomial&fp,int k) + :TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fp,k){} + FTensorPolynomial(const UTensorPolynomial&up); + FTensorPolynomial(const FTensorPolynomial&tp,const Vector&xval) + :TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (tp,xval){} + FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp) + :TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (first_row,num,tp){} }; /*:14*/ @@ -291,61 +291,61 @@ FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp) template<class _Ttype,class _TGStype,class _Stype> class CompactPolynomial:public _Ttype{ public: -/*16:*/ - -CompactPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &pol) -:_Ttype(pol.nrows(),pol.nvars()+1,pol.getMaxDim()) -{ -_Ttype::zeros(); - -IntSequence dumnvs(2); -dumnvs[0]= 1; -dumnvs[1]= pol.nvars(); - -int offset= 0; -_Ttype dum(0,2,_Ttype::dimen()); -for(Tensor::index i= dum.begin();i!=dum.end();++i){ -int d= i.getCoor().sum(); -Symmetry symrun(_Ttype::dimen()-d,d); -_TGStype dumgs(0,TensorDimens(symrun,dumnvs)); -if(pol.check(Symmetry(d))){ -TwoDMatrix subt(*this,offset,dumgs.ncols()); -subt.add(1.0,*(pol.get(Symmetry(d)))); -} -offset+= dumgs.ncols(); -} -} - - -/*:16*/ -; -/*17:*/ - -void eval(Vector&out,const ConstVector&v)const -{ -TL_RAISE_IF(v.length()+1!=_Ttype::nvar(), -"Wrong input vector length in CompactPolynomial::eval"); -TL_RAISE_IF(out.length()!=_Ttype::nrows(), -"Wrong output vector length in CompactPolynomial::eval"); - -Vector x1(v.length()+1); -Vector x1p(x1,1,v.length()); -x1p= v; -x1[0]= 1.0; - -if(_Ttype::dimen()==0) -out= ConstVector(*this,0); -else{ -PowerProvider pp(x1); -const _Stype&xpow= pp.getNext((const _Stype*)NULL); -for(int i= 1;i<_Ttype::dimen();i++) -xpow= pp.getNext((const _Stype*)NULL); -multVec(0.0,out,1.0,xpow); -} -} - -/*:17*/ -; + /*16:*/ + + CompactPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &pol) + :_Ttype(pol.nrows(),pol.nvars()+1,pol.getMaxDim()) + { + _Ttype::zeros(); + + IntSequence dumnvs(2); + dumnvs[0]= 1; + dumnvs[1]= pol.nvars(); + + int offset= 0; + _Ttype dum(0,2,_Ttype::dimen()); + for(Tensor::index i= dum.begin();i!=dum.end();++i){ + int d= i.getCoor().sum(); + Symmetry symrun(_Ttype::dimen()-d,d); + _TGStype dumgs(0,TensorDimens(symrun,dumnvs)); + if(pol.check(Symmetry(d))){ + TwoDMatrix subt(*this,offset,dumgs.ncols()); + subt.add(1.0,*(pol.get(Symmetry(d)))); + } + offset+= dumgs.ncols(); + } + } + + + /*:16*/ + ; + /*17:*/ + + void eval(Vector&out,const ConstVector&v)const + { + TL_RAISE_IF(v.length()+1!=_Ttype::nvar(), + "Wrong input vector length in CompactPolynomial::eval"); + TL_RAISE_IF(out.length()!=_Ttype::nrows(), + "Wrong output vector length in CompactPolynomial::eval"); + + Vector x1(v.length()+1); + Vector x1p(x1,1,v.length()); + x1p= v; + x1[0]= 1.0; + + if(_Ttype::dimen()==0) + out= ConstVector(*this,0); + else{ + PowerProvider pp(x1); + const _Stype&xpow= pp.getNext((const _Stype*)NULL); + for(int i= 1;i<_Ttype::dimen();i++) + xpow= pp.getNext((const _Stype*)NULL); + multVec(0.0,out,1.0,xpow); + } + } + + /*:17*/ + ; }; /*:15*/ @@ -354,8 +354,8 @@ multVec(0.0,out,1.0,xpow); class UCompactPolynomial:public CompactPolynomial<UFSTensor,UGSTensor,URSingleTensor> { public: -UCompactPolynomial(const UTensorPolynomial&upol) -:CompactPolynomial<UFSTensor,UGSTensor,URSingleTensor> (upol){} + UCompactPolynomial(const UTensorPolynomial&upol) + :CompactPolynomial<UFSTensor,UGSTensor,URSingleTensor> (upol){} }; /*:18*/ @@ -364,8 +364,8 @@ UCompactPolynomial(const UTensorPolynomial&upol) class FCompactPolynomial:public CompactPolynomial<FFSTensor,FGSTensor,FRSingleTensor> { public: -FCompactPolynomial(const FTensorPolynomial&fpol) -:CompactPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fpol){} + FCompactPolynomial(const FTensorPolynomial&fpol) + :CompactPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fpol){} }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/tensor.h b/mex/sources/korderpert/Dyn_pp/tl/cc/tensor.h index 44acdba41e29745c09247a937519827dfd4e4031..7e8fb0164a49904871ddb7690d3f3be8f426aa42 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/tensor.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/tensor.h @@ -10,36 +10,36 @@ /*2:*/ template<class _Tptr> class _index{ -typedef _index<_Tptr> _Self; -_Tptr tensor; -int offset; -IntSequence coor; + typedef _index<_Tptr> _Self; + _Tptr tensor; + int offset; + IntSequence coor; public: -_index(_Tptr t,int n) -:tensor(t),offset(0),coor(n,0){} -_index(_Tptr t,const IntSequence&cr,int c) -:tensor(t),offset(c),coor(cr){} -_index(_Tptr t,const IntSequence&cr) -:tensor(t),offset(tensor->getOffset(cr)),coor(cr){} -_index(const _index&ind) -:tensor(ind.tensor),offset(ind.offset),coor(ind.coor){} -const _Self&operator= (const _Self&in) -{tensor= in.tensor;offset= in.offset;coor= in.coor; -return*this;} -_Self&operator++() -{tensor->increment(coor);offset++;return*this;} -_Self&operator--() -{tensor->decrement(coor);offset--;return*this;} -int operator*()const -{return offset;} -bool operator==(const _index&n)const -{return offset==n.offset;} -bool operator!=(const _index&n)const -{return offset!=n.offset;} -const IntSequence&getCoor()const -{return coor;} -void print()const -{printf("%4d: ",offset);coor.print();} + _index(_Tptr t,int n) + :tensor(t),offset(0),coor(n,0){} + _index(_Tptr t,const IntSequence&cr,int c) + :tensor(t),offset(c),coor(cr){} + _index(_Tptr t,const IntSequence&cr) + :tensor(t),offset(tensor->getOffset(cr)),coor(cr){} + _index(const _index&ind) + :tensor(ind.tensor),offset(ind.offset),coor(ind.coor){} + const _Self&operator= (const _Self&in) + {tensor= in.tensor;offset= in.offset;coor= in.coor; + return*this;} + _Self&operator++() + {tensor->increment(coor);offset++;return*this;} + _Self&operator--() + {tensor->decrement(coor);offset--;return*this;} + int operator*()const + {return offset;} + bool operator==(const _index&n)const + {return offset==n.offset;} + bool operator!=(const _index&n)const + {return offset!=n.offset;} + const IntSequence&getCoor()const + {return coor;} + void print()const + {printf("%4d: ",offset);coor.print();} }; /*:2*/ @@ -48,55 +48,55 @@ void print()const class Tensor:public TwoDMatrix{ public: -enum indor{along_row,along_col}; -typedef _index<const Tensor*> index; + enum indor{along_row,along_col}; + typedef _index<const Tensor*> index; protected: -const index in_beg; -const index in_end; -int dim; + const index in_beg; + const index in_end; + int dim; public: -Tensor(indor io,const IntSequence&last,int r,int c,int d) -:TwoDMatrix(r,c), -in_beg(this,d), -in_end(this,last,(io==along_row)?r:c), -dim(d){} -Tensor(indor io,const IntSequence&first,const IntSequence&last, -int r,int c,int d) -:TwoDMatrix(r,c), -in_beg(this,first,0), -in_end(this,last,(io==along_row)?r:c), -dim(d){} -Tensor(int first_row,int num,Tensor&t) -:TwoDMatrix(first_row,num,t), -in_beg(t.in_beg), -in_end(t.in_end), -dim(t.dim){} -Tensor(const Tensor&t) -:TwoDMatrix(t), -in_beg(this,t.in_beg.getCoor(),*(t.in_beg)), -in_end(this,t.in_end.getCoor(),*(t.in_end)), -dim(t.dim){} -virtual~Tensor(){} -virtual void increment(IntSequence&v)const= 0; -virtual void decrement(IntSequence&v)const= 0; -virtual int getOffset(const IntSequence&v)const= 0; -int dimen()const -{return dim;} - -const index&begin()const -{return in_beg;} -const index&end()const -{return in_end;} - -static int noverk(int n,int k); -static int power(int a,int b); -static int noverseq(const IntSequence&s) -{ -IntSequence seq(s); -return noverseq_ip((IntSequence&)s); -} + Tensor(indor io,const IntSequence&last,int r,int c,int d) + :TwoDMatrix(r,c), + in_beg(this,d), + in_end(this,last,(io==along_row)?r:c), + dim(d){} + Tensor(indor io,const IntSequence&first,const IntSequence&last, + int r,int c,int d) + :TwoDMatrix(r,c), + in_beg(this,first,0), + in_end(this,last,(io==along_row)?r:c), + dim(d){} + Tensor(int first_row,int num,Tensor&t) + :TwoDMatrix(first_row,num,t), + in_beg(t.in_beg), + in_end(t.in_end), + dim(t.dim){} + Tensor(const Tensor&t) + :TwoDMatrix(t), + in_beg(this,t.in_beg.getCoor(),*(t.in_beg)), + in_end(this,t.in_end.getCoor(),*(t.in_end)), + dim(t.dim){} + virtual~Tensor(){} + virtual void increment(IntSequence&v)const= 0; + virtual void decrement(IntSequence&v)const= 0; + virtual int getOffset(const IntSequence&v)const= 0; + int dimen()const + {return dim;} + + const index&begin()const + {return in_beg;} + const index&end()const + {return in_end;} + + static int noverk(int n,int k); + static int power(int a,int b); + static int noverseq(const IntSequence&s) + { + IntSequence seq(s); + return noverseq_ip((IntSequence&)s); + } private: -static int noverseq_ip(IntSequence&s); + static int noverseq_ip(IntSequence&s); }; /*:3*/ @@ -106,21 +106,21 @@ static int noverseq_ip(IntSequence&s); class FTensor; class UTensor:public Tensor{ public: -UTensor(indor io,const IntSequence&last,int r,int c,int d) -:Tensor(io,last,r,c,d){} -UTensor(const UTensor&ut) -:Tensor(ut){} -UTensor(int first_row,int num,UTensor&t) -:Tensor(first_row,num,t){} -virtual~UTensor(){} -virtual FTensor&fold()const= 0; - -static void increment(IntSequence&v,int nv); -static void decrement(IntSequence&v,int nv); -static void increment(IntSequence&v,const IntSequence&nvmx); -static void decrement(IntSequence&v,const IntSequence&nvmx); -static int getOffset(const IntSequence&v,int nv); -static int getOffset(const IntSequence&v,const IntSequence&nvmx); + UTensor(indor io,const IntSequence&last,int r,int c,int d) + :Tensor(io,last,r,c,d){} + UTensor(const UTensor&ut) + :Tensor(ut){} + UTensor(int first_row,int num,UTensor&t) + :Tensor(first_row,num,t){} + virtual~UTensor(){} + virtual FTensor&fold()const= 0; + + static void increment(IntSequence&v,int nv); + static void decrement(IntSequence&v,int nv); + static void increment(IntSequence&v,const IntSequence&nvmx); + static void decrement(IntSequence&v,const IntSequence&nvmx); + static int getOffset(const IntSequence&v,int nv); + static int getOffset(const IntSequence&v,const IntSequence&nvmx); }; /*:4*/ @@ -129,20 +129,20 @@ static int getOffset(const IntSequence&v,const IntSequence&nvmx); class FTensor:public Tensor{ public: -FTensor(indor io,const IntSequence&last,int r,int c,int d) -:Tensor(io,last,r,c,d){} -FTensor(const FTensor&ft) -:Tensor(ft){} -FTensor(int first_row,int num,FTensor&t) -:Tensor(first_row,num,t){} -virtual~FTensor(){} -virtual UTensor&unfold()const= 0; - -static void decrement(IntSequence&v,int nv); -static int getOffset(const IntSequence&v,int nv) -{IntSequence vtmp(v);return getOffsetRecurse(vtmp,nv);} + FTensor(indor io,const IntSequence&last,int r,int c,int d) + :Tensor(io,last,r,c,d){} + FTensor(const FTensor&ft) + :Tensor(ft){} + FTensor(int first_row,int num,FTensor&t) + :Tensor(first_row,num,t){} + virtual~FTensor(){} + virtual UTensor&unfold()const= 0; + + static void decrement(IntSequence&v,int nv); + static int getOffset(const IntSequence&v,int nv) + {IntSequence vtmp(v);return getOffsetRecurse(vtmp,nv);} private: -static int getOffsetRecurse(IntSequence&v,int nv); + static int getOffsetRecurse(IntSequence&v,int nv); }; /*:5*/ diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/tl_exception.h b/mex/sources/korderpert/Dyn_pp/tl/cc/tl_exception.h index 2fd08885fef49485f47447056e6b8d7abb4c77e0..71d4cc8a4f2f6368034894fc8ea271ad62ce506c 100644 --- a/mex/sources/korderpert/Dyn_pp/tl/cc/tl_exception.h +++ b/mex/sources/korderpert/Dyn_pp/tl/cc/tl_exception.h @@ -25,19 +25,19 @@ if (TL_DEBUG >= TL_DEBUG_EXCEPTION && (expr)) throw TLException(__FILE__, __LINE /*3:*/ class TLException{ -char fname[50]; -int lnum; -char message[500]; + char fname[50]; + int lnum; + char message[500]; public: -TLException(const char*f,int l,const char*mes) -{ -strncpy(fname,f,50);fname[49]= '\0'; -strncpy(message,mes,500);message[499]= '\0'; -lnum= l; -} -virtual~TLException(){} -virtual void print()const -{printf("At %s:%d:%s\n",fname,lnum,message);} + TLException(const char*f,int l,const char*mes) + { + strncpy(fname,f,50);fname[49]= '\0'; + strncpy(message,mes,500);message[499]= '\0'; + lnum= l; + } + virtual~TLException(){} + virtual void print()const + {printf("At %s:%d:%s\n",fname,lnum,message);} };