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
time-shift
george 2008-10-20 16:19:58 +00:00
parent 5f5705de80
commit 0337bae230
14 changed files with 1972 additions and 1972 deletions

View File

@ -9,10 +9,10 @@
class UNormalMoments:public TensorContainer<URSingleTensor> {
public:
UNormalMoments(int maxdim,const TwoDMatrix&v);
UNormalMoments(int maxdim,const TwoDMatrix&v);
private:
void generateMoments(int maxdim,const TwoDMatrix&v);
static bool selectEquiv(const Equivalence&e);
void generateMoments(int maxdim,const TwoDMatrix&v);
static bool selectEquiv(const Equivalence&e);
};
/*:2*/
@ -21,7 +21,7 @@ static bool selectEquiv(const Equivalence&e);
class FNormalMoments:public TensorContainer<FRSingleTensor> {
public:
FNormalMoments(const UNormalMoments&moms);
FNormalMoments(const UNormalMoments&moms);
};

View File

@ -12,40 +12,40 @@
class Permutation{
protected:
IntSequence permap;
IntSequence permap;
public:
Permutation(int len)
:permap(len){for(int i= 0;i<len;i++)permap[i]= i;}
Permutation(const Equivalence&e)
:permap(e.getN()){e.trace(permap);}
Permutation(const Equivalence&e,const Permutation&per)
:permap(e.getN()){e.trace(permap,per);}
Permutation(const IntSequence&s)
:permap(s.size()){computeSortingMap(s);};
Permutation(const Permutation&p)
:permap(p.permap){}
Permutation(const Permutation&p1,const Permutation&p2)
:permap(p2.permap){p1.apply(permap);}
Permutation(const Permutation&p,int i)
:permap(p.size(),p.permap,i){}
const Permutation&operator= (const Permutation&p)
{permap= p.permap;return*this;}
bool operator==(const Permutation&p)
{return permap==p.permap;}
int size()const
{return permap.size();}
void print()const
{permap.print();}
void apply(const IntSequence&src,IntSequence&tar)const;
void apply(IntSequence&tar)const;
void inverse();
int tailIdentity()const;
const IntSequence&getMap()const
{return permap;}
IntSequence&getMap()
{return permap;}
Permutation(int len)
:permap(len){for(int i= 0;i<len;i++)permap[i]= i;}
Permutation(const Equivalence&e)
:permap(e.getN()){e.trace(permap);}
Permutation(const Equivalence&e,const Permutation&per)
:permap(e.getN()){e.trace(permap,per);}
Permutation(const IntSequence&s)
:permap(s.size()){computeSortingMap(s);};
Permutation(const Permutation&p)
:permap(p.permap){}
Permutation(const Permutation&p1,const Permutation&p2)
:permap(p2.permap){p1.apply(permap);}
Permutation(const Permutation&p,int i)
:permap(p.size(),p.permap,i){}
const Permutation&operator= (const Permutation&p)
{permap= p.permap;return*this;}
bool operator==(const Permutation&p)
{return permap==p.permap;}
int size()const
{return permap.size();}
void print()const
{permap.print();}
void apply(const IntSequence&src,IntSequence&tar)const;
void apply(IntSequence&tar)const;
void inverse();
int tailIdentity()const;
const IntSequence&getMap()const
{return permap;}
IntSequence&getMap()
{return permap;}
protected:
void computeSortingMap(const IntSequence&s);
void computeSortingMap(const IntSequence&s);
};
@ -54,18 +54,18 @@ void computeSortingMap(const IntSequence&s);
/*3:*/
class PermutationSet{
int order;
int size;
const Permutation**const pers;
int order;
int size;
const Permutation**const pers;
public:
PermutationSet();
PermutationSet(const PermutationSet&ps,int n);
~PermutationSet();
int getNum()const
{return size;}
const Permutation&get(int i)const
{return*(pers[i]);}
vector<const Permutation*> getPreserving(const IntSequence&s)const;
PermutationSet();
PermutationSet(const PermutationSet&ps,int n);
~PermutationSet();
int getNum()const
{return size;}
const Permutation&get(int i)const
{return*(pers[i]);}
vector<const Permutation*> getPreserving(const IntSequence&s)const;
};
@ -74,12 +74,12 @@ vector<const Permutation*> getPreserving(const IntSequence&s)const;
/*4:*/
class PermutationBundle{
vector<PermutationSet*> bundle;
vector<PermutationSet*> bundle;
public:
PermutationBundle(int nmax);
~PermutationBundle();
const PermutationSet&get(int n)const;
void generateUpTo(int nmax);
PermutationBundle(int nmax);
~PermutationBundle();
const PermutationSet&get(int n)const;
void generateUpTo(int nmax);
};
/*:4*/

View File

@ -15,8 +15,8 @@
class SortIntSequence:public IntSequence{
public:
SortIntSequence(const IntSequence&s)
:IntSequence(s){sort();}
SortIntSequence(const IntSequence&s)
:IntSequence(s){sort();}
};
@ -26,31 +26,31 @@ SortIntSequence(const IntSequence&s)
class PerTensorDimens:public TensorDimens{
protected:
Permutation per;
Permutation per;
public:
PerTensorDimens(const Symmetry&s,const IntSequence&nvars,
const Equivalence&e)
:TensorDimens(s,nvars),per(e)
{per.apply(nvmax);}
PerTensorDimens(const TensorDimens&td,const Equivalence&e)
:TensorDimens(td),per(e)
{per.apply(nvmax);}
PerTensorDimens(const TensorDimens&td,const Permutation&p)
:TensorDimens(td),per(p)
{per.apply(nvmax);}
PerTensorDimens(const IntSequence&ss,const IntSequence&coor)
:TensorDimens(ss,SortIntSequence(coor)),per(coor)
{per.apply(nvmax);}
PerTensorDimens(const PerTensorDimens&td)
:TensorDimens(td),per(td.per){}
const PerTensorDimens&operator= (const PerTensorDimens&td)
{TensorDimens::operator= (td);per= td.per;return*this;}
bool operator==(const PerTensorDimens&td)
{return TensorDimens::operator==(td)&&per==td.per;}
int tailIdentity()const
{return per.tailIdentity();}
const Permutation&getPer()const
{return per;}
PerTensorDimens(const Symmetry&s,const IntSequence&nvars,
const Equivalence&e)
:TensorDimens(s,nvars),per(e)
{per.apply(nvmax);}
PerTensorDimens(const TensorDimens&td,const Equivalence&e)
:TensorDimens(td),per(e)
{per.apply(nvmax);}
PerTensorDimens(const TensorDimens&td,const Permutation&p)
:TensorDimens(td),per(p)
{per.apply(nvmax);}
PerTensorDimens(const IntSequence&ss,const IntSequence&coor)
:TensorDimens(ss,SortIntSequence(coor)),per(coor)
{per.apply(nvmax);}
PerTensorDimens(const PerTensorDimens&td)
:TensorDimens(td),per(td.per){}
const PerTensorDimens&operator= (const PerTensorDimens&td)
{TensorDimens::operator= (td);per= td.per;return*this;}
bool operator==(const PerTensorDimens&td)
{return TensorDimens::operator==(td)&&per==td.per;}
int tailIdentity()const
{return per.tailIdentity();}
const Permutation&getPer()const
{return per;}
};
/*:3*/
@ -58,54 +58,54 @@ const Permutation&getPer()const
/*4:*/
class UPSTensor:public UTensor{
const PerTensorDimens tdims;
const PerTensorDimens tdims;
public:
/*5:*/
UPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:UTensor(along_col,PerTensorDimens(td,e).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e)
{kp.mult(a,*this);}
UPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:UTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,kp.getPer()))
{kp.mult(a,*this);}
UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:UTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,p))
{kp.mult(a,*this);}
UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:UTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,Permutation(p,kp.getPer())))
{kp.mult(a,*this);}
/*:5*/
;
UPSTensor(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor,const PerTensorDimens&ptd);
UPSTensor(const UPSTensor&ut)
:UTensor(ut),tdims(ut.tdims){}
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
FTensor&fold()const;
int getOffset(const IntSequence&v)const;
void addTo(FGSTensor&out)const;
void addTo(UGSTensor&out)const;
enum fill_method{first,second};
static fill_method decideFillMethod(const FSSparseTensor&t);
/*5:*/
UPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:UTensor(along_col,PerTensorDimens(td,e).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e)
{kp.mult(a,*this);}
UPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:UTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,kp.getPer()))
{kp.mult(a,*this);}
UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:UTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,p))
{kp.mult(a,*this);}
UPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:UTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,Permutation(e,Permutation(p,kp.getPer())))
{kp.mult(a,*this);}
/*:5*/
;
UPSTensor(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor,const PerTensorDimens&ptd);
UPSTensor(const UPSTensor&ut)
:UTensor(ut),tdims(ut.tdims){}
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
FTensor&fold()const;
int getOffset(const IntSequence&v)const;
void addTo(FGSTensor&out)const;
void addTo(UGSTensor&out)const;
enum fill_method{first,second};
static fill_method decideFillMethod(const FSSparseTensor&t);
private:
int tailIdentitySize()const;
void fillFromSparseOne(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor);
void fillFromSparseTwo(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor);
int tailIdentitySize()const;
void fillFromSparseOne(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor);
void fillFromSparseTwo(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor);
};
/*:4*/
@ -113,30 +113,30 @@ const IntSequence&coor);
/*6:*/
class PerTensorDimens2:public PerTensorDimens{
InducedSymmetries syms;
IntSequence ds;
InducedSymmetries syms;
IntSequence ds;
public:
PerTensorDimens2(const TensorDimens&td,const Equivalence&e,
const Permutation&p)
:PerTensorDimens(td,Permutation(e,p)),
syms(e,p,td.getSym()),
ds(syms.size())
{setDimensionSizes();}
PerTensorDimens2(const TensorDimens&td,const Equivalence&e)
:PerTensorDimens(td,e),
syms(e,td.getSym()),
ds(syms.size())
{setDimensionSizes();}
int numSyms()const
{return(int)syms.size();}
const Symmetry&getSym(int i)const
{return syms[i];}
int calcMaxOffset()const
{return ds.mult();}
int calcOffset(const IntSequence&coor)const;
void print()const;
PerTensorDimens2(const TensorDimens&td,const Equivalence&e,
const Permutation&p)
:PerTensorDimens(td,Permutation(e,p)),
syms(e,p,td.getSym()),
ds(syms.size())
{setDimensionSizes();}
PerTensorDimens2(const TensorDimens&td,const Equivalence&e)
:PerTensorDimens(td,e),
syms(e,td.getSym()),
ds(syms.size())
{setDimensionSizes();}
int numSyms()const
{return(int)syms.size();}
const Symmetry&getSym(int i)const
{return syms[i];}
int calcMaxOffset()const
{return ds.mult();}
int calcOffset(const IntSequence&coor)const;
void print()const;
protected:
void setDimensionSizes();
void setDimensionSizes();
};
/*:6*/
@ -146,46 +146,46 @@ void setDimensionSizes();
template<typename _Ttype> class StackProduct;
class FPSTensor:public FTensor{
const PerTensorDimens2 tdims;
const PerTensorDimens2 tdims;
public:
/*8:*/
FPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:FTensor(along_col,PerTensorDimens(td,e).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e)
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:FTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,kp.getPer())
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:FTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,p)
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:FTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,Permutation(p,kp.getPer()))
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const GSSparseTensor&t,const KronProdAll&kp);
FPSTensor(const FPSTensor&ft)
:FTensor(ft),tdims(ft.tdims){}
/*:8*/
;
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
UTensor&unfold()const;
int getOffset(const IntSequence&v)const;
void addTo(FGSTensor&out)const;
/*8:*/
FPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:FTensor(along_col,PerTensorDimens(td,e).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e)
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:FTensor(along_col,PerTensorDimens(td,Permutation(e,kp.getPer())).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,kp.getPer())
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAll&kp)
:FTensor(along_col,PerTensorDimens(td,Permutation(e,p)).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,p)
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const ConstTwoDMatrix&a,const KronProdAllOptim&kp)
:FTensor(along_col,PerTensorDimens(td,Permutation(e,Permutation(p,kp.getPer()))).getNVX(),
a.nrows(),kp.ncols(),td.dimen()),tdims(td,e,Permutation(p,kp.getPer()))
{kp.mult(a,*this);}
FPSTensor(const TensorDimens&td,const Equivalence&e,const Permutation&p,
const GSSparseTensor&t,const KronProdAll&kp);
FPSTensor(const FPSTensor&ft)
:FTensor(ft),tdims(ft.tdims){}
/*:8*/
;
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
UTensor&unfold()const;
int getOffset(const IntSequence&v)const;
void addTo(FGSTensor&out)const;
};
/*:7*/

