From 0337bae230a884326f4c0bb6e8182c74310dab81 Mon Sep 17 00:00:00 2001 From: george Date: Mon, 20 Oct 2008 16:19:58 +0000 Subject: [PATCH] More Dynare_pp reformatted and cleared from CWEB references git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2183 ac1d8469-bf42-47a9-8791-bf33cf982152 --- .../korderpert/Dyn_pp/tl/cc/normal_moments.h | 8 +- .../korderpert/Dyn_pp/tl/cc/permutation.h | 96 +- .../korderpert/Dyn_pp/tl/cc/ps_tensor.h | 266 ++--- .../korderpert/Dyn_pp/tl/cc/pyramid_prod.h | 8 +- .../korderpert/Dyn_pp/tl/cc/pyramid_prod2.h | 42 +- .../korderpert/Dyn_pp/tl/cc/rfs_tensor.h | 122 +-- .../korderpert/Dyn_pp/tl/cc/sparse_tensor.h | 130 +-- .../korderpert/Dyn_pp/tl/cc/stack_container.h | 934 +++++++++--------- mex/sources/korderpert/Dyn_pp/tl/cc/sthread.h | 896 ++++++++--------- .../korderpert/Dyn_pp/tl/cc/symmetry.h | 134 +-- .../korderpert/Dyn_pp/tl/cc/t_container.h | 468 ++++----- .../korderpert/Dyn_pp/tl/cc/t_polynomial.h | 610 ++++++------ mex/sources/korderpert/Dyn_pp/tl/cc/tensor.h | 206 ++-- .../korderpert/Dyn_pp/tl/cc/tl_exception.h | 24 +- 14 files changed, 1972 insertions(+), 1972 deletions(-) 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 f179f78a3..8a2cf9194 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 { 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 { 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 0ce3345be..3388aa2a2 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 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 getPreserving(const IntSequence&s)const; }; @@ -74,12 +74,12 @@ vector getPreserving(const IntSequence&s)const; /*4:*/ class PermutationBundle{ -vector bundle; + vector 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 b2df33b39..081d6e19c 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 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 f3a519273..795e30c11 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 &lst); -void addKronColumn(int i,const vector &ts, -const IntSequence&pindex); + USubTensor(const TensorDimens&bdims,const TensorDimens&hdims, + const FGSContainer&cont,const vector &lst); + void addKronColumn(int i,const vector &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 ee0b80657..266056931 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 &sp,const IntSequence&c); -~IrregTensorHeader(); -int dimen()const -{return unit_flag.size();} -void increment(IntSequence&v)const; -int calcMaxOffset()const; + IrregTensorHeader(const StackProduct &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 10fbf2479..00b3072bb 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 &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 &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 66976d7ce..110ee3950 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 Item; -typedef multimap Map; -typedef Map::const_iterator const_iterator; + typedef pair Item; + typedef multimap 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 55fb96039..e9bc5ad67 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 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 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=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(), + "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=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 { -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 &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 &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 { -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 &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 &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 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 , public FoldedStackContainer{ public: -typedef TensorContainer _Ctype; -FoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, -int ny,int nu) -:ZContainer (gss,ngss,g,ng,ny,nu){} + typedef TensorContainer _Ctype; + FoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, + int ny,int nu) + :ZContainer (gss,ngss,g,ng,ny,nu){} }; /*:14*/ @@ -307,10 +307,10 @@ int ny,int nu) class UnfoldedZContainer:public ZContainer , public UnfoldedStackContainer{ public: -typedef TensorContainer _Ctype; -UnfoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, -int ny,int nu) -:ZContainer (gss,ngss,g,ng,ny,nu){} + typedef TensorContainer _Ctype; + UnfoldedZContainer(const _Ctype*gss,int ngss,const _Ctype*g,int ng, + int ny,int nu) + :ZContainer (gss,ngss,g,ng,ny,nu){} }; /*:15*/ @@ -320,48 +320,48 @@ int ny,int nu) template 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 , public FoldedStackContainer{ public: -typedef TensorContainer _Ctype; -FoldedGContainer(const _Ctype*gs,int ngs,int nu) -:GContainer (gs,ngs,nu){} + typedef TensorContainer _Ctype; + FoldedGContainer(const _Ctype*gs,int ngs,int nu) + :GContainer (gs,ngs,nu){} }; /*:18*/ @@ -383,9 +383,9 @@ FoldedGContainer(const _Ctype*gs,int ngs,int nu) class UnfoldedGContainer:public GContainer , public UnfoldedStackContainer{ public: -typedef TensorContainer _Ctype; -UnfoldedGContainer(const _Ctype*gs,int ngs,int nu) -:GContainer (gs,ngs,nu){} + typedef TensorContainer _Ctype; + UnfoldedGContainer(const _Ctype*gs,int ngs,int nu) + :GContainer (gs,ngs,nu){} }; @@ -396,112 +396,112 @@ UnfoldedGContainer(const _Ctype*gs,int ngs,int nu) template 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=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=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 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;inrows()!=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;inrows()!=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 2e82be645..27964291e 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 namespace sthread{ -using namespace std; - -class Empty{}; -/*2:*/ - -template -struct IF{ -typedef Then RET; -}; - -template -struct IF { -typedef Else RET; -}; - - - -/*:2*/ -; -enum{posix,empty}; -template class thread_traits; -template class detach_thread; -/*3:*/ - -template -class thread{ -typedef thread_traits _Ttraits; -typedef typename _Ttraits::_Tthread _Tthread; -_Tthread th; + using namespace std; + + class Empty{}; + /*2:*/ + + template + struct IF{ + typedef Then RET; + }; + + template + struct IF { + typedef Else RET; + }; + + + + /*:2*/ + ; + enum{posix,empty}; + template class thread_traits; + template class detach_thread; + /*3:*/ + + template + class thread{ + typedef thread_traits _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 + class thread_group{ + typedef thread_traits _Ttraits; + typedef thread _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;crun(); + } + iterator ret; + c= 0; + for(ret= start;c + struct thread_traits{ + typedef typename IF ::RET _Tthread; + typedef thread _Ctype; + typedef detach_thread _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 mmkey; + + template + struct mutex_traits{ + typedef typename IF ::RET _Tmutex; + typedef map ,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 + class mutex_map + :public mutex_traits ::mutex_int_map + { + typedef typename mutex_traits ::_Tmutex _Tmutex; + typedef mutex_traits _Mtraits; + typedef pair<_Tmutex,int> mmval; + typedef map _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 + class synchro{ + typedef typename mutex_traits ::_Tmutex _Tmutex; + typedef mutex_traits _Mtraits; + public: + typedef mutex_map 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 + struct cond_traits{ + typedef typename IF ::RET _Tcond; + typedef typename mutex_traits ::_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 + class condition_counter{ + typedef typename mutex_traits ::_Tmutex _Tmutex; + typedef typename cond_traits ::_Tcond _Tcond; + int counter; + _Tmutex mut; + _Tcond cond; + bool changed; 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 -class thread_group{ -typedef thread_traits _Ttraits; -typedef thread _Ctype; -list<_Ctype*> tlist; -typedef typename list<_Ctype*> ::iterator iterator; + /*18:*/ + + condition_counter() + :counter(0),changed(true) + { + mutex_traits ::init(mut); + cond_traits ::init(cond); + } + + /*:18*/ + ; + /*19:*/ + + ~condition_counter() + { + cond_traits ::destroy(cond); + } + + /*:19*/ + ; + /*20:*/ + + void increase() + { + mutex_traits ::lock(mut); + counter++; + changed= true; + cond_traits ::broadcast(cond); + mutex_traits ::unlock(mut); + } + + /*:20*/ + ; + /*21:*/ + + void decrease() + { + mutex_traits ::lock(mut); + counter--; + changed= true; + cond_traits ::broadcast(cond); + mutex_traits ::unlock(mut); + } + + /*:21*/ + ; + /*22:*/ + + int waitForChange() + { + mutex_traits ::lock(mut); + if(!changed){ + cond_traits ::wait(cond,mut); + } + changed= false; + int res= counter; + mutex_traits ::unlock(mut); + return res; + } + + + /*:22*/ + ; + }; + + /*:17*/ + ; + /*23:*/ + + template + class detach_thread:public thread { 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;crun(); -} -iterator ret; -c= 0; -for(ret= start;c -struct thread_traits{ -typedef typename IF ::RET _Tthread; -typedef thread _Ctype; -typedef detach_thread _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 mmkey; - -template -struct mutex_traits{ -typedef typename IF ::RET _Tmutex; -typedef map ,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 -class mutex_map -:public mutex_traits ::mutex_int_map -{ -typedef typename mutex_traits ::_Tmutex _Tmutex; -typedef mutex_traits _Mtraits; -typedef pair<_Tmutex,int> mmval; -typedef map _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 -class synchro{ -typedef typename mutex_traits ::_Tmutex _Tmutex; -typedef mutex_traits _Mtraits; -public: -typedef mutex_map 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 -struct cond_traits{ -typedef typename IF ::RET _Tcond; -typedef typename mutex_traits ::_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 -class condition_counter{ -typedef typename mutex_traits ::_Tmutex _Tmutex; -typedef typename cond_traits ::_Tcond _Tcond; -int counter; -_Tmutex mut; -_Tcond cond; -bool changed; -public: -/*18:*/ - -condition_counter() -:counter(0),changed(true) -{ -mutex_traits ::init(mut); -cond_traits ::init(cond); -} - -/*:18*/ -; -/*19:*/ - -~condition_counter() -{ -cond_traits ::destroy(cond); -} - -/*:19*/ -; -/*20:*/ - -void increase() -{ -mutex_traits ::lock(mut); -counter++; -changed= true; -cond_traits ::broadcast(cond); -mutex_traits ::unlock(mut); -} - -/*:20*/ -; -/*21:*/ - -void decrease() -{ -mutex_traits ::lock(mut); -counter--; -changed= true; -cond_traits ::broadcast(cond); -mutex_traits ::unlock(mut); -} - -/*:21*/ -; -/*22:*/ - -int waitForChange() -{ -mutex_traits ::lock(mut); -if(!changed){ -cond_traits ::wait(cond,mut); -} -changed= false; -int res= counter; -mutex_traits ::unlock(mut); -return res; -} - - -/*:22*/ -; -}; - -/*:17*/ -; -/*23:*/ - -template -class detach_thread:public thread { -public: -condition_counter *counter; -detach_thread():counter(NULL){} -void installCounter(condition_counter *c) -{counter= c;} -void run() -{thread_traits ::detach_run(this);} -}; - -/*:23*/ -; -/*24:*/ - -template -class detach_thread_group{ -typedef thread_traits _Ttraits; -typedef cond_traits _Ctraits; -typedef detach_thread _Ctype; -list<_Ctype*> tlist; -typedef typename list<_Ctype*> ::iterator iterator; -condition_counter 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()run(); -++it; -} -} -while(counter.waitForChange()> 0){} -} - - -/*:27*/ -; -}; - -/*:24*/ -; + condition_counter *counter; + detach_thread():counter(NULL){} + void installCounter(condition_counter *c) + {counter= c;} + void run() + {thread_traits ::detach_run(this);} + }; + + /*:23*/ + ; + /*24:*/ + + template + class detach_thread_group{ + typedef thread_traits _Ttraits; + typedef cond_traits _Ctraits; + typedef detach_thread _Ctype; + list<_Ctype*> tlist; + typedef typename list<_Ctype*> ::iterator iterator; + condition_counter 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()run(); + ++it; + } + } + while(counter.waitForChange()> 0){} + } + + + /*:27*/ + ; + }; + + /*:24*/ + ; #ifdef POSIX_THREADS -/*28:*/ - -typedef detach_thread PosixThread; -typedef detach_thread_group PosixThreadGroup; -typedef synchro posix_synchro; -class PosixSynchro:public posix_synchro{ -public: -PosixSynchro(const void*c,const char*id); -}; - + /*28:*/ + + typedef detach_thread PosixThread; + typedef detach_thread_group PosixThreadGroup; + typedef synchro 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 NoThread; -typedef thread_group NoThreadGroup; -typedef synchro no_synchro; -class NoSynchro{ -public: -NoSynchro(const void*c,const char*id){} -~NoSynchro(){} -}; - + /*29:*/ + + typedef thread NoThread; + typedef thread_group NoThreadGroup; + typedef synchro 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 8700b33f4..3a5c8447c 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 { 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 392c27be1..3ec7f86d9 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 class TensorContainer{ protected: -typedef const _Ttype*_const_ptr; -typedef _Ttype*_ptr; -typedef map _Map; -typedef typename _Map::value_type _mvtype; + typedef const _Ttype*_const_ptr; + typedef _Ttype*_ptr; + typedef map _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 &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 { public: -UGSContainer(int nn) -:TensorContainer (nn){} -UGSContainer(const UGSContainer&uc) -:TensorContainer (uc){} -UGSContainer(const FGSContainer&c); -void multAndAdd(const UGSTensor&t,UGSTensor&out)const; + UGSContainer(int nn) + :TensorContainer (nn){} + UGSContainer(const UGSContainer&uc) + :TensorContainer (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 { -static const int num_one_time; + static const int num_one_time; public: -FGSContainer(int nn) -:TensorContainer (nn){} -FGSContainer(const FGSContainer&fc) -:TensorContainer (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 (nn){} + FGSContainer(const FGSContainer&fc) + :TensorContainer (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 &out, -const Tensor::index&start, -const Tensor::index&end); + static Tensor::index + getIndices(int num,vector &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 535ad7818..84e8739df 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 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(maxdimdimen()) -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;jadd(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(maxdimdimen()) + 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;jadd(1.0/dfac,*last); + delete last; + } + } + + return res; + } + + /*:12*/ + ; }; @@ -256,15 +256,15 @@ return res; class FTensorPolynomial; class UTensorPolynomial:public TensorPolynomial { public: -UTensorPolynomial(int rows,int vars) -:TensorPolynomial (rows,vars){} -UTensorPolynomial(const UTensorPolynomial&up,int k) -:TensorPolynomial (up,k){} -UTensorPolynomial(const FTensorPolynomial&fp); -UTensorPolynomial(const UTensorPolynomial&tp,const Vector&xval) -:TensorPolynomial (tp,xval){} -UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp) -:TensorPolynomial (first_row,num,tp){} + UTensorPolynomial(int rows,int vars) + :TensorPolynomial (rows,vars){} + UTensorPolynomial(const UTensorPolynomial&up,int k) + :TensorPolynomial (up,k){} + UTensorPolynomial(const FTensorPolynomial&fp); + UTensorPolynomial(const UTensorPolynomial&tp,const Vector&xval) + :TensorPolynomial (tp,xval){} + UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp) + :TensorPolynomial (first_row,num,tp){} }; /*:13*/ @@ -273,15 +273,15 @@ UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp) class FTensorPolynomial:public TensorPolynomial { public: -FTensorPolynomial(int rows,int vars) -:TensorPolynomial (rows,vars){} -FTensorPolynomial(const FTensorPolynomial&fp,int k) -:TensorPolynomial (fp,k){} -FTensorPolynomial(const UTensorPolynomial&up); -FTensorPolynomial(const FTensorPolynomial&tp,const Vector&xval) -:TensorPolynomial (tp,xval){} -FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp) -:TensorPolynomial (first_row,num,tp){} + FTensorPolynomial(int rows,int vars) + :TensorPolynomial (rows,vars){} + FTensorPolynomial(const FTensorPolynomial&fp,int k) + :TensorPolynomial (fp,k){} + FTensorPolynomial(const UTensorPolynomial&up); + FTensorPolynomial(const FTensorPolynomial&tp,const Vector&xval) + :TensorPolynomial (tp,xval){} + FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp) + :TensorPolynomial (first_row,num,tp){} }; /*:14*/ @@ -291,61 +291,61 @@ FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp) template 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 { public: -UCompactPolynomial(const UTensorPolynomial&upol) -:CompactPolynomial (upol){} + UCompactPolynomial(const UTensorPolynomial&upol) + :CompactPolynomial (upol){} }; /*:18*/ @@ -364,8 +364,8 @@ UCompactPolynomial(const UTensorPolynomial&upol) class FCompactPolynomial:public CompactPolynomial { public: -FCompactPolynomial(const FTensorPolynomial&fpol) -:CompactPolynomial (fpol){} + FCompactPolynomial(const FTensorPolynomial&fpol) + :CompactPolynomial (fpol){} }; diff --git a/mex/sources/korderpert/Dyn_pp/tl/cc/tensor.h b/mex/sources/korderpert/Dyn_pp/tl/cc/tensor.h index 44acdba41..7e8fb0164 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 _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 index; + enum indor{along_row,along_col}; + typedef _index 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 2fd08885f..71d4cc8a4 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);} };