View File

@ -16,10 +16,10 @@ using namespace std;
class USubTensor:public URTensor{
public:
USubTensor(const TensorDimens&bdims,const TensorDimens&hdims,
const FGSContainer&cont,const vector<IntSequence> &lst);
void addKronColumn(int i,const vector<const FGSTensor*> &ts,
const IntSequence&pindex);
USubTensor(const TensorDimens&bdims,const TensorDimens&hdims,
const FGSContainer&cont,const vector<IntSequence> &lst);
void addKronColumn(int i,const vector<const FGSTensor*> &ts,
const IntSequence&pindex);
};
/*:2*/

View File

@ -15,20 +15,20 @@
class IrregTensor;
class IrregTensorHeader{
friend class IrregTensor;
int nv;
IntSequence unit_flag;
Vector**const cols;
IntSequence end_seq;
friend class IrregTensor;
int nv;
IntSequence unit_flag;
Vector**const cols;
IntSequence end_seq;
public:
IrregTensorHeader(const StackProduct<FGSTensor> &sp,const IntSequence&c);
~IrregTensorHeader();
int dimen()const
{return unit_flag.size();}
void increment(IntSequence&v)const;
int calcMaxOffset()const;
IrregTensorHeader(const StackProduct<FGSTensor> &sp,const IntSequence&c);
~IrregTensorHeader();
int dimen()const
{return unit_flag.size();}
void increment(IntSequence&v)const;
int calcMaxOffset()const;
private:
IrregTensorHeader(const IrregTensorHeader&);
IrregTensorHeader(const IrregTensorHeader&);
};
@ -37,16 +37,16 @@ IrregTensorHeader(const IrregTensorHeader&);
/*3:*/
class IrregTensor:public Tensor{
const IrregTensorHeader&header;
const IrregTensorHeader&header;
public:
IrregTensor(const IrregTensorHeader&h);
void addTo(FRSingleTensor&out)const;
void increment(IntSequence&v)const
{header.increment(v);}
void decrement(IntSequence&v)const
{TL_RAISE("Not implemented error in IrregTensor::decrement");}
int getOffset(const IntSequence&v)const
{TL_RAISE("Not implemented error in IrregTensor::getOffset");return 0;}
IrregTensor(const IrregTensorHeader&h);
void addTo(FRSingleTensor&out)const;
void increment(IntSequence&v)const
{header.increment(v);}
void decrement(IntSequence&v)const
{TL_RAISE("Not implemented error in IrregTensor::decrement");}
int getOffset(const IntSequence&v)const
{TL_RAISE("Not implemented error in IrregTensor::getOffset");return 0;}
};
/*:3*/

View File

@ -11,30 +11,30 @@
class FRTensor;
class URTensor:public UTensor{
int nv;
int nv;
public:
/*3:*/
URTensor(int c,int nvar,int d)
:UTensor(along_row,IntSequence(d,nvar),
UFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){}
URTensor(const URTensor&ut)
:UTensor(ut),nv(ut.nv){}
URTensor(const FRTensor&ft);
/*:3*/
;
virtual~URTensor(){}
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
FTensor&fold()const;
int getOffset(const IntSequence&v)const;
int nvar()const
{return nv;}
Symmetry getSym()const
{return Symmetry(dimen());}
/*3:*/
URTensor(int c,int nvar,int d)
:UTensor(along_row,IntSequence(d,nvar),
UFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){}
URTensor(const URTensor&ut)
:UTensor(ut),nv(ut.nv){}
URTensor(const FRTensor&ft);
/*:3*/
;
virtual~URTensor(){}
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
FTensor&fold()const;
int getOffset(const IntSequence&v)const;
int nvar()const
{return nv;}
Symmetry getSym()const
{return Symmetry(dimen());}
};
/*:2*/
@ -42,31 +42,31 @@ Symmetry getSym()const
/*4:*/
class FRTensor:public FTensor{
int nv;
int nv;
public:
/*5:*/
FRTensor(int c,int nvar,int d)
:FTensor(along_row,IntSequence(d,nvar),
FFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){}
FRTensor(const FRTensor&ft)
:FTensor(ft),nv(ft.nv){}
FRTensor(const URTensor&ut);
/*:5*/
;
virtual~FRTensor(){}
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
UTensor&unfold()const;
int nvar()const
{return nv;}
int getOffset(const IntSequence&v)const
{return FTensor::getOffset(v,nv);}
Symmetry getSym()const
{return Symmetry(dimen());}
/*5:*/
FRTensor(int c,int nvar,int d)
:FTensor(along_row,IntSequence(d,nvar),
FFSTensor::calcMaxOffset(nvar,d),c,d),nv(nvar){}
FRTensor(const FRTensor&ft)
:FTensor(ft),nv(ft.nv){}
FRTensor(const URTensor&ut);
/*:5*/
;
virtual~FRTensor(){}
void increment(IntSequence&v)const;
void decrement(IntSequence&v)const;
UTensor&unfold()const;
int nvar()const
{return nv;}
int getOffset(const IntSequence&v)const
{return FTensor::getOffset(v,nv);}
Symmetry getSym()const
{return Symmetry(dimen());}
};
/*:4*/
@ -75,14 +75,14 @@ Symmetry getSym()const
class URSingleTensor:public URTensor{
public:
URSingleTensor(int nvar,int d)
:URTensor(1,nvar,d){}
URSingleTensor(const vector<ConstVector> &cols);
URSingleTensor(const ConstVector&v,int d);
URSingleTensor(const URSingleTensor&ut)
:URTensor(ut){}
virtual~URSingleTensor(){}
FTensor&fold()const;
URSingleTensor(int nvar,int d)
:URTensor(1,nvar,d){}
URSingleTensor(const vector<ConstVector> &cols);
URSingleTensor(const ConstVector&v,int d);
URSingleTensor(const URSingleTensor&ut)
:URTensor(ut){}
virtual~URSingleTensor(){}
FTensor&fold()const;
};
/*:6*/
@ -91,12 +91,12 @@ FTensor&fold()const;
class FRSingleTensor:public FRTensor{
public:
FRSingleTensor(int nvar,int d)
:FRTensor(1,nvar,d){}
FRSingleTensor(const URSingleTensor&ut);
FRSingleTensor(const FRSingleTensor&ft)
:FRTensor(ft){}
virtual~FRSingleTensor(){}
FRSingleTensor(int nvar,int d)
:FRTensor(1,nvar,d){}
FRSingleTensor(const URSingleTensor&ut);
FRSingleTensor(const FRSingleTensor&ft)
:FRTensor(ft){}
virtual~FRSingleTensor(){}
};

View File

@ -15,8 +15,8 @@ using namespace std;
/*2:*/
struct ltseq{
bool operator()(const IntSequence&s1,const IntSequence&s2)const
{return s1<s2;}
bool operator()(const IntSequence&s1,const IntSequence&s2)const
{return s1<s2;}
};
/*:2*/
@ -25,46 +25,46 @@ bool operator()(const IntSequence&s1,const IntSequence&s2)const
class SparseTensor{
public:
typedef pair<int,double> Item;
typedef multimap<IntSequence,Item,ltseq> Map;
typedef Map::const_iterator const_iterator;
typedef pair<int,double> Item;
typedef multimap<IntSequence,Item,ltseq> Map;
typedef Map::const_iterator const_iterator;
protected:
typedef Map::iterator iterator;
Map m;
const int dim;
const int nr;
const int nc;
int first_nz_row;
int last_nz_row;
typedef Map::iterator iterator;
Map m;
const int dim;
const int nr;
const int nc;
int first_nz_row;
int last_nz_row;
public:
SparseTensor(int d,int nnr,int nnc)
:dim(d),nr(nnr),nc(nnc),first_nz_row(nr),last_nz_row(-1){}
SparseTensor(const SparseTensor&t)
:m(t.m),dim(t.dim),nr(t.nr),nc(t.nc){}
virtual~SparseTensor(){}
void insert(const IntSequence&s,int r,double c);
const Map&getMap()const
{return m;}
int dimen()const
{return dim;}
int nrows()const
{return nr;}
int ncols()const
{return nc;}
double getFillFactor()const
{return((double)m.size())/(nrows()*ncols());}
double getFoldIndexFillFactor()const;
double getUnfoldIndexFillFactor()const;
int getNumNonZero()const
{return m.size();}
int getFirstNonZeroRow()const
{return first_nz_row;}
int getLastNonZeroRow()const
{return last_nz_row;}
virtual const Symmetry&getSym()const= 0;
void print()const;
bool isFinite()const;
SparseTensor(int d,int nnr,int nnc)
:dim(d),nr(nnr),nc(nnc),first_nz_row(nr),last_nz_row(-1){}
SparseTensor(const SparseTensor&t)
:m(t.m),dim(t.dim),nr(t.nr),nc(t.nc){}
virtual~SparseTensor(){}
void insert(const IntSequence&s,int r,double c);
const Map&getMap()const
{return m;}
int dimen()const
{return dim;}
int nrows()const
{return nr;}
int ncols()const
{return nc;}
double getFillFactor()const
{return((double)m.size())/(nrows()*ncols());}
double getFoldIndexFillFactor()const;
double getUnfoldIndexFillFactor()const;
int getNumNonZero()const
{return m.size();}
int getFirstNonZeroRow()const
{return first_nz_row;}
int getLastNonZeroRow()const
{return last_nz_row;}
virtual const Symmetry&getSym()const= 0;
void print()const;
bool isFinite()const;
}
/*:3*/
@ -73,20 +73,20 @@ bool isFinite()const;
class FSSparseTensor:public SparseTensor{
public:
typedef SparseTensor::const_iterator const_iterator;
typedef SparseTensor::const_iterator const_iterator;
private:
const int nv;
const Symmetry sym;
const int nv;
const Symmetry sym;
public:
FSSparseTensor(int d,int nvar,int r);
FSSparseTensor(const FSSparseTensor&t);
void insert(const IntSequence&s,int r,double c);
void multColumnAndAdd(const Tensor&t,Vector&v)const;
const Symmetry&getSym()const
{return sym;}
int nvar()const
{return nv;}
void print()const;
FSSparseTensor(int d,int nvar,int r);
FSSparseTensor(const FSSparseTensor&t);
void insert(const IntSequence&s,int r,double c);
void multColumnAndAdd(const Tensor&t,Vector&v)const;
const Symmetry&getSym()const
{return sym;}
int nvar()const
{return nv;}
void print()const;
};
@ -96,21 +96,21 @@ void print()const;
class GSSparseTensor:public SparseTensor{
public:
typedef SparseTensor::const_iterator const_iterator;
typedef SparseTensor::const_iterator const_iterator;
private:
const TensorDimens tdims;
const TensorDimens tdims;
public:
GSSparseTensor(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor,const TensorDimens&td);
GSSparseTensor(const GSSparseTensor&t)
:SparseTensor(t),tdims(t.tdims){}
void insert(const IntSequence&s,int r,double c);
const Symmetry&getSym()const
{return tdims.getSym();}
const TensorDimens&getDims()const
{return tdims;}
void print()const;
GSSparseTensor(const FSSparseTensor&t,const IntSequence&ss,
const IntSequence&coor,const TensorDimens&td);
GSSparseTensor(const GSSparseTensor&t)
:SparseTensor(t),tdims(t.tdims){}
void insert(const IntSequence&s,int r,double c);
const Symmetry&getSym()const
{return tdims.getSym();}
const TensorDimens&getDims()const
{return tdims;}
void print()const;
};
/*:5*/

File diff suppressed because it is too large Load Diff

View File

@ -12,464 +12,464 @@
#include <map>
namespace sthread{
using namespace std;
class Empty{};
/*2:*/
template<bool condition,class Then,class Else>
struct IF{
typedef Then RET;
};
template<class Then,class Else>
struct IF<false,Then,Else> {
typedef Else RET;
};
/*:2*/
;
enum{posix,empty};
template<int> class thread_traits;
template<int> class detach_thread;
/*3:*/
template<int thread_impl>
class thread{
typedef thread_traits<thread_impl> _Ttraits;
typedef typename _Ttraits::_Tthread _Tthread;
_Tthread th;
using namespace std;
class Empty{};
/*2:*/
template<bool condition,class Then,class Else>
struct IF{
typedef Then RET;
};
template<class Then,class Else>
struct IF<false,Then,Else> {
typedef Else RET;
};
/*:2*/
;
enum{posix,empty};
template<int> class thread_traits;
template<int> class detach_thread;
/*3:*/
template<int thread_impl>
class thread{
typedef thread_traits<thread_impl> _Ttraits;
typedef typename _Ttraits::_Tthread _Tthread;
_Tthread th;
public:
virtual~thread(){}
_Tthread&getThreadIden()
{return th;}
const _Tthread&getThreadIden()const
{return th;}
virtual void operator()()= 0;
void run()
{_Ttraits::run(this);}
void detach_run()
{_Ttraits::detach_run(this);}
void exit()
{_Ttraits::exit();}
};
/*:3*/
;
/*4:*/
template<int thread_impl>
class thread_group{
typedef thread_traits<thread_impl> _Ttraits;
typedef thread<thread_impl> _Ctype;
list<_Ctype*> tlist;
typedef typename list<_Ctype*> ::iterator iterator;
public:
static int max_parallel_threads;
void insert(_Ctype*c)
{tlist.push_back(c);}
/*5:*/
~thread_group()
{
while(!tlist.empty()){
delete tlist.front();
tlist.pop_front();
}
}
/*:5*/
;
/*7:*/
void run()
{
int rem= tlist.size();
iterator pfirst= tlist.begin();
while(rem> 2*max_parallel_threads){
pfirst= run_portion(pfirst,max_parallel_threads);
rem-= max_parallel_threads;
}
if(rem> max_parallel_threads){
pfirst= run_portion(pfirst,rem/2);
rem-= rem/2;
}
run_portion(pfirst,rem);
}
/*:7*/
;
private:
/*6:*/
iterator run_portion(iterator start,int n)
{
int c= 0;
for(iterator i= start;c<n;++i,c++){
(*i)->run();
}
iterator ret;
c= 0;
for(ret= start;c<n;++ret,c++){
_Ttraits::join(*ret);
}
return ret;
}
/*:6*/
;
};
/*:4*/
;
/*8:*/
template<int thread_impl>
struct thread_traits{
typedef typename IF<thread_impl==posix,pthread_t,Empty> ::RET _Tthread;
typedef thread<thread_impl> _Ctype;
typedef detach_thread<thread_impl> _Dtype;
static void run(_Ctype*c);
static void detach_run(_Dtype*c);
static void exit();
static void join(_Ctype*c);
};
/*:8*/
;
/*9:*/
struct ltmmkey;
typedef pair<const void*,const char*> mmkey;
template<int thread_impl>
struct mutex_traits{
typedef typename IF<thread_impl==posix,pthread_mutex_t,Empty> ::RET _Tmutex;
typedef map<mmkey,pair<_Tmutex,int> ,ltmmkey> mutex_int_map;
static void init(_Tmutex&m);
static void lock(_Tmutex&m);
static void unlock(_Tmutex&m);
};
/*:9*/
;
/*10:*/
struct ltmmkey{
bool operator()(const mmkey&k1,const mmkey&k2)const
{return k1.first<k2.first||
(k1.first==k2.first&&strcmp(k1.second,k2.second)<0);}
};
template<int thread_impl>
class mutex_map
:public mutex_traits<thread_impl> ::mutex_int_map
{
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
typedef mutex_traits<thread_impl> _Mtraits;
typedef pair<_Tmutex,int> mmval;
typedef map<mmkey,mmval,ltmmkey> _Tparent;
typedef typename _Tparent::iterator iterator;
typedef typename _Tparent::value_type _mvtype;
_Tmutex m;
public:
mutex_map()
{_Mtraits::init(m);}
void insert(const void*c,const char*id,const _Tmutex&m)
{_Tparent::insert(_mvtype(mmkey(c,id),mmval(m,0)));}
bool check(const void*c,const char*id)const
{return _Tparent::find(mmkey(c,id))!=_Tparent::end();}
/*11:*/
mmval*get(const void*c,const char*id)
{
iterator it= _Tparent::find(mmkey(c,id));
if(it==_Tparent::end())
return NULL;
return&((*it).second);
}
/*:11*/
;
/*12:*/
void remove(const void*c,const char*id)
{
iterator it= _Tparent::find(mmkey(c,id));
if(it!=_Tparent::end())
erase(it);
}
/*:12*/
;
void lock_map()
{_Mtraits::lock(m);}
void unlock_map()
{_Mtraits::unlock(m);}
};
/*:10*/
;
/*13:*/
template<int thread_impl>
class synchro{
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
typedef mutex_traits<thread_impl> _Mtraits;
public:
typedef mutex_map<thread_impl> mutex_map_t;
private:
const void*caller;
const char*iden;
mutex_map_t&mutmap;
public:
synchro(const void*c,const char*id,mutex_map_t&mmap)
:caller(c),iden(id),mutmap(mmap)
{lock();}
~synchro()
{unlock();}
private:
/*14:*/
void lock(){
mutmap.lock_map();
if(!mutmap.check(caller,iden)){
_Tmutex mut;
_Mtraits::init(mut);
mutmap.insert(caller,iden,mut);
}
mutmap.get(caller,iden)->second++;
mutmap.unlock_map();
_Mtraits::lock(mutmap.get(caller,iden)->first);
}
/*:14*/
;
/*15:*/
void unlock(){
mutmap.lock_map();
if(mutmap.check(caller,iden)){
_Mtraits::unlock(mutmap.get(caller,iden)->first);
mutmap.get(caller,iden)->second--;
if(mutmap.get(caller,iden)->second==0)
mutmap.remove(caller,iden);
}
mutmap.unlock_map();
}
/*:15*/
;
};
/*:13*/
;
/*16:*/
template<int thread_impl>
struct cond_traits{
typedef typename IF<thread_impl==posix,pthread_cond_t,Empty> ::RET _Tcond;
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
static void init(_Tcond&cond);
static void broadcast(_Tcond&cond);
static void wait(_Tcond&cond,_Tmutex&mutex);
static void destroy(_Tcond&cond);
};
/*:16*/
;
/*17:*/
template<int thread_impl>
class condition_counter{
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
typedef typename cond_traits<thread_impl> ::_Tcond _Tcond;
int counter;
_Tmutex mut;
_Tcond cond;
bool changed;
public:
virtual~thread(){}
_Tthread&getThreadIden()
{return th;}
const _Tthread&getThreadIden()const
{return th;}
virtual void operator()()= 0;
void run()
{_Ttraits::run(this);}
void detach_run()
{_Ttraits::detach_run(this);}
void exit()
{_Ttraits::exit();}
};
/*:3*/
;
/*4:*/
template<int thread_impl>
class thread_group{
typedef thread_traits<thread_impl> _Ttraits;
typedef thread<thread_impl> _Ctype;
list<_Ctype*> tlist;
typedef typename list<_Ctype*> ::iterator iterator;
/*18:*/
condition_counter()
:counter(0),changed(true)
{
mutex_traits<thread_impl> ::init(mut);
cond_traits<thread_impl> ::init(cond);
}
/*:18*/
;
/*19:*/
~condition_counter()
{
cond_traits<thread_impl> ::destroy(cond);
}
/*:19*/
;
/*20:*/
void increase()
{
mutex_traits<thread_impl> ::lock(mut);
counter++;
changed= true;
cond_traits<thread_impl> ::broadcast(cond);
mutex_traits<thread_impl> ::unlock(mut);
}
/*:20*/
;
/*21:*/
void decrease()
{
mutex_traits<thread_impl> ::lock(mut);
counter--;
changed= true;
cond_traits<thread_impl> ::broadcast(cond);
mutex_traits<thread_impl> ::unlock(mut);
}
/*:21*/
;
/*22:*/
int waitForChange()
{
mutex_traits<thread_impl> ::lock(mut);
if(!changed){
cond_traits<thread_impl> ::wait(cond,mut);
}
changed= false;
int res= counter;
mutex_traits<thread_impl> ::unlock(mut);
return res;
}
/*:22*/
;
};
/*:17*/
;
/*23:*/
template<int thread_impl>
class detach_thread:public thread<thread_impl> {
public:
static int max_parallel_threads;
void insert(_Ctype*c)
{tlist.push_back(c);}
/*5:*/
~thread_group()
{
while(!tlist.empty()){
delete tlist.front();
tlist.pop_front();
}
}
/*:5*/
;
/*7:*/
void run()
{
int rem= tlist.size();
iterator pfirst= tlist.begin();
while(rem> 2*max_parallel_threads){
pfirst= run_portion(pfirst,max_parallel_threads);
rem-= max_parallel_threads;
}
if(rem> max_parallel_threads){
pfirst= run_portion(pfirst,rem/2);
rem-= rem/2;
}
run_portion(pfirst,rem);
}
/*:7*/
;
private:
/*6:*/
iterator run_portion(iterator start,int n)
{
int c= 0;
for(iterator i= start;c<n;++i,c++){
(*i)->run();
}
iterator ret;
c= 0;
for(ret= start;c<n;++ret,c++){
_Ttraits::join(*ret);
}
return ret;
}
/*:6*/
;
};
/*:4*/
;
/*8:*/
template<int thread_impl>
struct thread_traits{
typedef typename IF<thread_impl==posix,pthread_t,Empty> ::RET _Tthread;
typedef thread<thread_impl> _Ctype;
typedef detach_thread<thread_impl> _Dtype;
static void run(_Ctype*c);
static void detach_run(_Dtype*c);
static void exit();
static void join(_Ctype*c);
};
/*:8*/
;
/*9:*/
struct ltmmkey;
typedef pair<const void*,const char*> mmkey;
template<int thread_impl>
struct mutex_traits{
typedef typename IF<thread_impl==posix,pthread_mutex_t,Empty> ::RET _Tmutex;
typedef map<mmkey,pair<_Tmutex,int> ,ltmmkey> mutex_int_map;
static void init(_Tmutex&m);
static void lock(_Tmutex&m);
static void unlock(_Tmutex&m);
};
/*:9*/
;
/*10:*/
struct ltmmkey{
bool operator()(const mmkey&k1,const mmkey&k2)const
{return k1.first<k2.first||
(k1.first==k2.first&&strcmp(k1.second,k2.second)<0);}
};
template<int thread_impl>
class mutex_map
:public mutex_traits<thread_impl> ::mutex_int_map
{
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
typedef mutex_traits<thread_impl> _Mtraits;
typedef pair<_Tmutex,int> mmval;
typedef map<mmkey,mmval,ltmmkey> _Tparent;
typedef typename _Tparent::iterator iterator;
typedef typename _Tparent::value_type _mvtype;
_Tmutex m;
public:
mutex_map()
{_Mtraits::init(m);}
void insert(const void*c,const char*id,const _Tmutex&m)
{_Tparent::insert(_mvtype(mmkey(c,id),mmval(m,0)));}
bool check(const void*c,const char*id)const
{return _Tparent::find(mmkey(c,id))!=_Tparent::end();}
/*11:*/
mmval*get(const void*c,const char*id)
{
iterator it= _Tparent::find(mmkey(c,id));
if(it==_Tparent::end())
return NULL;
return&((*it).second);
}
/*:11*/
;
/*12:*/
void remove(const void*c,const char*id)
{
iterator it= _Tparent::find(mmkey(c,id));
if(it!=_Tparent::end())
erase(it);
}
/*:12*/
;
void lock_map()
{_Mtraits::lock(m);}
void unlock_map()
{_Mtraits::unlock(m);}
};
/*:10*/
;
/*13:*/
template<int thread_impl>
class synchro{
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
typedef mutex_traits<thread_impl> _Mtraits;
public:
typedef mutex_map<thread_impl> mutex_map_t;
private:
const void*caller;
const char*iden;
mutex_map_t&mutmap;
public:
synchro(const void*c,const char*id,mutex_map_t&mmap)
:caller(c),iden(id),mutmap(mmap)
{lock();}
~synchro()
{unlock();}
private:
/*14:*/
void lock(){
mutmap.lock_map();
if(!mutmap.check(caller,iden)){
_Tmutex mut;
_Mtraits::init(mut);
mutmap.insert(caller,iden,mut);
}
mutmap.get(caller,iden)->second++;
mutmap.unlock_map();
_Mtraits::lock(mutmap.get(caller,iden)->first);
}
/*:14*/
;
/*15:*/
void unlock(){
mutmap.lock_map();
if(mutmap.check(caller,iden)){
_Mtraits::unlock(mutmap.get(caller,iden)->first);
mutmap.get(caller,iden)->second--;
if(mutmap.get(caller,iden)->second==0)
mutmap.remove(caller,iden);
}
mutmap.unlock_map();
}
/*:15*/
;
};
/*:13*/
;
/*16:*/
template<int thread_impl>
struct cond_traits{
typedef typename IF<thread_impl==posix,pthread_cond_t,Empty> ::RET _Tcond;
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
static void init(_Tcond&cond);
static void broadcast(_Tcond&cond);
static void wait(_Tcond&cond,_Tmutex&mutex);
static void destroy(_Tcond&cond);
};
/*:16*/
;
/*17:*/
template<int thread_impl>
class condition_counter{
typedef typename mutex_traits<thread_impl> ::_Tmutex _Tmutex;
typedef typename cond_traits<thread_impl> ::_Tcond _Tcond;
int counter;
_Tmutex mut;
_Tcond cond;
bool changed;
public:
/*18:*/
condition_counter()
:counter(0),changed(true)
{
mutex_traits<thread_impl> ::init(mut);
cond_traits<thread_impl> ::init(cond);
}
/*:18*/
;
/*19:*/
~condition_counter()
{
cond_traits<thread_impl> ::destroy(cond);
}
/*:19*/
;
/*20:*/
void increase()
{
mutex_traits<thread_impl> ::lock(mut);
counter++;
changed= true;
cond_traits<thread_impl> ::broadcast(cond);
mutex_traits<thread_impl> ::unlock(mut);
}
/*:20*/
;
/*21:*/
void decrease()
{
mutex_traits<thread_impl> ::lock(mut);
counter--;
changed= true;
cond_traits<thread_impl> ::broadcast(cond);
mutex_traits<thread_impl> ::unlock(mut);
}
/*:21*/
;
/*22:*/
int waitForChange()
{
mutex_traits<thread_impl> ::lock(mut);
if(!changed){
cond_traits<thread_impl> ::wait(cond,mut);
}
changed= false;
int res= counter;
mutex_traits<thread_impl> ::unlock(mut);
return res;
}
/*:22*/
;
};
/*:17*/
;
/*23:*/
template<int thread_impl>
class detach_thread:public thread<thread_impl> {
public:
condition_counter<thread_impl> *counter;
detach_thread():counter(NULL){}
void installCounter(condition_counter<thread_impl> *c)
{counter= c;}
void run()
{thread_traits<thread_impl> ::detach_run(this);}
};
/*:23*/
;
/*24:*/
template<int thread_impl>
class detach_thread_group{
typedef thread_traits<thread_impl> _Ttraits;
typedef cond_traits<thread_impl> _Ctraits;
typedef detach_thread<thread_impl> _Ctype;
list<_Ctype*> tlist;
typedef typename list<_Ctype*> ::iterator iterator;
condition_counter<thread_impl> counter;
public:
static int max_parallel_threads;
/*25:*/
void insert(_Ctype*c)
{
tlist.push_back(c);
c->installCounter(&counter);
}
/*:25*/
;
/*26:*/
~detach_thread_group()
{
while(!tlist.empty()){
delete tlist.front();
tlist.pop_front();
}
}
/*:26*/
;
/*27:*/
void run()
{
int mpt= max_parallel_threads;
iterator it= tlist.begin();
while(it!=tlist.end()){
if(counter.waitForChange()<mpt){
counter.increase();
(*it)->run();
++it;
}
}
while(counter.waitForChange()> 0){}
}
/*:27*/
;
};
/*:24*/
;
condition_counter<thread_impl> *counter;
detach_thread():counter(NULL){}
void installCounter(condition_counter<thread_impl> *c)
{counter= c;}
void run()
{thread_traits<thread_impl> ::detach_run(this);}
};
/*:23*/
;
/*24:*/
template<int thread_impl>
class detach_thread_group{
typedef thread_traits<thread_impl> _Ttraits;
typedef cond_traits<thread_impl> _Ctraits;
typedef detach_thread<thread_impl> _Ctype;
list<_Ctype*> tlist;
typedef typename list<_Ctype*> ::iterator iterator;
condition_counter<thread_impl> counter;
public:
static int max_parallel_threads;
/*25:*/
void insert(_Ctype*c)
{
tlist.push_back(c);
c->installCounter(&counter);
}
/*:25*/
;
/*26:*/
~detach_thread_group()
{
while(!tlist.empty()){
delete tlist.front();
tlist.pop_front();
}
}
/*:26*/
;
/*27:*/
void run()
{
int mpt= max_parallel_threads;
iterator it= tlist.begin();
while(it!=tlist.end()){
if(counter.waitForChange()<mpt){
counter.increase();
(*it)->run();
++it;
}
}
while(counter.waitForChange()> 0){}
}
/*:27*/
;
};
/*:24*/
;
#ifdef POSIX_THREADS
/*28:*/
typedef detach_thread<posix> PosixThread;
typedef detach_thread_group<posix> PosixThreadGroup;
typedef synchro<posix> posix_synchro;
class PosixSynchro:public posix_synchro{
public:
PosixSynchro(const void*c,const char*id);
};
/*28:*/
typedef detach_thread<posix> PosixThread;
typedef detach_thread_group<posix> PosixThreadGroup;
typedef synchro<posix> posix_synchro;
class PosixSynchro:public posix_synchro{
public:
PosixSynchro(const void*c,const char*id);
};
#define THREAD sthread::PosixThread
#define THREAD_GROUP sthread::PosixThreadGroup
#define SYNCHRO sthread::PosixSynchro
/*:28*/
;
/*:28*/
;
#else
/*29:*/
typedef thread<empty> NoThread;
typedef thread_group<empty> NoThreadGroup;
typedef synchro<empty> no_synchro;
class NoSynchro{
public:
NoSynchro(const void*c,const char*id){}
~NoSynchro(){}
};
/*29:*/
typedef thread<empty> NoThread;
typedef thread_group<empty> NoThreadGroup;
typedef synchro<empty> no_synchro;
class NoSynchro{
public:
NoSynchro(const void*c,const char*id){}
~NoSynchro(){}
};
#define THREAD sthread::NoThread
#define THREAD_GROUP sthread::NoThreadGroup
#define SYNCHRO sthread::NoSynchro
/*:29*/
;
/*:29*/
;
#endif
};

View File

@ -13,45 +13,45 @@
class Symmetry:public IntSequence{
public:
/*3:*/
Symmetry(int len,const char*dummy)
:IntSequence(len,0){}
Symmetry(int i1)
:IntSequence(1,i1){}
Symmetry(int i1,int i2)
:IntSequence(2){operator[](0)= i1;operator[](1)= i2;}
Symmetry(int i1,int i2,int i3)
:IntSequence(3)
{
operator[](0)= i1;
operator[](1)= i2;
operator[](2)= i3;
}
Symmetry(int i1,int i2,int i3,int i4)
:IntSequence(4)
{
operator[](0)= i1;
operator[](1)= i2;
operator[](2)= i3;
operator[](3)= i4;
}
Symmetry(const Symmetry&s)
:IntSequence(s){}
Symmetry(const Symmetry&s,const OrdSequence&cl)
:IntSequence(s,cl.getData()){}
Symmetry(Symmetry&s,int len)
:IntSequence(s,s.size()-len,s.size()){}
Symmetry(const IntSequence&s);
/*:3*/
;
int num()const
{return size();}
int dimen()const
{return sum();}
int findClass(int i)const;
bool isFull()const;
/*3:*/
Symmetry(int len,const char*dummy)
:IntSequence(len,0){}
Symmetry(int i1)
:IntSequence(1,i1){}
Symmetry(int i1,int i2)
:IntSequence(2){operator[](0)= i1;operator[](1)= i2;}
Symmetry(int i1,int i2,int i3)
:IntSequence(3)
{
operator[](0)= i1;
operator[](1)= i2;
operator[](2)= i3;
}
Symmetry(int i1,int i2,int i3,int i4)
:IntSequence(4)
{
operator[](0)= i1;
operator[](1)= i2;
operator[](2)= i3;
operator[](3)= i4;
}
Symmetry(const Symmetry&s)
:IntSequence(s){}
Symmetry(const Symmetry&s,const OrdSequence&cl)
:IntSequence(s,cl.getData()){}
Symmetry(Symmetry&s,int len)
:IntSequence(s,s.size()-len,s.size()){}
Symmetry(const IntSequence&s);
/*:3*/
;
int num()const
{return size();}
int dimen()const
{return sum();}
int findClass(int i)const;
bool isFull()const;
};
/*:2*/
@ -59,21 +59,21 @@ bool isFull()const;
/*4:*/
class SymmetrySet{
Symmetry run;
int dim;
Symmetry run;
int dim;
public:
SymmetrySet(int d,int length)
:run(length,""),dim(d){}
SymmetrySet(SymmetrySet&s,int d)
:run(s.run,s.size()-1),dim(d){}
int dimen()const
{return dim;}
const Symmetry&sym()const
{return run;}
Symmetry&sym()
{return run;}
int size()const
{return run.size();}
SymmetrySet(int d,int length)
:run(length,""),dim(d){}
SymmetrySet(SymmetrySet&s,int d)
:run(s.run,s.size()-1),dim(d){}
int dimen()const
{return dim;}
const Symmetry&sym()const
{return run;}
Symmetry&sym()
{return run;}
int size()const
{return run.size();}
};
/*:4*/
@ -81,18 +81,18 @@ int size()const
/*5:*/
class symiterator{
SymmetrySet&s;
symiterator*subit;
SymmetrySet*subs;
bool end_flag;
SymmetrySet&s;
symiterator*subit;
SymmetrySet*subs;
bool end_flag;
public:
symiterator(SymmetrySet&ss);
~symiterator();
symiterator&operator++();
bool isEnd()const
{return end_flag;}
const Symmetry&operator*()const
{return s.sym();}
symiterator(SymmetrySet&ss);
~symiterator();
symiterator&operator++();
bool isEnd()const
{return end_flag;}
const Symmetry&operator*()const
{return s.sym();}
};
@ -102,9 +102,9 @@ const Symmetry&operator*()const
class InducedSymmetries:public vector<Symmetry> {
public:
InducedSymmetries(const Equivalence&e,const Symmetry&s);
InducedSymmetries(const Equivalence&e,const Permutation&p,const Symmetry&s);
void print()const;
InducedSymmetries(const Equivalence&e,const Symmetry&s);
InducedSymmetries(const Equivalence&e,const Permutation&p,const Symmetry&s);
void print()const;
};

View File

@ -17,8 +17,8 @@
/*2:*/
struct ltsym{
bool operator()(const Symmetry&s1,const Symmetry&s2)const
{return s1<s2;}
bool operator()(const Symmetry&s1,const Symmetry&s2)const
{return s1<s2;}
};
/*:2*/
@ -27,224 +27,224 @@ bool operator()(const Symmetry&s1,const Symmetry&s2)const
template<class _Ttype> class TensorContainer{
protected:
typedef const _Ttype*_const_ptr;
typedef _Ttype*_ptr;
typedef map<Symmetry,_ptr,ltsym> _Map;
typedef typename _Map::value_type _mvtype;
typedef const _Ttype*_const_ptr;
typedef _Ttype*_ptr;
typedef map<Symmetry,_ptr,ltsym> _Map;
typedef typename _Map::value_type _mvtype;
public:
typedef typename _Map::iterator iterator;
typedef typename _Map::const_iterator const_iterator;
typedef typename _Map::iterator iterator;
typedef typename _Map::const_iterator const_iterator;
private:
int n;
_Map m;
int n;
_Map m;
protected:
const EquivalenceBundle&ebundle;
const EquivalenceBundle&ebundle;
public:
TensorContainer(int nn)
:n(nn),ebundle(*(tls.ebundle)){}
/*5:*/
TensorContainer(const TensorContainer<_Ttype> &c)
:n(c.n),m(),ebundle(c.ebundle)
{
for(const_iterator it= c.m.begin();it!=c.m.end();++it){
_Ttype*ten= new _Ttype(*((*it).second));
insert(ten);
}
}
/*:5*/
;
/*6:*/
TensorContainer(int first_row,int num,TensorContainer<_Ttype> &c)
:n(c.n),ebundle(*(tls.ebundle))
{
for(iterator it= c.m.begin();it!=c.m.end();++it){
_Ttype*t= new _Ttype(first_row,num,*((*it).second));
insert(t);
}
}
/*:6*/
;
/*7:*/
_const_ptr get(const Symmetry&s)const
{
TL_RAISE_IF(s.num()!=num(),
"Incompatible symmetry lookup in TensorContainer::get");
const_iterator it= m.find(s);
if(it==m.end()){
TL_RAISE("Symmetry not found in TensorContainer::get");
return NULL;
}else{
return(*it).second;
}
}
_ptr get(const Symmetry&s)
{
TL_RAISE_IF(s.num()!=num(),
"Incompatible symmetry lookup in TensorContainer::get");
iterator it= m.find(s);
if(it==m.end()){
TL_RAISE("Symmetry not found in TensorContainer::get");
return NULL;
}else{
return(*it).second;
}
}
/*:7*/
;
/*8:*/
bool check(const Symmetry&s)const
{
TL_RAISE_IF(s.num()!=num(),
"Incompatible symmetry lookup in TensorContainer::check");
const_iterator it= m.find(s);
return it!=m.end();
}
/*:8*/
;
/*9:*/
void insert(_ptr t)
{
TL_RAISE_IF(t->getSym().num()!=num(),
"Incompatible symmetry insertion in TensorContainer::insert");
TL_RAISE_IF(check(t->getSym()),
"Tensor already in container in TensorContainer::insert");
m.insert(_mvtype(t->getSym(),t));
if(!t->isFinite()){
throw TLException(__FILE__,__LINE__,"NaN or Inf asserted in TensorContainer::insert");
}
}
/*:9*/
;
/*10:*/
void remove(const Symmetry&s)
{
iterator it= m.find(s);
if(it!=m.end()){
_ptr t= (*it).second;
m.erase(it);
delete t;
}
}
/*:10*/
;
/*11:*/
void clear()
{
while(!m.empty()){
delete(*(m.begin())).second;
m.erase(m.begin());
}
}
/*:11*/
;
/*15:*/
vector<_const_ptr>
fetchTensors(const Symmetry&rsym,const Equivalence&e)const
{
vector<_const_ptr> res(e.numClasses());
int i= 0;
for(Equivalence::const_seqit it= e.begin();
it!=e.end();++it,i++){
Symmetry s(rsym,*it);
res[i]= get(s);
}
return res;
}
/*:15*/
;
/*12:*/
int getMaxDim()const
{
int res= -1;
for(const_iterator run= m.begin();run!=m.end();++run){
int dim= (*run).first.dimen();
if(dim> res)
res= dim;
}
return res;
}
/*:12*/
;
/*13:*/
void print()const
{
printf("Tensor container: nvars=%d, tensors=%d\n",n,m.size());
for(const_iterator it= m.begin();it!=m.end();++it){
printf("Symmetry: ");
(*it).first.print();
((*it).second)->print();
}
}
/*:13*/
;
/*14:*/
void writeMat4(FILE*fd,const char*prefix)const
{
for(const_iterator it= begin();it!=end();++it){
char lname[100];
sprintf(lname,"%s_g",prefix);
const Symmetry&sym= (*it).first;
for(int i= 0;i<sym.num();i++){
char tmp[10];
sprintf(tmp,"_%d",sym[i]);
strcat(lname,tmp);
}
ConstTwoDMatrix m(*((*it).second));
m.writeMat4(fd,lname);
}
}
/*:14*/
;
virtual~TensorContainer()
{clear();}
/*4:*/
int num()const
{return n;}
const EquivalenceBundle&getEqBundle()const
{return ebundle;}
const_iterator begin()const
{return m.begin();}
const_iterator end()const
{return m.end();}
iterator begin()
{return m.begin();}
iterator end()
{return m.end();}
/*:4*/
;
TensorContainer(int nn)
:n(nn),ebundle(*(tls.ebundle)){}
/*5:*/
TensorContainer(const TensorContainer<_Ttype> &c)
:n(c.n),m(),ebundle(c.ebundle)
{
for(const_iterator it= c.m.begin();it!=c.m.end();++it){
_Ttype*ten= new _Ttype(*((*it).second));
insert(ten);
}
}
/*:5*/
;
/*6:*/
TensorContainer(int first_row,int num,TensorContainer<_Ttype> &c)
:n(c.n),ebundle(*(tls.ebundle))
{
for(iterator it= c.m.begin();it!=c.m.end();++it){
_Ttype*t= new _Ttype(first_row,num,*((*it).second));
insert(t);
}
}
/*:6*/
;
/*7:*/
_const_ptr get(const Symmetry&s)const
{
TL_RAISE_IF(s.num()!=num(),
"Incompatible symmetry lookup in TensorContainer::get");
const_iterator it= m.find(s);
if(it==m.end()){
TL_RAISE("Symmetry not found in TensorContainer::get");
return NULL;
}else{
return(*it).second;
}
}
_ptr get(const Symmetry&s)
{
TL_RAISE_IF(s.num()!=num(),
"Incompatible symmetry lookup in TensorContainer::get");
iterator it= m.find(s);
if(it==m.end()){
TL_RAISE("Symmetry not found in TensorContainer::get");
return NULL;
}else{
return(*it).second;
}
}
/*:7*/
;
/*8:*/
bool check(const Symmetry&s)const
{
TL_RAISE_IF(s.num()!=num(),
"Incompatible symmetry lookup in TensorContainer::check");
const_iterator it= m.find(s);
return it!=m.end();
}
/*:8*/
;
/*9:*/
void insert(_ptr t)
{
TL_RAISE_IF(t->getSym().num()!=num(),
"Incompatible symmetry insertion in TensorContainer::insert");
TL_RAISE_IF(check(t->getSym()),
"Tensor already in container in TensorContainer::insert");
m.insert(_mvtype(t->getSym(),t));
if(!t->isFinite()){
throw TLException(__FILE__,__LINE__,"NaN or Inf asserted in TensorContainer::insert");
}
}
/*:9*/
;
/*10:*/
void remove(const Symmetry&s)
{
iterator it= m.find(s);
if(it!=m.end()){
_ptr t= (*it).second;
m.erase(it);
delete t;
}
}
/*:10*/
;
/*11:*/
void clear()
{
while(!m.empty()){
delete(*(m.begin())).second;
m.erase(m.begin());
}
}
/*:11*/
;
/*15:*/
vector<_const_ptr>
fetchTensors(const Symmetry&rsym,const Equivalence&e)const
{
vector<_const_ptr> res(e.numClasses());
int i= 0;
for(Equivalence::const_seqit it= e.begin();
it!=e.end();++it,i++){
Symmetry s(rsym,*it);
res[i]= get(s);
}
return res;
}
/*:15*/
;
/*12:*/
int getMaxDim()const
{
int res= -1;
for(const_iterator run= m.begin();run!=m.end();++run){
int dim= (*run).first.dimen();
if(dim> res)
res= dim;
}
return res;
}
/*:12*/
;
/*13:*/
void print()const
{
printf("Tensor container: nvars=%d, tensors=%d\n",n,m.size());
for(const_iterator it= m.begin();it!=m.end();++it){
printf("Symmetry: ");
(*it).first.print();
((*it).second)->print();
}
}
/*:13*/
;
/*14:*/
void writeMat4(FILE*fd,const char*prefix)const
{
for(const_iterator it= begin();it!=end();++it){
char lname[100];
sprintf(lname,"%s_g",prefix);
const Symmetry&sym= (*it).first;
for(int i= 0;i<sym.num();i++){
char tmp[10];
sprintf(tmp,"_%d",sym[i]);
strcat(lname,tmp);
}
ConstTwoDMatrix m(*((*it).second));
m.writeMat4(fd,lname);
}
}
/*:14*/
;
virtual~TensorContainer()
{clear();}
/*4:*/
int num()const
{return n;}
const EquivalenceBundle&getEqBundle()const
{return ebundle;}
const_iterator begin()const
{return m.begin();}
const_iterator end()const
{return m.end();}
iterator begin()
{return m.begin();}
iterator end()
{return m.end();}
/*:4*/
;
};
/*:3*/
@ -254,12 +254,12 @@ iterator end()
class FGSContainer;
class UGSContainer:public TensorContainer<UGSTensor> {
public:
UGSContainer(int nn)
:TensorContainer<UGSTensor> (nn){}
UGSContainer(const UGSContainer&uc)
:TensorContainer<UGSTensor> (uc){}
UGSContainer(const FGSContainer&c);
void multAndAdd(const UGSTensor&t,UGSTensor&out)const;
UGSContainer(int nn)
:TensorContainer<UGSTensor> (nn){}
UGSContainer(const UGSContainer&uc)
:TensorContainer<UGSTensor> (uc){}
UGSContainer(const FGSContainer&c);
void multAndAdd(const UGSTensor&t,UGSTensor&out)const;
};
@ -268,20 +268,20 @@ void multAndAdd(const UGSTensor&t,UGSTensor&out)const;
/*17:*/
class FGSContainer:public TensorContainer<FGSTensor> {
static const int num_one_time;
static const int num_one_time;
public:
FGSContainer(int nn)
:TensorContainer<FGSTensor> (nn){}
FGSContainer(const FGSContainer&fc)
:TensorContainer<FGSTensor> (fc){}
FGSContainer(const UGSContainer&c);
void multAndAdd(const FGSTensor&t,FGSTensor&out)const;
void multAndAdd(const UGSTensor&t,FGSTensor&out)const;
FGSContainer(int nn)
:TensorContainer<FGSTensor> (nn){}
FGSContainer(const FGSContainer&fc)
:TensorContainer<FGSTensor> (fc){}
FGSContainer(const UGSContainer&c);
void multAndAdd(const FGSTensor&t,FGSTensor&out)const;
void multAndAdd(const UGSTensor&t,FGSTensor&out)const;
private:
static Tensor::index
getIndices(int num,vector<IntSequence> &out,
const Tensor::index&start,
const Tensor::index&end);
static Tensor::index
getIndices(int num,vector<IntSequence> &out,
const Tensor::index&start,
const Tensor::index&end);
};

View File

@ -8,16 +8,16 @@
/*2:*/
class PowerProvider{
Vector origv;
URSingleTensor*ut;
FRSingleTensor*ft;
int nv;
Vector origv;
URSingleTensor*ut;
FRSingleTensor*ft;
int nv;
public:
PowerProvider(const ConstVector&v)
:origv(v),ut(NULL),ft(NULL),nv(v.length()){}
~PowerProvider();
const URSingleTensor&getNext(const URSingleTensor*dummy);
const FRSingleTensor&getNext(const FRSingleTensor*dummy);
PowerProvider(const ConstVector&v)
:origv(v),ut(NULL),ft(NULL),nv(v.length()){}
~PowerProvider();
const URSingleTensor&getNext(const URSingleTensor*dummy);
const FRSingleTensor&getNext(const FRSingleTensor*dummy);
};
/*:2*/
@ -26,226 +26,226 @@ const FRSingleTensor&getNext(const FRSingleTensor*dummy);
template<class _Ttype,class _TGStype,class _Stype>
class TensorPolynomial:public TensorContainer<_Ttype> {
int nr;
int nv;
int maxdim;
typedef TensorContainer<_Ttype> _Tparent;
typedef typename _Tparent::_ptr _ptr;
int nr;
int nv;
int maxdim;
typedef TensorContainer<_Ttype> _Tparent;
typedef typename _Tparent::_ptr _ptr;
public:
TensorPolynomial(int rows,int vars)
:TensorContainer<_Ttype> (1),
nr(rows),nv(vars),maxdim(0){}
TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,int k)
:TensorContainer<_Ttype> (tp),
nr(tp.nr),nv(tp.nv),maxdim(0){derivative(k);}
TensorPolynomial(int first_row,int num,TensorPolynomial<_Ttype,_TGStype,_Stype> &tp)
:TensorContainer<_Ttype> (first_row,num,tp),
nr(num),nv(tp.nv),maxdim(tp.maxdim){}
/*4:*/
TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,const Vector&xval)
:TensorContainer<_Ttype> (1),
nr(tp.nrows()),nv(tp.nvars()-xval.length()),maxdim(0)
{
TL_RAISE_IF(nvars()<0,
"Length of xval too big in TensorPolynomial contract constructor");
IntSequence ss(2);ss[0]= xval.length();ss[1]= nvars();
IntSequence pp(2);pp[0]= 0;pp[1]= 1;
/*5:*/
PowerProvider pwp(xval);
for(int i= 1;i<=tp.maxdim;i++){
const _Stype&xpow= pwp.getNext((const _Stype*)NULL);
for(int j= 0;j<=tp.maxdim-i;j++){
if(tp.check(Symmetry(i+j))){
/*7:*/
_Ttype*ten;
if(_Tparent::check(Symmetry(j))){
ten= _Tparent::get(Symmetry(j));
}else{
ten= new _Ttype(nrows(),nvars(),j);
ten->zeros();
insert(ten);
}
/*:7*/
;
Symmetry sym(i,j);
IntSequence coor(sym,pp);
_TGStype slice(*(tp.get(Symmetry(i+j))),ss,coor,TensorDimens(sym,ss));
slice.mult(Tensor::noverk(i+j,j));
_TGStype tmp(*ten);
slice.contractAndAdd(0,tmp,xpow);
}
}
}
/*:5*/
;
/*6:*/
for(int j= 0;j<=tp.maxdim;j++){
if(tp.check(Symmetry(j))){
/*7:*/
_Ttype*ten;
if(_Tparent::check(Symmetry(j))){
ten= _Tparent::get(Symmetry(j));
}else{
ten= new _Ttype(nrows(),nvars(),j);
ten->zeros();
insert(ten);
}
/*:7*/
;
Symmetry sym(0,j);
IntSequence coor(sym,pp);
_TGStype slice(*(tp.get(Symmetry(j))),ss,coor,TensorDimens(sym,ss));
ten->add(1.0,slice);
}
}
/*:6*/
;
}
/*:4*/
;
TensorPolynomial(const TensorPolynomial&tp)
:TensorContainer<_Ttype> (tp),nr(tp.nr),nv(tp.nv),maxdim(tp.maxdim){}
int nrows()const
{return nr;}
int nvars()const
{return nv;}
/*8:*/
void evalTrad(Vector&out,const ConstVector&v)const
{
if(_Tparent::check(Symmetry(0)))
out= _Tparent::get(Symmetry(0))->getData();
else
out.zeros();
PowerProvider pp(v);
for(int d= 1;d<=maxdim;d++){
const _Stype&p= pp.getNext((const _Stype*)NULL);
Symmetry cs(d);
if(_Tparent::check(cs)){
const _Ttype*t= _Tparent::get(cs);
t->multaVec(out,p.getData());
}
}
}
/*:8*/
;
/*9:*/
void evalHorner(Vector&out,const ConstVector&v)const
{
if(_Tparent::check(Symmetry(0)))
out= _Tparent::get(Symmetry(0))->getData();
else
out.zeros();
if(maxdim==0)
return;
_Ttype*last;
if(maxdim==1)
last= new _Ttype(*(_Tparent::get(Symmetry(1))));
else
last= new _Ttype(*(_Tparent::get(Symmetry(maxdim))),v);
for(int d= maxdim-1;d>=1;d--){
Symmetry cs(d);
if(_Tparent::check(cs)){
const _Ttype*nt= _Tparent::get(cs);
last->add(1.0,ConstTwoDMatrix(*nt));
}
if(d> 1){
_Ttype*new_last= new _Ttype(*last,v);
delete last;
last= new_last;
}
}
last->multaVec(out,v);
delete last;
}
/*:9*/
;
/*10:*/
void insert(_ptr t)
{
TL_RAISE_IF(t->nrows()!=nr,
"Wrong number of rows in TensorPolynomial::insert");
TL_RAISE_IF(t->nvar()!=nv,
"Wrong number of variables in TensorPolynomial::insert");
TensorContainer<_Ttype> ::insert(t);
if(maxdim<t->dimen())
maxdim= t->dimen();
}
/*:10*/
;
/*11:*/
void derivative(int k)
{
for(int d= 1;d<=maxdim;d++){
if(_Tparent::check(Symmetry(d))){
_Ttype*ten= _Tparent::get(Symmetry(d));
ten->mult((double)max((d-k),0));
}
}
}
/*:11*/
;
/*12:*/
_Ttype*evalPartially(int s,const ConstVector&v)
{
TL_RAISE_IF(v.length()!=nvars(),
"Wrong length of vector for TensorPolynomial::evalPartially");
_Ttype*res= new _Ttype(nrows(),nvars(),s);
res->zeros();
int sfac= 1;
for(int i= 1;i<=s;i++)
sfac*= i;
if(_Tparent::check(Symmetry(s)))
res->add(1.0/sfac,*(_Tparent::get(Symmetry(s))));
int dfac= sfac*(s+1);
for(int d= s+1;d<=maxdim;d++,dfac*= d){
if(_Tparent::check(Symmetry(d))){
const _Ttype&ltmp= *(_Tparent::get(Symmetry(d)));
_Ttype*last= new _Ttype(ltmp);
for(int j= 0;j<d-s;j++){
_Ttype*newlast= new _Ttype(*last,v);
delete last;
last= newlast;
}
res->add(1.0/dfac,*last);
delete last;
}
}
return res;
}
/*:12*/
;
TensorPolynomial(int rows,int vars)
:TensorContainer<_Ttype> (1),
nr(rows),nv(vars),maxdim(0){}
TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,int k)
:TensorContainer<_Ttype> (tp),
nr(tp.nr),nv(tp.nv),maxdim(0){derivative(k);}
TensorPolynomial(int first_row,int num,TensorPolynomial<_Ttype,_TGStype,_Stype> &tp)
:TensorContainer<_Ttype> (first_row,num,tp),
nr(num),nv(tp.nv),maxdim(tp.maxdim){}
/*4:*/
TensorPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &tp,const Vector&xval)
:TensorContainer<_Ttype> (1),
nr(tp.nrows()),nv(tp.nvars()-xval.length()),maxdim(0)
{
TL_RAISE_IF(nvars()<0,
"Length of xval too big in TensorPolynomial contract constructor");
IntSequence ss(2);ss[0]= xval.length();ss[1]= nvars();
IntSequence pp(2);pp[0]= 0;pp[1]= 1;
/*5:*/
PowerProvider pwp(xval);
for(int i= 1;i<=tp.maxdim;i++){
const _Stype&xpow= pwp.getNext((const _Stype*)NULL);
for(int j= 0;j<=tp.maxdim-i;j++){
if(tp.check(Symmetry(i+j))){
/*7:*/
_Ttype*ten;
if(_Tparent::check(Symmetry(j))){
ten= _Tparent::get(Symmetry(j));
}else{
ten= new _Ttype(nrows(),nvars(),j);
ten->zeros();
insert(ten);
}
/*:7*/
;
Symmetry sym(i,j);
IntSequence coor(sym,pp);
_TGStype slice(*(tp.get(Symmetry(i+j))),ss,coor,TensorDimens(sym,ss));
slice.mult(Tensor::noverk(i+j,j));
_TGStype tmp(*ten);
slice.contractAndAdd(0,tmp,xpow);
}
}
}
/*:5*/
;
/*6:*/
for(int j= 0;j<=tp.maxdim;j++){
if(tp.check(Symmetry(j))){
/*7:*/
_Ttype*ten;
if(_Tparent::check(Symmetry(j))){
ten= _Tparent::get(Symmetry(j));
}else{
ten= new _Ttype(nrows(),nvars(),j);
ten->zeros();
insert(ten);
}
/*:7*/
;
Symmetry sym(0,j);
IntSequence coor(sym,pp);
_TGStype slice(*(tp.get(Symmetry(j))),ss,coor,TensorDimens(sym,ss));
ten->add(1.0,slice);
}
}
/*:6*/
;
}
/*:4*/
;
TensorPolynomial(const TensorPolynomial&tp)
:TensorContainer<_Ttype> (tp),nr(tp.nr),nv(tp.nv),maxdim(tp.maxdim){}
int nrows()const
{return nr;}
int nvars()const
{return nv;}
/*8:*/
void evalTrad(Vector&out,const ConstVector&v)const
{
if(_Tparent::check(Symmetry(0)))
out= _Tparent::get(Symmetry(0))->getData();
else
out.zeros();
PowerProvider pp(v);
for(int d= 1;d<=maxdim;d++){
const _Stype&p= pp.getNext((const _Stype*)NULL);
Symmetry cs(d);
if(_Tparent::check(cs)){
const _Ttype*t= _Tparent::get(cs);
t->multaVec(out,p.getData());
}
}
}
/*:8*/
;
/*9:*/
void evalHorner(Vector&out,const ConstVector&v)const
{
if(_Tparent::check(Symmetry(0)))
out= _Tparent::get(Symmetry(0))->getData();
else
out.zeros();
if(maxdim==0)
return;
_Ttype*last;
if(maxdim==1)
last= new _Ttype(*(_Tparent::get(Symmetry(1))));
else
last= new _Ttype(*(_Tparent::get(Symmetry(maxdim))),v);
for(int d= maxdim-1;d>=1;d--){
Symmetry cs(d);
if(_Tparent::check(cs)){
const _Ttype*nt= _Tparent::get(cs);
last->add(1.0,ConstTwoDMatrix(*nt));
}
if(d> 1){
_Ttype*new_last= new _Ttype(*last,v);
delete last;
last= new_last;
}
}
last->multaVec(out,v);
delete last;
}
/*:9*/
;
/*10:*/
void insert(_ptr t)
{
TL_RAISE_IF(t->nrows()!=nr,
"Wrong number of rows in TensorPolynomial::insert");
TL_RAISE_IF(t->nvar()!=nv,
"Wrong number of variables in TensorPolynomial::insert");
TensorContainer<_Ttype> ::insert(t);
if(maxdim<t->dimen())
maxdim= t->dimen();
}
/*:10*/
;
/*11:*/
void derivative(int k)
{
for(int d= 1;d<=maxdim;d++){
if(_Tparent::check(Symmetry(d))){
_Ttype*ten= _Tparent::get(Symmetry(d));
ten->mult((double)max((d-k),0));
}
}
}
/*:11*/
;
/*12:*/
_Ttype*evalPartially(int s,const ConstVector&v)
{
TL_RAISE_IF(v.length()!=nvars(),
"Wrong length of vector for TensorPolynomial::evalPartially");
_Ttype*res= new _Ttype(nrows(),nvars(),s);
res->zeros();
int sfac= 1;
for(int i= 1;i<=s;i++)
sfac*= i;
if(_Tparent::check(Symmetry(s)))
res->add(1.0/sfac,*(_Tparent::get(Symmetry(s))));
int dfac= sfac*(s+1);
for(int d= s+1;d<=maxdim;d++,dfac*= d){
if(_Tparent::check(Symmetry(d))){
const _Ttype&ltmp= *(_Tparent::get(Symmetry(d)));
_Ttype*last= new _Ttype(ltmp);
for(int j= 0;j<d-s;j++){
_Ttype*newlast= new _Ttype(*last,v);
delete last;
last= newlast;
}
res->add(1.0/dfac,*last);
delete last;
}
}
return res;
}
/*:12*/
;
};
@ -256,15 +256,15 @@ return res;
class FTensorPolynomial;
class UTensorPolynomial:public TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> {
public:
UTensorPolynomial(int rows,int vars)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (rows,vars){}
UTensorPolynomial(const UTensorPolynomial&up,int k)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (up,k){}
UTensorPolynomial(const FTensorPolynomial&fp);
UTensorPolynomial(const UTensorPolynomial&tp,const Vector&xval)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (tp,xval){}
UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (first_row,num,tp){}
UTensorPolynomial(int rows,int vars)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (rows,vars){}
UTensorPolynomial(const UTensorPolynomial&up,int k)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (up,k){}
UTensorPolynomial(const FTensorPolynomial&fp);
UTensorPolynomial(const UTensorPolynomial&tp,const Vector&xval)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (tp,xval){}
UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp)
:TensorPolynomial<UFSTensor,UGSTensor,URSingleTensor> (first_row,num,tp){}
};
/*:13*/
@ -273,15 +273,15 @@ UTensorPolynomial(int first_row,int num,UTensorPolynomial&tp)
class FTensorPolynomial:public TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> {
public:
FTensorPolynomial(int rows,int vars)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (rows,vars){}
FTensorPolynomial(const FTensorPolynomial&fp,int k)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fp,k){}
FTensorPolynomial(const UTensorPolynomial&up);
FTensorPolynomial(const FTensorPolynomial&tp,const Vector&xval)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (tp,xval){}
FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (first_row,num,tp){}
FTensorPolynomial(int rows,int vars)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (rows,vars){}
FTensorPolynomial(const FTensorPolynomial&fp,int k)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fp,k){}
FTensorPolynomial(const UTensorPolynomial&up);
FTensorPolynomial(const FTensorPolynomial&tp,const Vector&xval)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (tp,xval){}
FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp)
:TensorPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (first_row,num,tp){}
};
/*:14*/
@ -291,61 +291,61 @@ FTensorPolynomial(int first_row,int num,FTensorPolynomial&tp)
template<class _Ttype,class _TGStype,class _Stype>
class CompactPolynomial:public _Ttype{
public:
/*16:*/
CompactPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &pol)
:_Ttype(pol.nrows(),pol.nvars()+1,pol.getMaxDim())
{
_Ttype::zeros();
IntSequence dumnvs(2);
dumnvs[0]= 1;
dumnvs[1]= pol.nvars();
int offset= 0;
_Ttype dum(0,2,_Ttype::dimen());
for(Tensor::index i= dum.begin();i!=dum.end();++i){
int d= i.getCoor().sum();
Symmetry symrun(_Ttype::dimen()-d,d);
_TGStype dumgs(0,TensorDimens(symrun,dumnvs));
if(pol.check(Symmetry(d))){
TwoDMatrix subt(*this,offset,dumgs.ncols());
subt.add(1.0,*(pol.get(Symmetry(d))));
}
offset+= dumgs.ncols();
}
}
/*:16*/
;
/*17:*/
void eval(Vector&out,const ConstVector&v)const
{
TL_RAISE_IF(v.length()+1!=_Ttype::nvar(),
"Wrong input vector length in CompactPolynomial::eval");
TL_RAISE_IF(out.length()!=_Ttype::nrows(),
"Wrong output vector length in CompactPolynomial::eval");
Vector x1(v.length()+1);
Vector x1p(x1,1,v.length());
x1p= v;
x1[0]= 1.0;
if(_Ttype::dimen()==0)
out= ConstVector(*this,0);
else{
PowerProvider pp(x1);
const _Stype&xpow= pp.getNext((const _Stype*)NULL);
for(int i= 1;i<_Ttype::dimen();i++)
xpow= pp.getNext((const _Stype*)NULL);
multVec(0.0,out,1.0,xpow);
}
}
/*:17*/
;
/*16:*/
CompactPolynomial(const TensorPolynomial<_Ttype,_TGStype,_Stype> &pol)
:_Ttype(pol.nrows(),pol.nvars()+1,pol.getMaxDim())
{
_Ttype::zeros();
IntSequence dumnvs(2);
dumnvs[0]= 1;
dumnvs[1]= pol.nvars();
int offset= 0;
_Ttype dum(0,2,_Ttype::dimen());
for(Tensor::index i= dum.begin();i!=dum.end();++i){
int d= i.getCoor().sum();
Symmetry symrun(_Ttype::dimen()-d,d);
_TGStype dumgs(0,TensorDimens(symrun,dumnvs));
if(pol.check(Symmetry(d))){
TwoDMatrix subt(*this,offset,dumgs.ncols());
subt.add(1.0,*(pol.get(Symmetry(d))));
}
offset+= dumgs.ncols();
}
}
/*:16*/
;
/*17:*/
void eval(Vector&out,const ConstVector&v)const
{
TL_RAISE_IF(v.length()+1!=_Ttype::nvar(),
"Wrong input vector length in CompactPolynomial::eval");
TL_RAISE_IF(out.length()!=_Ttype::nrows(),
"Wrong output vector length in CompactPolynomial::eval");
Vector x1(v.length()+1);
Vector x1p(x1,1,v.length());
x1p= v;
x1[0]= 1.0;
if(_Ttype::dimen()==0)
out= ConstVector(*this,0);
else{
PowerProvider pp(x1);
const _Stype&xpow= pp.getNext((const _Stype*)NULL);
for(int i= 1;i<_Ttype::dimen();i++)
xpow= pp.getNext((const _Stype*)NULL);
multVec(0.0,out,1.0,xpow);
}
}
/*:17*/
;
};
/*:15*/
@ -354,8 +354,8 @@ multVec(0.0,out,1.0,xpow);
class UCompactPolynomial:public CompactPolynomial<UFSTensor,UGSTensor,URSingleTensor> {
public:
UCompactPolynomial(const UTensorPolynomial&upol)
:CompactPolynomial<UFSTensor,UGSTensor,URSingleTensor> (upol){}
UCompactPolynomial(const UTensorPolynomial&upol)
:CompactPolynomial<UFSTensor,UGSTensor,URSingleTensor> (upol){}
};
/*:18*/
@ -364,8 +364,8 @@ UCompactPolynomial(const UTensorPolynomial&upol)
class FCompactPolynomial:public CompactPolynomial<FFSTensor,FGSTensor,FRSingleTensor> {
public:
FCompactPolynomial(const FTensorPolynomial&fpol)
:CompactPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fpol){}
FCompactPolynomial(const FTensorPolynomial&fpol)
:CompactPolynomial<FFSTensor,FGSTensor,FRSingleTensor> (fpol){}
};

View File

@ -10,36 +10,36 @@
/*2:*/
template<class _Tptr> class _index{
typedef _index<_Tptr> _Self;
_Tptr tensor;
int offset;
IntSequence coor;
typedef _index<_Tptr> _Self;
_Tptr tensor;
int offset;
IntSequence coor;
public:
_index(_Tptr t,int n)
:tensor(t),offset(0),coor(n,0){}
_index(_Tptr t,const IntSequence&cr,int c)
:tensor(t),offset(c),coor(cr){}
_index(_Tptr t,const IntSequence&cr)
:tensor(t),offset(tensor->getOffset(cr)),coor(cr){}
_index(const _index&ind)
:tensor(ind.tensor),offset(ind.offset),coor(ind.coor){}
const _Self&operator= (const _Self&in)
{tensor= in.tensor;offset= in.offset;coor= in.coor;
return*this;}
_Self&operator++()
{tensor->increment(coor);offset++;return*this;}
_Self&operator--()
{tensor->decrement(coor);offset--;return*this;}
int operator*()const
{return offset;}
bool operator==(const _index&n)const
{return offset==n.offset;}
bool operator!=(const _index&n)const
{return offset!=n.offset;}
const IntSequence&getCoor()const
{return coor;}
void print()const
{printf("%4d: ",offset);coor.print();}
_index(_Tptr t,int n)
:tensor(t),offset(0),coor(n,0){}
_index(_Tptr t,const IntSequence&cr,int c)
:tensor(t),offset(c),coor(cr){}
_index(_Tptr t,const IntSequence&cr)
:tensor(t),offset(tensor->getOffset(cr)),coor(cr){}
_index(const _index&ind)
:tensor(ind.tensor),offset(ind.offset),coor(ind.coor){}
const _Self&operator= (const _Self&in)
{tensor= in.tensor;offset= in.offset;coor= in.coor;
return*this;}
_Self&operator++()
{tensor->increment(coor);offset++;return*this;}
_Self&operator--()
{tensor->decrement(coor);offset--;return*this;}
int operator*()const
{return offset;}
bool operator==(const _index&n)const
{return offset==n.offset;}
bool operator!=(const _index&n)const
{return offset!=n.offset;}
const IntSequence&getCoor()const
{return coor;}
void print()const
{printf("%4d: ",offset);coor.print();}
};
/*:2*/
@ -48,55 +48,55 @@ void print()const
class Tensor:public TwoDMatrix{
public:
enum indor{along_row,along_col};
typedef _index<const Tensor*> index;
enum indor{along_row,along_col};
typedef _index<const Tensor*> index;
protected:
const index in_beg;
const index in_end;
int dim;
const index in_beg;
const index in_end;
int dim;
public:
Tensor(indor io,const IntSequence&last,int r,int c,int d)
:TwoDMatrix(r,c),
in_beg(this,d),
in_end(this,last,(io==along_row)?r:c),
dim(d){}
Tensor(indor io,const IntSequence&first,const IntSequence&last,
int r,int c,int d)
:TwoDMatrix(r,c),
in_beg(this,first,0),
in_end(this,last,(io==along_row)?r:c),
dim(d){}
Tensor(int first_row,int num,Tensor&t)
:TwoDMatrix(first_row,num,t),
in_beg(t.in_beg),
in_end(t.in_end),
dim(t.dim){}
Tensor(const Tensor&t)
:TwoDMatrix(t),
in_beg(this,t.in_beg.getCoor(),*(t.in_beg)),
in_end(this,t.in_end.getCoor(),*(t.in_end)),
dim(t.dim){}
virtual~Tensor(){}
virtual void increment(IntSequence&v)const= 0;
virtual void decrement(IntSequence&v)const= 0;
virtual int getOffset(const IntSequence&v)const= 0;
int dimen()const
{return dim;}
const index&begin()const
{return in_beg;}
const index&end()const
{return in_end;}
static int noverk(int n,int k);
static int power(int a,int b);
static int noverseq(const IntSequence&s)
{
IntSequence seq(s);
return noverseq_ip((IntSequence&)s);
}
Tensor(indor io,const IntSequence&last,int r,int c,int d)
:TwoDMatrix(r,c),
in_beg(this,d),
in_end(this,last,(io==along_row)?r:c),
dim(d){}
Tensor(indor io,const IntSequence&first,const IntSequence&last,
int r,int c,int d)
:TwoDMatrix(r,c),
in_beg(this,first,0),
in_end(this,last,(io==along_row)?r:c),
dim(d){}
Tensor(int first_row,int num,Tensor&t)
:TwoDMatrix(first_row,num,t),
in_beg(t.in_beg),
in_end(t.in_end),
dim(t.dim){}
Tensor(const Tensor&t)
:TwoDMatrix(t),
in_beg(this,t.in_beg.getCoor(),*(t.in_beg)),
in_end(this,t.in_end.getCoor(),*(t.in_end)),
dim(t.dim){}
virtual~Tensor(){}
virtual void increment(IntSequence&v)const= 0;
virtual void decrement(IntSequence&v)const= 0;
virtual int getOffset(const IntSequence&v)const= 0;
int dimen()const
{return dim;}
const index&begin()const
{return in_beg;}
const index&end()const
{return in_end;}
static int noverk(int n,int k);
static int power(int a,int b);
static int noverseq(const IntSequence&s)
{
IntSequence seq(s);
return noverseq_ip((IntSequence&)s);
}
private:
static int noverseq_ip(IntSequence&s);
static int noverseq_ip(IntSequence&s);
};
/*:3*/
@ -106,21 +106,21 @@ static int noverseq_ip(IntSequence&s);
class FTensor;
class UTensor:public Tensor{
public:
UTensor(indor io,const IntSequence&last,int r,int c,int d)
:Tensor(io,last,r,c,d){}
UTensor(const UTensor&ut)
:Tensor(ut){}
UTensor(int first_row,int num,UTensor&t)
:Tensor(first_row,num,t){}
virtual~UTensor(){}
virtual FTensor&fold()const= 0;
static void increment(IntSequence&v,int nv);
static void decrement(IntSequence&v,int nv);
static void increment(IntSequence&v,const IntSequence&nvmx);
static void decrement(IntSequence&v,const IntSequence&nvmx);
static int getOffset(const IntSequence&v,int nv);
static int getOffset(const IntSequence&v,const IntSequence&nvmx);
UTensor(indor io,const IntSequence&last,int r,int c,int d)
:Tensor(io,last,r,c,d){}
UTensor(const UTensor&ut)
:Tensor(ut){}
UTensor(int first_row,int num,UTensor&t)
:Tensor(first_row,num,t){}
virtual~UTensor(){}
virtual FTensor&fold()const= 0;
static void increment(IntSequence&v,int nv);
static void decrement(IntSequence&v,int nv);
static void increment(IntSequence&v,const IntSequence&nvmx);
static void decrement(IntSequence&v,const IntSequence&nvmx);
static int getOffset(const IntSequence&v,int nv);
static int getOffset(const IntSequence&v,const IntSequence&nvmx);
};
/*:4*/
@ -129,20 +129,20 @@ static int getOffset(const IntSequence&v,const IntSequence&nvmx);
class FTensor:public Tensor{
public:
FTensor(indor io,const IntSequence&last,int r,int c,int d)
:Tensor(io,last,r,c,d){}
FTensor(const FTensor&ft)
:Tensor(ft){}
FTensor(int first_row,int num,FTensor&t)
:Tensor(first_row,num,t){}
virtual~FTensor(){}
virtual UTensor&unfold()const= 0;
static void decrement(IntSequence&v,int nv);
static int getOffset(const IntSequence&v,int nv)
{IntSequence vtmp(v);return getOffsetRecurse(vtmp,nv);}
FTensor(indor io,const IntSequence&last,int r,int c,int d)
:Tensor(io,last,r,c,d){}
FTensor(const FTensor&ft)
:Tensor(ft){}
FTensor(int first_row,int num,FTensor&t)
:Tensor(first_row,num,t){}
virtual~FTensor(){}
virtual UTensor&unfold()const= 0;
static void decrement(IntSequence&v,int nv);
static int getOffset(const IntSequence&v,int nv)
{IntSequence vtmp(v);return getOffsetRecurse(vtmp,nv);}
private:
static int getOffsetRecurse(IntSequence&v,int nv);
static int getOffsetRecurse(IntSequence&v,int nv);
};
/*:5*/

View File

@ -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);}
};