stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Dynare++: improve comments using Unicode (first batch) 

[skip ci]
time-shift
Sébastien Villemot 2019-03-27 19:22:35 +01:00
parent aa16d46006
commit 5958848145
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
35 changed files with 307 additions and 303 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
dynare++/kord/faa_di_bruno.hh
View File

@ -1,8 +1,8 @@
// Copyright 2005, Ondra Kamenik
// Faa Di Bruno evaluator
// Faà Di Bruno evaluator
/* This defines a class which implements Faa Di Bruno Formula
/* This defines a class which implements Faà Di Bruno Formula
$$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_l}=\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
\sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s^k(c_m)}\right]^{\beta_m}_{c_m(\alpha)}$$
where $s^k$ is a general symmetry of dimension $k$ and $z$ is a stack of

14
dynare++/kord/korder.hh
View File

@ -214,9 +214,9 @@ public:
containers
- |sylvesterSolve|: solve the sylvester equation (templated fold, and
unfold)
- |faaDiBrunoZ|: calculates derivatives of $F$ by Faa Di Bruno for the
- |faaDiBrunoZ|: calculates derivatives of $F$ by Faà Di Bruno for the
sparse container of system derivatives and $Z$ stack container
- |faaDiBrunoG|: calculates derivatives of $G$ by Faa Di Bruno for the
- |faaDiBrunoG|: calculates derivatives of $G$ by Faà Di Bruno for the
dense container $g^{**}$ and $G$ stack
- |recover_y|: recovers $g_{y^{*i}}$
- |recover_yu|: recovers $g_{y^{*i}u^j}$
@ -387,7 +387,7 @@ KOrder::insertDerivative(std::unique_ptr<typename ctraits<t>::Ttensor> der)
ypart.nyss(), *der_ptr));
}
/* Here we implement Faa Di Bruno formula
/* Here we implement Faà Di Bruno formula
$$\sum_{l=1}^k\left[f_{z^l}\right]_{\gamma_1\ldots\gamma_l}
\sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s(c_m)}\right]^{\gamma_m},
$$
@ -399,7 +399,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
KOrder::faaDiBrunoZ(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno Z container for " << sym << endrec;
pa << "Faà Di Bruno Z container for " << sym << endrec;
auto res = std::make_unique<typename ctraits<t>::Ttensor>(ny, TensorDimens(sym, nvs));
FaaDiBruno bruno(journal);
bruno.calculate(Zstack<t>(), f, *res);
@ -414,7 +414,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
KOrder::faaDiBrunoG(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno G container for " << sym << endrec;
pa << "Faà Di Bruno G container for " << sym << endrec;
TensorDimens tdims(sym, nvs);
auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.nyss(), tdims);
FaaDiBruno bruno(journal);
@ -762,9 +762,9 @@ KOrder::calcE_k(int k) const
}
/* Here is the core routine. It calls methods recovering derivatives in
the right order. Recall, that the code, namely Faa Di Bruno's formula,
the right order. Recall, that the code, namely Faà Di Bruno's formula,
is implemented as to be run conditionally on the current contents of
containers. So, if some call of Faa Di Bruno evaluates derivatives,
containers. So, if some call of Faà Di Bruno evaluates derivatives,
and some derivatives are not present in the container, then it is
considered to be zero. So, we have to be very careful to put
everything in the right order. The order here can be derived from

4
dynare++/kord/korder_stoch.hh
View File

@ -479,7 +479,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno ZX container for " << sym << endrec;
pa << "Faà Di Bruno ZX container for " << sym << endrec;
auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.ny(), TensorDimens(sym, nvs));
FaaDiBruno bruno(journal);
bruno.calculate(Zstack<t>(), f, *res);
@ -494,7 +494,7 @@ std::unique_ptr<typename ctraits<t>::Ttensor>
KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno GX container for " << sym << endrec;
pa << "Faà Di Bruno GX container for " << sym << endrec;
TensorDimens tdims(sym, nvs);
auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.nyss(), tdims);
FaaDiBruno bruno(journal);

18
dynare++/sylv/cc/BlockDiagonal.cc
View File

@ -36,8 +36,8 @@ BlockDiagonal::BlockDiagonal(int p, const BlockDiagonal &b)
{
}
/* put zeroes to right upper submatrix whose first column is defined
* by 'edge' */
/* Put zeroes to right upper submatrix whose first column is defined
by ‘edge’ */
void
BlockDiagonal::setZerosToRU(diag_iter edge)
{
@ -47,14 +47,14 @@ BlockDiagonal::setZerosToRU(diag_iter edge)
get(i, j) = 0.0;
}
/* Updates row_len and col_len so that there are zeroes in upper right part, this
* |T1 0 |
* |0 T2|. The first column of T2 is given by diagonal iterator 'edge'.
/* Updates row_len and col_len so that there are zeroes in upper right part, i.e.
⎡ T1 0 ⎤
⎣ 0 T2 ⎦. The first column of T2 is given by diagonal iterator ‘edge’.
* Note the semantics of row_len and col_len. row_len[i] is distance
* of the right-most non-zero element of i-th row from the left, and
* col_len[j] is distance of top-most non-zero element of j-th column
* to the top. (First element has distance 1).
Note the semantics of row_len and col_len. row_len[i] is distance
of the right-most non-zero element of i-th row from the left, and
col_len[j] is distance of top-most non-zero element of j-th column
to the top. (First element has distance 1).
*/
void
BlockDiagonal::setZeroBlockEdge(diag_iter edge)

8
dynare++/sylv/cc/GeneralMatrix.cc
View File

@ -329,13 +329,13 @@ GeneralMatrix::gemm_partial_right(const std::string &trans, const ConstGeneralMa
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols)
: data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
{
// can check that the submatirx is fully in the matrix
// FIXME: check that the submatrix is fully in the matrix
}
ConstGeneralMatrix::ConstGeneralMatrix(const ConstGeneralMatrix &m, int i, int j, int nrows, int ncols)
: data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
{
// can check that the submatirx is fully in the matrix
// FIXME: check that the submatrix is fully in the matrix
}
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix &m)
@ -559,8 +559,8 @@ SVDDecomp::solve(const ConstGeneralMatrix &B, GeneralMatrix &X) const
if (B.numRows() != U.numRows())
throw SYLV_MES_EXCEPTION("Incompatible number of rows ");
// reciprocal condition number for determination of zeros in the
// end of sigma
/* Reciprocal condition number for determination of zeros in the
end of sigma */
constexpr double rcond = 1e-13;
// determine nz=number of zeros in the end of sigma

67
dynare++/sylv/cc/GeneralMatrix.hh
View File

@ -282,7 +282,7 @@ public:
place(ConstGeneralMatrix(m), i, j);
}
/* this = a*b */
// this = a·b
void mult(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b);
void
mult(const GeneralMatrix &a, const GeneralMatrix &b)
@ -290,7 +290,7 @@ public:
mult(ConstGeneralMatrix(a), ConstGeneralMatrix(b));
}
/* this = this + scalar*a*b */
// this = this + scalar·a·b
void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
double mult = 1.0);
void
@ -300,7 +300,7 @@ public:
multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), mult);
}
/* this = this + scalar*a*b' */
// this = this + scalar·a·bᵀ
void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b,
const std::string &dum, double mult = 1.0);
void
@ -310,7 +310,7 @@ public:
multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult);
}
/* this = this + scalar*a'*b */
// this = this + scalar·aᵀ·b
void multAndAdd(const ConstGeneralMatrix &a, const std::string &dum, const ConstGeneralMatrix &b,
double mult = 1.0);
void
@ -320,7 +320,7 @@ public:
multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult);
}
/* this = this + scalar*a'*b' */
// this = this + scalar·aᵀ·bᵀ
void multAndAdd(const ConstGeneralMatrix &a, const std::string &dum1,
const ConstGeneralMatrix &b, const std::string &dum2, double mult = 1.0);
void
@ -330,7 +330,7 @@ public:
multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b), dum2, mult);
}
/* this = this + scalar*a*a' */
// this = this + scalar·a·aᵀ
void addOuter(const ConstVector &a, double mult = 1.0);
void
addOuter(const Vector &a, double mult = 1.0)
@ -338,7 +338,7 @@ public:
addOuter(ConstVector(a), mult);
}
/* this = this * m */
// this = this·m
void multRight(const ConstGeneralMatrix &m);
void
multRight(const GeneralMatrix &m)
@ -346,7 +346,7 @@ public:
multRight(ConstGeneralMatrix(m));
}
/* this = m * this */
// this = m·this
void multLeft(const ConstGeneralMatrix &m);
void
multLeft(const GeneralMatrix &m)
@ -354,7 +354,7 @@ public:
multLeft(ConstGeneralMatrix(m));
}
/* this = this * m' */
// this = this·mᵀ
void multRightTrans(const ConstGeneralMatrix &m);
void
multRightTrans(const GeneralMatrix &m)
@ -362,7 +362,7 @@ public:
multRightTrans(ConstGeneralMatrix(m));
}
/* this = m' * this */
// this = mᵀ·this
void multLeftTrans(const ConstGeneralMatrix &m);
void
multLeftTrans(const GeneralMatrix &m)
@ -370,64 +370,64 @@ public:
multLeftTrans(ConstGeneralMatrix(m));
}
/* x = scalar(a)*x + scalar(b)*this*d */
// x = scalar(a)·x + scalar(b)·this·d
void
multVec(double a, Vector &x, double b, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multVec(a, x, b, d);
}
/* x = scalar(a)*x + scalar(b)*this'*d */
// x = scalar(a)·x + scalar(b)·thisᵀ·d
void
multVecTrans(double a, Vector &x, double b, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multVecTrans(a, x, b, d);
}
/* x = x + this*d */
// x = x + this·d
void
multaVec(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multaVec(x, d);
}
/* x = x + this'*d */
// x = x + thisᵀ·d */
void
multaVecTrans(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multaVecTrans(x, d);
}
/* x = x - this*d */
// x = x - this·d
void
multsVec(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multsVec(x, d);
}
/* x = x - this'*d */
// x = x - thisᵀ·d
void
multsVecTrans(Vector &x, const ConstVector &d) const
{
ConstGeneralMatrix(*this).multsVecTrans(x, d);
}
/* this = zero */
// this = zero
void zeros();
/** this = unit (on main diagonal) */
// this = unit (on main diagonal)
void unit();
/* this = NaN */
// this = NaN
void nans();
/* this = Inf */
// this = ∞
void infs();
/* this = scalar*this */
// this = scalar·this
void mult(double a);
/* this = this + scalar*m */
// this = this + scalar·m
void add(double a, const ConstGeneralMatrix &m);
void
add(double a, const GeneralMatrix &m)
@ -435,7 +435,7 @@ public:
add(a, ConstGeneralMatrix(m));
}
/* this = this + scalar*m' */
// this = this + scalar·mᵀ
void add(double a, const ConstGeneralMatrix &m, const std::string &dum);
void
add(double a, const GeneralMatrix &m, const std::string &dum)
@ -490,7 +490,7 @@ private:
gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta);
}
/* this = op(m) *this (without whole copy of this) */
// this = op(m)·this (without whole copy of ‘this’)
void gemm_partial_left(const std::string &trans, const ConstGeneralMatrix &m,
double alpha, double beta);
void
@ -504,7 +504,7 @@ private:
static constexpr int md_length = 23;
};
// Computes a*b
// Computes a·b
inline GeneralMatrix
operator*(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b)
{
@ -513,7 +513,7 @@ operator*(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b)
return m;
}
// Computes a*b'
// Computes a·bᵀ
template<class T>
GeneralMatrix
operator*(const ConstGeneralMatrix &a, const TransposedMatrix<T> &b)
@ -523,7 +523,7 @@ operator*(const ConstGeneralMatrix &a, const TransposedMatrix<T> &b)
return m;
}
// Computes a'*b
// Computes aᵀ·b
template<class T>
GeneralMatrix
operator*(const TransposedMatrix<T> &a, const ConstGeneralMatrix &b)
@ -533,7 +533,7 @@ operator*(const TransposedMatrix<T> &a, const ConstGeneralMatrix &b)
return m;
}
// Computes a'*b'
// Computes aᵀ·bᵀ
template<class T1, class T2>
GeneralMatrix
operator*(const TransposedMatrix<T1> &a, const TransposedMatrix<T2> &b)
@ -546,16 +546,15 @@ operator*(const TransposedMatrix<T1> &a, const TransposedMatrix<T2> &b)
class SVDDecomp
{
protected:
/** Minimum of number of rows and columns of the decomposed
* matrix. */
// Minimum of number of rows and columns of the decomposed matrix
const int minmn;
/** Singular values. */
// Singular values
Vector sigma;
/** Orthogonal matrix U. */
// Orthogonal matrix U
GeneralMatrix U;
/** Orthogonal matrix V^T. */
// Orthogonal matrix Vᵀ
GeneralMatrix VT;
/** Convered flag. */
// Convered flag
bool conv;
public:
SVDDecomp(const GeneralMatrix &A)

2
dynare++/sylv/cc/GeneralSylvester.cc
View File

@ -107,7 +107,7 @@ GeneralSylvester::check(const ConstVector &ds)
if (!solved)
throw SYLV_MES_EXCEPTION("Cannot run check on system, which is not solved yet.");
// calculate xcheck = AX+BXC^i-D
// calculate xcheck = A·X+B·X·⊗ⁱC-D
SylvMatrix dcheck(d.numRows(), d.numCols());
dcheck.multLeft(b.numRows()-b.numCols(), b, d);
dcheck.multRightKron(c, order);

4
dynare++/sylv/cc/GeneralSylvester.hh
View File

@ -24,7 +24,7 @@ class GeneralSylvester
std::unique_ptr<SimilarityDecomp> cdecomp;
std::unique_ptr<SylvesterSolver> sylv;
public:
/* construct with my copy of d*/
// Construct with my copy of d
GeneralSylvester(int ord, int n, int m, int zero_cols,
const ConstVector &da, const ConstVector &db,
const ConstVector &dc, const ConstVector &dd,
@ -33,7 +33,7 @@ public:
const ConstVector &da, const ConstVector &db,
const ConstVector &dc, const ConstVector &dd,
bool alloc_for_check = false);
/* construct with provided storage for d */
// Construct with provided storage for d
GeneralSylvester(int ord, int n, int m, int zero_cols,
const ConstVector &da, const ConstVector &db,
const ConstVector &dc, Vector &dd,

18
dynare++/sylv/cc/KronUtils.cc
View File

@ -10,13 +10,11 @@ KronUtils::multAtLevel(int level, const QuasiTriangular &t,
KronVector &x)
{
if (0 < level && level < x.getDepth())
{
for (int i = 0; i < x.getM(); i++)
{
KronVector xi(x, i);
multAtLevel(level, t, xi);
}
}
for (int i = 0; i < x.getM(); i++)
{
KronVector xi(x, i);
multAtLevel(level, t, xi);
}
else if (0 == level && 0 < x.getDepth())
{
GeneralMatrix tmp(x, x.getN(), power(x.getM(), x.getDepth()));
@ -61,8 +59,6 @@ KronUtils::multKron(const QuasiTriangular &f, const QuasiTriangular &k,
{
multAtLevel(0, k, x);
if (x.getDepth() > 0)
{
for (int level = 1; level <= x.getDepth(); level++)
multAtLevelTrans(level, f, x);
}
for (int level = 1; level <= x.getDepth(); level++)
multAtLevelTrans(level, f, x);
}

11
dynare++/sylv/cc/KronUtils.hh
View File

@ -11,17 +11,16 @@
class KronUtils
{
public:
/* multiplies I_m\otimes..\I_m\otimes T\otimes I_m...I_m\otimes I_n
with given b and returns x. T must be (m,m), number of
\otimes is b.getDepth(), level is a number of I_m's between T
and I_n plus 1. If level=0, then we multiply
\I_m\otimes ..\otimes I_m\otimes T, T is (n,n) */
/* Computes x = (Iₘ⊗…⊗Iₘ⊗T⊗Iₘ⊗…⊗Iₘ⊗Iₙ)·x, where x is n×mᵈ.
T must be m×m, number of ⊗ is d, level is the number of Iₘ’s
between T and Iₙ plus 1. If level=0, then we multiply by Iₘ⊗…⊗Iₘ⊗T,
T must be n×n. */
static void multAtLevel(int level, const QuasiTriangular &t,
KronVector &x);
static void multAtLevelTrans(int level, const QuasiTriangular &t,
KronVector &x);
/* multiplies x=(F'\otimes F'\otimes..\otimes K)x */
// Computes x=(Fᵀ⊗Fᵀ⊗…⊗K)·x
static void multKron(const QuasiTriangular &f, const QuasiTriangular &k,
KronVector &x);
};

22
dynare++/sylv/cc/QuasiTriangular.cc
View File

@ -34,9 +34,9 @@ DiagonalBlock::getSize() const
return std::sqrt(getDeterminant());
}
// this function makes Diagonal inconsistent, it should only be used
// on temorary matrices, which will not be used any more, e.g. in
// QuasiTriangular::solve (we need fast performance)
/* This function makes Diagonal inconsistent, it should only be used
on temorary matrices, which will not be used any more, e.g. in
QuasiTriangular::solve() (we need fast performance) */
void
DiagonalBlock::setReal()
{
@ -64,7 +64,7 @@ DiagonalBlock::checkBlock(const double *d, int d_size)
Diagonal::Diagonal(double *data, int d_size)
{
int nc = getNumComplex(data, d_size); // return nc <= d_size/2
int nc = getNumComplex(data, d_size); // return nc ≤ d_size/2
num_all = d_size - nc;
num_real = d_size - 2*nc;
@ -184,9 +184,9 @@ Diagonal::getEigenValues(Vector &eig) const
}
}
/* swaps logically blocks 'it', and '++it'. remember to move also
* addresses, alpha, beta1, beta2. This is a dirty (but most
* effective) way how to do it. */
/* Swaps logically blocks ‘it’, and ‘++it’. remember to move also
addresses, alpha, beta1, beta2. This is a dirty (but most
effective) way how to do it. */
void
Diagonal::swapLogically(diag_iter it)
{
@ -432,7 +432,7 @@ QuasiTriangular::QuasiTriangular(const SchurDecomp &decomp)
{
}
/* this pads matrix with intial columns with zeros */
// This pads matrix with intial columns with zeros
QuasiTriangular::QuasiTriangular(const SchurDecompZero &decomp)
: SqSylvMatrix(decomp.getDim())
{
@ -590,7 +590,7 @@ QuasiTriangular::solvePreTrans(Vector &x, double &eig_min)
dtrsv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
}
/* calculates x = Tb */
// Calculates x = T·b
void
QuasiTriangular::multVec(Vector &x, const ConstVector &b) const
{
@ -639,7 +639,7 @@ QuasiTriangular::multaVecTrans(Vector &x, const ConstVector &b) const
x.add(1.0, tmp);
}
/* calculates x=x+(T\otimes I)b, where size of I is given by b (KronVector) */
// Calculates x=x+(T⊗I)·b, where size of I is given by b (KronVector)
void
QuasiTriangular::multaKron(KronVector &x, const ConstKronVector &b) const
{
@ -649,7 +649,7 @@ QuasiTriangular::multaKron(KronVector &x, const ConstKronVector &b) const
x_resh.multAndAdd(b_resh, ConstGeneralMatrix(*this), "trans");
}
/* calculates x=x+(T'\otimes I)b, where size of I is given by b (KronVector) */
// Calculates x=x+(T'⊗I)·b, where size of I is given by b (KronVector)
void
QuasiTriangular::multaKronTrans(KronVector &x, const ConstKronVector &b) const
{

24
dynare++/sylv/cc/QuasiTriangular.hh
View File

@ -335,29 +335,29 @@ public:
diag_iter findClosestDiagBlock(diag_iter start, diag_iter end, double a);
diag_iter findNextLargerBlock(diag_iter start, diag_iter end, double a);
/* (I+T)y = x, y-->x */
/* (I+T)·y = x, y→x */
virtual void solvePre(Vector &x, double &eig_min);
/* (I+T')y = x, y-->x */
/* (I+T')·y = x, y→x */
virtual void solvePreTrans(Vector &x, double &eig_min);
/* (I+T)x = b */
virtual void solve(Vector &x, const ConstVector &b, double &eig_min);
/* (I+T')x = b */
virtual void solveTrans(Vector &x, const ConstVector &b, double &eig_min);
/* x = Tb */
/* x = T·b */
virtual void multVec(Vector &x, const ConstVector &b) const;
/* x = T'b */
/* x = T'·b */
virtual void multVecTrans(Vector &x, const ConstVector &b) const;
/* x = x + Tb */
/* x = x + T·b */
virtual void multaVec(Vector &x, const ConstVector &b) const;
/* x = x + T'b */
/* x = x + T'·b */
virtual void multaVecTrans(Vector &x, const ConstVector &b) const;
/* x = (T\otimes I)x */
/* x = (T⊗I)·x */
virtual void multKron(KronVector &x) const;
/* x = (T'\otimes I)x */
/* x = (T'⊗I)·x */
virtual void multKronTrans(KronVector &x) const;
/* A = T*A */
/* A = T·A */
virtual void multLeftOther(GeneralMatrix &a) const;
/* A = T'*A */
/* A = T'·A */
virtual void multLeftOtherTrans(GeneralMatrix &a) const;
const_diag_iter
@ -417,9 +417,9 @@ protected:
void addMatrix(double r, const QuasiTriangular &t);
private:
void addUnit();
/* x = x + (T\otimes I)b */
/* x = x + (T⊗I)·b */
void multaKron(KronVector &x, const ConstKronVector &b) const;
/* x = x + (T'\otimes I)b */
/* x = x + (T'⊗)·b */
void multaKronTrans(KronVector &x, const ConstKronVector &b) const;
/* implementation via iterators, useful for large matrices */
void setMatrixViaIter(double r, const QuasiTriangular &t);

2
dynare++/sylv/cc/SchurDecomp.hh
View File

@ -47,7 +47,7 @@ public:
class SchurDecompZero : public SchurDecomp
{
GeneralMatrix ru; /* right upper matrix */
GeneralMatrix ru; // right upper matrix
public:
SchurDecompZero(const GeneralMatrix &m);
ConstGeneralMatrix

40
dynare++/sylv/cc/SchurDecompEig.cc
View File

@ -5,11 +5,11 @@
#include <vector>
/* bubble diagonal 1-1, or 2-2 block from position 'from' to position
* 'to'. If an eigenvalue cannot be swapped with its neighbour, the
* neighbour is bubbled also in front. The method returns a new
* position 'to', where the original block pointed by 'to' happens to
* appear at the end. 'from' must be greater than 'to'.
/* Bubble diagonal 1-1, or 2-2 block from position ‘from’ to position
‘to’. If an eigenvalue cannot be swapped with its neighbour, the
neighbour is bubbled also in front. The method returns a new
position ‘to’, where the original block pointed by ‘to’ happens to
appear at the end. ‘from’ must be greater than ‘to’.
*/
SchurDecompEig::diag_iter
SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
@ -22,9 +22,9 @@ SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
++to;
else
{
// bubble all eigenvalues from runm(incl.) to run(excl.),
// this includes either bubbling generated eigenvalues due
// to split, or an eigenvalue which couldn't be swapped
/* Bubble all eigenvalues from runm(incl.) to run(excl.),
this includes either bubbling generated eigenvalues due
to split, or an eigenvalue which couldn't be swapped */
while (runm != run)
{
to = bubbleEigen(runm, to);
@ -35,16 +35,16 @@ SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
return to;
}
/* this tries to swap two neighbouring eigenvalues, 'it' and '--it',
* and returns 'itadd'. If the blocks can be swapped, new eigenvalues
* can emerge due to possible 2-2 block splits. 'it' then points to
* the last eigenvalue coming from block pointed by 'it' at the
* begining, and 'itadd' points to the first. On swap failure, 'it' is
* not changed, and 'itadd' points to previous eignevalue (which must
* be moved backwards before). In either case, it is necessary to
* resolve eigenvalues from 'itadd' to 'it', before the 'it' can be
* resolved.
* The success is signaled by returned true.
/* This tries to swap two neighbouring eigenvalues, ‘it’ and ‘--it’,
and returns ‘itadd’. If the blocks can be swapped, new eigenvalues
can emerge due to possible 2-2 block splits. ‘it’ then points to
the last eigenvalue coming from block pointed by ‘it’ at the
begining, and ‘itadd’ points to the first. On swap failure, ‘it’ is
not changed, and ‘itadd’ points to previous eignevalue (which must
be moved backwards before). In either case, it is necessary to
resolve eigenvalues from ‘itadd’ to ‘it’, before the ‘it’ can be
resolved.
The success is signaled by returned true.
*/
bool
SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
@ -66,10 +66,10 @@ SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
{
// swap successful
getT().swapDiagLogically(itadd);
//check for 2-2 block splits
// check for 2-2 block splits
getT().checkDiagConsistency(it);
getT().checkDiagConsistency(itadd);
// and go back by 'it' in NEW eigenvalue set
// and go back by ‘it’ in NEW eigenvalue set
--it;
return true;
}

15
dynare++/sylv/cc/SimilarityDecomp.cc
View File

@ -33,9 +33,9 @@ SimilarityDecomp::getXDim(diag_iter start, diag_iter end,
rows = ei - si;
}
/* find solution of X for diagonal block given by start(incl.) and
* end(excl.). If the solution cannot be found, or it is greater than
* norm, X is not changed and flase is returned.
/* Find solution of X for diagonal block given by start(incl.) and
end(excl.). If the solution cannot be found, or it is greater than
norm, X is not changed and flase is returned.
*/
bool
SimilarityDecomp::solveX(diag_iter start, diag_iter end,
@ -68,7 +68,8 @@ SimilarityDecomp::solveX(diag_iter start, diag_iter end,
return true;
}
/* multiply Q and invQ with (I -X; 0 I), and (I X; 0 I). This also sets X=-X. */
/* ⎡ I -X ⎤ ⎡ I X ⎤
Multiply Q and invQ with ⎣ 0 I ⎦ and ⎣ 0 I ⎦. This also sets X=-X. */
void
SimilarityDecomp::updateTransform(diag_iter start, diag_iter end,
GeneralMatrix &X)
@ -127,21 +128,21 @@ SimilarityDecomp::diagonalize(double norm)
void
SimilarityDecomp::check(SylvParams &pars, const GeneralMatrix &m) const
{
// M - Q*B*inv(Q)
// M - Q·B·Q⁻¹
SqSylvMatrix c(getQ() * getB());
c.multRight(getInvQ());
c.add(-1.0, m);
pars.f_err1 = c.getNorm1();
pars.f_errI = c.getNormInf();
// I - Q*inv(Q)
// I - Q·Q⁻¹
c.setUnit();
c.mult(-1);
c.multAndAdd(getQ(), getInvQ());
pars.viv_err1 = c.getNorm1();
pars.viv_errI = c.getNormInf();
// I - inv(Q)*Q
// I - Q⁻¹·Q
c.setUnit();
c.mult(-1);
c.multAndAdd(getInvQ(), getQ());

6
dynare++/sylv/cc/SylvMatrix.cc
View File

@ -42,9 +42,9 @@ SylvMatrix::multLeft(int zero_cols, const GeneralMatrix &a, const GeneralMatrix
|| rows != b.numRows() || cols != b.numCols())
throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeft.");
// here we cannot call SylvMatrix::gemm since it would require
// another copy of (usually big) b (we are not able to do inplace
// submatrix of const GeneralMatrix)
/* Here we cannot call SylvMatrix::gemm() since it would require
another copy of (usually big) b (we are not able to do inplace
submatrix of const GeneralMatrix) */
if (a.getLD() > 0 && ld > 0)
{
blas_int mm = a.numRows();

30
dynare++/sylv/cc/SylvMatrix.hh
View File

@ -40,23 +40,23 @@ public:
SylvMatrix &operator=(const SylvMatrix &m) = default;
SylvMatrix &operator=(SylvMatrix &&m) = default;
/* this = |I 0|* this
|0 m| */
/* ⎡ I 0 ⎤
this = ⎣ 0 m ⎦ · this */
void multLeftI(const SqSylvMatrix &m);
/* this = |I 0|* this
|0 m'| */
/* ⎡ I 0 ⎤
this = ⎣ 0 m' ⎦ · this */
void multLeftITrans(const SqSylvMatrix &m);
/* this = |0 a|*b, so that |0 a| is square */
// this = [0 a]·b, so that [0 a] is square
void multLeft(int zero_cols, const GeneralMatrix &a, const GeneralMatrix &b);
/* this = this * (m\otimes m..\otimes m) */
// this = this·(m⊗m⊗…⊗m)
void multRightKron(const SqSylvMatrix &m, int order);
/* this = this * (m'\otimes m'..\otimes m') */
// this = this·(mᵀ⊗mᵀ⊗…⊗mᵀ)
void multRightKronTrans(const SqSylvMatrix &m, int order);
/* this = P*this, x = P*x, where P is gauss transformation setting
* a given element to zero */
/* this = P·this, x = P·x, where P is gauss transformation setting
a given element to zero */
void eliminateLeft(int row, int col, Vector &x);
/* this = this*P, x = P'*x, where P is gauss transformation setting
* a given element to zero */
/* this = this·P, x = Pᵀ·x, where P is gauss transformation setting
a given element to zero */
void eliminateRight(int row, int col, Vector &x);
};
@ -84,14 +84,14 @@ public:
}
SqSylvMatrix &operator=(const SqSylvMatrix &m) = default;
SqSylvMatrix &operator=(SqSylvMatrix &&m) = default;
/* x = (this \otimes this..\otimes this)*d */
// x = (this⊗this⊗…⊗this)·d
void multVecKron(KronVector &x, const ConstKronVector &d) const;
/* x = (this' \otimes this'..\otimes this')*d */
// x = (thisᵀ⊗thisᵀ⊗…⊗thisᵀ)·d
void multVecKronTrans(KronVector &x, const ConstKronVector &d) const;
/* a = inv(this)*a, b=inv(this)*b */
// a = this⁻¹·a, b=this⁻¹·b */
void multInvLeft2(GeneralMatrix &a, GeneralMatrix &b,
double &rcond1, double &rcondinf) const;
/* this = I */
// this = I
void setUnit();
};

4
dynare++/sylv/cc/SylvesterSolver.hh
View File

@ -20,8 +20,8 @@ protected:
const std::unique_ptr<const QuasiTriangular> matrixK;
const std::unique_ptr<const QuasiTriangular> matrixF;
private:
/* return true when it is more efficient to use QuasiTriangular
* than QuasiTriangularZero */
/* Return true when it is more efficient to use QuasiTriangular
than QuasiTriangularZero */
static bool
zeroPad(const SchurDecompZero &kdecomp)
{

4
dynare++/sylv/cc/SymSchurDecomp.cc
View File

@ -78,8 +78,8 @@ SymSchurDecomp::getFactor(GeneralMatrix &f) const
}
}
// LAPACK says that eigenvalues are ordered in ascending order, but we
// do not rely her on it
/* LAPACK says that eigenvalues are ordered in ascending order, but we
do not rely on it */
bool
SymSchurDecomp::isPositiveSemidefinite() const
{

15
dynare++/sylv/cc/SymSchurDecomp.hh
View File

@ -13,8 +13,8 @@ protected:
Vector lambda;
SqSylvMatrix q;
public:
/** Calculates A = Q*Lambda*Q^T, where A is assummed to be
* symmetric and Lambda real diagonal, hence a vector. */
/* Computes the factorization A = Q·Λ·Qᵀ, where A is assummed to be
symmetric and Λ real diagonal, hence a vector. */
SymSchurDecomp(const ConstGeneralMatrix &a);
SymSchurDecomp(const SymSchurDecomp &ssd) = default;
virtual ~SymSchurDecomp() = default;
@ -28,15 +28,14 @@ public:
{
return q;
}
/** Return factor F*F^T = A, raises and exception if A is not
* positive semidefinite, F must be square. */
/* Return factor F·Fᵀ = A, raises and exception if A is not
positive semidefinite, F must be square. */
void getFactor(GeneralMatrix &f) const;
/** Returns true if A is positive semidefinite. */
// Returns true if A is positive semidefinite.
bool isPositiveSemidefinite() const;
/** Correct definitness. This sets all eigenvalues between minus
* tolerance and zero to zero. */
/* Correct definitness. This sets all eigenvalues between minus
tolerance and zero to zero. */
void correctDefinitness(double tol);
};
#endif

18
dynare++/sylv/cc/Vector.hh
View File

@ -5,9 +5,9 @@
#ifndef VECTOR_H
#define VECTOR_H
/* NOTE! Vector and ConstVector have not common super class in order
* to avoid running virtual method invokation mechanism. Some
* members, and methods are thus duplicated */
/* NOTE: Vector and ConstVector have not common super class in order
to avoid running virtual method invokation mechanism. Some
members, and methods are thus duplicated */
#include <complex>
@ -23,7 +23,7 @@ class Vector
friend class ConstVector;
protected:
int len{0};
int s{1}; // stride (also called "skip" in some places)
int s{1}; // stride (also called “skip” in some places)
double *data;
bool destroy{true};
public:
@ -43,7 +43,7 @@ public:
v.data = nullptr;
v.destroy = false;
}
// We don't want implict conversion from ConstVector, since it's expensive
// We don't want implict conversion from ConstVector, since it’s expensive
explicit Vector(const ConstVector &v);
Vector(double *d, int l)
: len(l), data{d}, destroy{false}
@ -90,10 +90,10 @@ public:
return s;
}
/** Exact equality. */
// Exact equality.
bool operator==(const Vector &y) const;
bool operator!=(const Vector &y) const;
/** Lexicographic ordering. */
// Lexicographic ordering.
bool operator<(const Vector &y) const;
bool operator<=(const Vector &y) const;
bool operator>(const Vector &y) const;
@ -125,8 +125,8 @@ public:
void print() const;
/* multiplies | alpha -beta1| |b1| |x1|
| |\otimes I .| | = | |
| -beta2 alpha| |b2| |x2|
| |\otimes I .| | = | |
| -beta2 alpha| |b2| |x2|
*/
static void mult2(double alpha, double beta1, double beta2,
Vector &x1, Vector &x2,

44
dynare++/tl/cc/equivalence.cc
View File

@ -39,7 +39,7 @@ OrdSequence::operator==(const OrdSequence &s) const
return (i == length());
}
/* The first |add| adds a given integer to the class, the second
/* The first add() adds a given integer to the class, the second
iterates through a given sequence and adds everything found in the
given class. */
@ -65,7 +65,7 @@ OrdSequence::add(const OrdSequence &s)
}
}
/* Answers |true| if a given number is in the class. */
/* Answers true if a given number is in the class. */
bool
OrdSequence::has(int i) const
{
@ -183,7 +183,7 @@ Equivalence::findHaving(int i)
return si;
}
/* Find $j$-th class for a given $j$. */
/* Find j-th class for a given j. */
Equivalence::const_seqit
Equivalence::find(int j) const
@ -228,7 +228,7 @@ Equivalence::insert(const OrdSequence &s)
/* Trace the equivalence into the integer sequence. The classes are in
some order (described earlier), and items within classes are ordered,
so this implies, that the data can be linearized. This method
``prints'' them to the sequence. We allow for tracing only a given
“prints” them to the sequence. We allow for tracing only a given
number of classes from the beginning. */
void
@ -275,18 +275,16 @@ Equivalence::print(const std::string &prefix) const
}
}
/* Here we construct a set of all equivalences over $n$-element
set. The construction proceeds as follows. We maintain a list of added
equivalences. At each iteration we pop front of the list, try to add
all parents of the popped equivalence. This action adds new
equivalences to the object and also to the added list. We finish the
iterations when the added list is empty.
/* Here we construct a set of all equivalences over n-element set. The
construction proceeds as follows. We maintain a list of added equivalences.
At each iteration we pop front of the list, try to add all parents of the
popped equivalence. This action adds new equivalences to the object and also
to the added list. We finish the iterations when the added list is empty.
In the beginning we start with
$\{\{0\},\{1\},\ldots,\{n-1\}\}$. Adding of parents is an action which
for a given equivalence tries to glue all possible couples and checks
whether a new equivalence is already in the equivalence set. This is
not effective, but we will do the construction only ones.
In the beginning we start with { {0}, {1}, …, {n-1} }. Adding of parents is
an action which for a given equivalence tries to glue all possible couples
and checks whether a new equivalence is already in the equivalence set. This
is not effective, but we will do the construction only ones.
In this way we breath-first search a lattice of all equivalences. Note
that the lattice is modular, that is why the result of a construction
@ -311,12 +309,12 @@ EquivalenceSet::EquivalenceSet(int num)
equis.emplace_back(n, "");
}
/* This method is used in |addParents| and returns |true| if the object
/* This method is used in addParents() and returns true if the object
already has that equivalence. We trace list of equivalences in reverse
order since equivalences are ordered in the list from the most
primitive (nothing equivalent) to maximal (all is equivalent). Since
we will have much more results of |has| method as |true|, and
|operator==| between equivalences is quick if number of classes
we will have much more results of has() method as true, and
operator==() between equivalences is quick if number of classes
differ, and in time we will compare with equivalences with less
classes, then it is more efficient to trace the equivalences from less
classes to more classes. hence the reverse order. */
@ -334,7 +332,7 @@ EquivalenceSet::has(const Equivalence &e) const
/* Responsibility of this methods is to try to glue all possible
couples within a given equivalence and add those which are not in the
list yet. These are added also to the |added| list.
list yet. These are added also to the ‘added’ list.
If number of classes is 2 or 1, we exit, because there is nothing to
be added. */
@ -373,14 +371,14 @@ EquivalenceSet::print(const std::string &prefix) const
}
}
/* Construct the bundle. |nmax| is a maximum size of underlying set. */
/* Construct the bundle. nmax is a maximum size of underlying set. */
EquivalenceBundle::EquivalenceBundle(int nmax)
{
nmax = std::max(nmax, 1);
generateUpTo(nmax);
}
/* Remember, that the first item is |EquivalenceSet(1)|. */
/* Remember, that the first item is EquivalenceSet(1). */
const EquivalenceSet &
EquivalenceBundle::get(int n) const
{
@ -389,8 +387,8 @@ EquivalenceBundle::get(int n) const
return bundle[n-1];
}
/* Get |curmax| which is a maximum size in the bundle, and generate for
all sizes from |curmax+1| up to |nmax|. */
/* Get ‘curmax’ which is a maximum size in the bundle, and generate for
all sizes from curmax+1 up to nmax. */
void
EquivalenceBundle::generateUpTo(int nmax)

62
dynare++/tl/cc/equivalence.hh
View File

@ -2,34 +2,30 @@
// Equivalences.
/* Here we define an equivalence of a set of integers $\{0, 1, \ldots,
k-1\}$. The purpose is clear, in the tensor library we often iterate
/* Here we define an equivalence of a set of integers {0, 1, …, k-1}.
The purpose is clear, in the tensor library we often iterate
through all equivalences and sum matrices. We need an abstraction for
an equivalence class, equivalence and a set of all equivalences.
The equivalence class (which is basically a set of integers) is here
implemented as ordered integer sequence. The ordered sequence is not
implemented via |IntSequence|, but via |vector<int>| since we need
implemented via IntSequence, but via vector<int> since we need
insertions. The equivalence is implemented as an ordered list of
equivalence classes, and equivalence set is a list of equivalences.
The ordering of the equivalence classes within an equivalence is very
important. For instance, if we iterate through equivalences for $k=5$
and pickup some equivalence class, say $\{\{0,4\},\{1,2\},\{3\}\}$, we
important. For instance, if we iterate through equivalences for k=5
and pickup some equivalence class, say { {0,4}, {1,2}, {3} }, we
then evaluate something like:
$$\left[B_{y^2u^3}\right]_{\alpha_1\alpha_2\beta_1\beta_2\beta_3}=
\cdots+\left[g_{y^3}\right]_{\gamma_1\gamma_2\gamma_3}
\left[g_{yu}\right]^{\gamma_1}_{\alpha_1\beta_3}
\left[g_{yu}\right]^{\gamma_2}_{\alpha_2\beta_1}
\left[g_u\right]^{\gamma_3}_{\beta_2}+\cdots
$$
[B_y²u³]_α₁α₂β₁β₂β₃ = … + [f_z³]_γ₁γ₂γ₃ [z_yu]^γ₁_α₁β₃ [z_yu]^γ₂_α₂β₁ [zᵤ]^γ₃_β₂ + …
If the tensors are unfolded, we can evaluate this expression as
$$g_{y^3}\cdot\left(g_{yu}\otimes g_{yu}\otimes g_{u}\right)\cdot P,$$
where $P$ is a suitable permutation of columns of the expressions,
which permutes them so that the index
$(\alpha_1,\beta_3,\alpha_2,\beta_1,\beta_2)$ would go to
$(\alpha_1,\alpha_2,\beta_1,\beta_2,\beta_3)$.
The permutation $P$ can be very ineffective (copying great amount of
f_z³·(z_yu ⊗ z_yu ⊗ zᵤ)·P, where P is a suitable permutation of columns of
the expressions, which permutes them so that the index
(α₁,β₃,α₂,β₁,β₂) would go to (α₁,α₂,β₁,β₂,β₃).
The permutation P can be very ineffective (copying great amount of
small chunks of data) if the equivalence class ordering is chosen
badly. However, we do not provide any heuristic minimizing a total
time spent in all permutations. We choose an ordering which orders the
@ -46,10 +42,10 @@
#include <string>
/* Here is the abstraction for an equivalence class. We implement it as
|vector<int>|. We have a constructor for empty class, copy
vector<int>. We have a constructor for empty class, copy
constructor. What is important here is the ordering operator
|operator<| and methods for addition of an integer, and addition of
another sequence. Also we provide method |has| which returns true if a
operator<() and methods for addition of an integer, and addition of
another sequence. Also we provide method has() which returns true if a
given integer is contained. */
class OrdSequence
@ -81,10 +77,10 @@ private:
};
/* Here is the abstraction for the equivalence. It is a list of
equivalence classes. Also we remember |n|, which is a size of
underlying set $\{0, 1, \ldots, n-1\}$.
equivalence classes. Also we remember n, which is a size of
underlying set {0, 1, …, n-1}.
Method |trace| ``prints'' the equivalence into the integer sequence. */
Method trace() “prints” the equivalence into the integer sequence. */
class Permutation;
class Equivalence
@ -96,11 +92,11 @@ public:
using const_seqit = std::list<OrdSequence>::const_iterator;
using seqit = std::list<OrdSequence>::iterator;
// Constructs $\{\{0\},\{1\},\ldots,\{n-1\}\}$
// Constructs { {0}, {1}, …, {n-1} }
explicit Equivalence(int num);
// Constructs $\{\{0,1,\ldots,n-1\}\}$
// Constructs { {0,1,…,n-1 } }
Equivalence(int num, const std::string &dummy);
// Copy constructor plus gluing |i1| and |i2| in one class
// Copy constructor plus gluing i1 and i2 in one class
Equivalence(const Equivalence &e, int i1, int i2);
bool operator==(const Equivalence &e) const;
@ -154,7 +150,7 @@ protected:
given number or we can find an equivalence class of a given index within
the ordering.
We have also an |insert| method which inserts a given class
We have also an insert() method which inserts a given class
according to the class ordering. */
const_seqit findHaving(int i) const;
seqit findHaving(int i);
@ -162,12 +158,12 @@ protected:
};
/* The |EquivalenceSet| is a list of equivalences. The unique
constructor constructs a set of all equivalences over $n$-element
/* The EquivalenceSet is a list of equivalences. The unique
constructor constructs a set of all equivalences over an n-elements
set. The equivalences are sorted in the list so that equivalences with
fewer number of classes are in the end.
The two methods |has| and |addParents| are useful in the constructor. */
The two methods has() and addParents() are useful in the constructor. */
class EquivalenceSet
{
@ -191,12 +187,12 @@ private:
void addParents(const Equivalence &e, std::list<Equivalence> &added);
};
/* The equivalence bundle class only encapsulates |EquivalenceSet|s
/* The equivalence bundle class only encapsulates EquivalenceSet·s
from 1 up to a given number. It is able to retrieve the equivalence set
over $n$-element set for a given $n$, and also it can generate some more
over n-element set for a given n, and also it can generate some more
sets on request.
It is fully responsible for storage needed for |EquivalenceSet|s. */
It is fully responsible for storage needed for EquivalenceSet·s. */
class EquivalenceBundle
{

2
dynare++/tl/cc/fine_container.cc
View File

@ -5,7 +5,7 @@
#include <cmath>
/* Here we construct the vector of new sizes of containers (before
|nc|) and copy all remaining sizes behind |nc|. */
nc) and copy all remaining sizes behind nc. */
SizeRefinement::SizeRefinement(const IntSequence &s, int nc, int max)
{

38
dynare++/tl/cc/fine_container.hh
View File

@ -12,14 +12,14 @@
of each stack in the refined container. The resulting object is stack
container, so everything works seamlessly.
We define here a class for refinement of sizes |SizeRefinement|, this
We define here a class for refinement of sizes SizeRefinement, this
is purely an auxiliary class allowing us to write a code more
concisely. The main class of this file is |FineContainer|, which
corresponds to refining. The two more classes |FoldedFineContainer|
and |UnfoldedFineContainer| are its specializations.
concisely. The main class of this file is FineContainer, which
corresponds to refining. The two more classes FoldedFineContainer
and UnfoldedFineContainer are its specializations.
NOTE: This code was implemented with a hope that it will help to cut
down memory allocations during the Faa Di Bruno formula
down memory allocations during the Faà Di Bruno formula
evaluation. However, it seems that this needs to be accompanied with a
similar thing for tensor multidimensional index. Thus, the abstraction
is not currently used, but it might be useful in future. */
@ -32,9 +32,9 @@
#include <vector>
#include <memory>
/* This class splits the first |nc| elements of the given sequence |s|
to a sequence not having items greater than given |max|. The remaining
elements (those behind |nc|) are left untouched. It also remembers the
/* This class splits the first nc elements of the given sequence s
to a sequence not having items greater than given max. The remaining
elements (those behind nc) are left untouched. It also remembers the
mapping, i.e. for a given index in a new sequence, it is able to
return a corresponding index in old sequence. */
@ -68,7 +68,7 @@ public:
};
/* This main class of this class refines a given stack container, and
inherits from the stack container. It also defines the |getType|
inherits from the stack container. It also defines the getType()
method, which returns a type for a given stack as the type of the
corresponding (old) stack of the former stack container. */
@ -82,19 +82,19 @@ protected:
std::vector<std::unique_ptr<_Ctype>> ref_conts;
const _Stype &stack_cont;
public:
/* Here we construct the |SizeRefinement| and allocate space for the
/* Here we construct the SizeRefinement and allocate space for the
refined containers. Then, the containers are created and put to
|conts| array. Note that the containers do not claim any further
conts array. Note that the containers do not claim any further
space, since all the tensors of the created containers are in-place
submatrices.
Here we use a dirty trick of converting |const| pointer to non-|const|
Here we use a dirty trick of converting const pointer to non-const
pointer and passing it to a subtensor container constructor. The
containers are stored in |ref_conts| and then in |conts| from
|StackContainer|. However, this is safe since neither |ref_conts| nor
|conts| are used in non-|const| contexts. For example,
|StackContainer| has only a |const| method to return a member of
|conts|. */
containers are stored in ref_conts and then in conts from
StackContainer. However, this is safe since neither ref_conts nor
conts are used in non-const contexts. For example,
StackContainer has only a const method to return a member of
conts. */
FineContainer(const _Stype &sc, int max)
: SizeRefinement(sc.getStackSizes(), sc.numConts(), max),
@ -130,7 +130,7 @@ public:
};
/* Here is |FineContainer| specialization for folded tensors. */
/* Here is FineContainer specialization for folded tensors. */
class FoldedFineContainer : public FineContainer<FGSTensor>, public FoldedStackContainer
{
public:
@ -140,7 +140,7 @@ public:
}
};
/* Here is |FineContainer| specialization for unfolded tensors. */
/* Here is FineContainer specialization for unfolded tensors. */
class UnfoldedFineContainer : public FineContainer<UGSTensor>, public UnfoldedStackContainer
{
public:

36
dynare++/tl/cc/fs_tensor.cc
View File

@ -8,13 +8,13 @@
#include "pascal_triangle.hh"
/* This constructs a fully symmetric tensor as given by the contraction:
$$\left[g_{y^n}\right]_{\alpha_1\ldots\alpha_n}=
\left[t_{y^{n+1}}\right]_{\alpha_1\ldots\alpha_n\beta}[x]^\beta$$
We go through all columns of output tensor $[g]$ and for each column
[g_yⁿ]_α₁…αₙ = [t_yⁿ⁺¹]_α₁…αₙβ [x]^β
We go through all columns of output tensor [g] and for each column
we cycle through all variables, insert a variable to the column
coordinates obtaining a column of tensor $[t]$. the column is multiplied
by an appropriate item of |x| and added to the column of $[g]$ tensor. */
coordinates obtaining a column of tensor [t]. The column is multiplied
by an appropriate item of x and added to the column of [g] tensor. */
FFSTensor::FFSTensor(const FFSTensor &t, const ConstVector &x)
: FTensor(indor::along_col, IntSequence(t.dimen()-1, t.nvar()),
@ -38,8 +38,10 @@ FFSTensor::FFSTensor(const FFSTensor &t, const ConstVector &x)
}
/* This returns number of indices for folded tensor with full
symmetry. Let $n$ be a number of variables |nvar| and $d$ the
dimension |dim|. Then the number of indices is $\pmatrix{n+d-1\cr d}$. */
symmetry. Let n be a number of variables and d the
⎛n+d-1⎞
dimension dim. Then the number of indices is ⎝ d ⎠.
*/
int
FFSTensor::calcMaxOffset(int nvar, int d)
@ -91,7 +93,7 @@ FFSTensor::unfold() const
}
/* Incrementing is easy. We have to increment by calling static method
|UTensor::increment| first. In this way, we have coordinates of
UTensor::increment() first. In this way, we have coordinates of
unfolded tensor. Then we have to skip to the closest folded index
which corresponds to monotonizeing the integer sequence. */
@ -105,7 +107,7 @@ FFSTensor::increment(IntSequence &v) const
v.monotone();
}
/* Decrement calls static |FTensor::decrement|. */
/* Decrement calls static FTensor::decrement(). */
void
FFSTensor::decrement(IntSequence &v) const
@ -143,7 +145,7 @@ FFSTensor::addSubTensor(const FGSTensor &t)
TL_RAISE_IF(nvar() != t.getDims().getNVS().sum(),
"Wrong nvs for FFSTensor::addSubTensor");
// set shift for |addSubTensor|
// set shift for addSubTensor()
/* Code shared with UFSTensor::addSubTensor() */
IntSequence shift_pre(t.getSym().num(), 0);
for (int i = 1; i < t.getSym().num(); i++)
@ -160,8 +162,8 @@ FFSTensor::addSubTensor(const FGSTensor &t)
}
}
// |UFSTensor| contraction constructor
/* This is a bit more straightforward than |@<|FFSTensor| contraction constructor@>|.
// UFSTensor contraction constructor
/* This is a bit more straightforward than FFSTensor contraction constructor.
We do not add column by column but we do it by submatrices due to
regularity of the unfolded tensor. */
@ -186,7 +188,7 @@ UFSTensor::UFSTensor(const UFSTensor &t, const ConstVector &x)
}
/* Here we convert folded full symmetry tensor to unfolded. We copy all
columns of folded tensor, and then call |unfoldData()|. */
columns of folded tensor, and then call unfoldData(). */
UFSTensor::UFSTensor(const FFSTensor &ft)
: UTensor(indor::along_col, IntSequence(ft.dimen(), ft.nvar()),
@ -207,8 +209,8 @@ UFSTensor::fold() const
return std::make_unique<FFSTensor>(*this);
}
// |UFSTensor| increment and decrement
/* Here we just call |UTensor| respective static methods. */
// UFSTensor increment and decrement
/* Here we just call UTensor respective static methods. */
void
UFSTensor::increment(IntSequence &v) const
{
@ -236,7 +238,7 @@ UFSTensor::getOffset(const IntSequence &v) const
return UTensor::getOffset(v, nv);
}
/* This is very similar to |@<|FFSTensor::addSubTensor| code@>|. The
/* This is very similar to FFSTensor::addSubTensor(). The
only difference is the addition. We go through all columns in the full
symmetry tensor and cancel the shift. If the coordinates after the
cancellation are positive, we find the column in the general symmetry
@ -250,7 +252,7 @@ UFSTensor::addSubTensor(const UGSTensor &t)
TL_RAISE_IF(nvar() != t.getDims().getNVS().sum(),
"Wrong nvs for UFSTensor::addSubTensor");
// set shift for |addSubTensor|
// set shift for addSubTensor()
/* Code shared with FFSTensor::addSubTensor() */
IntSequence shift_pre(t.getSym().num(), 0);
for (int i = 1; i < t.getSym().num(); i++)

25
dynare++/tl/cc/fs_tensor.hh
View File

@ -19,21 +19,21 @@ class FRSingleTensor;
class FSSparseTensor;
/* Folded tensor with full symmetry maintains only information about
number of symmetrical variables |nv|. Further, we implement what is
left from the super class |FTensor|.
number of symmetrical variables nv. Further, we implement what is
left from the super class FTensor.
We implement |getOffset| which should be used with care since
We implement getOffset() which should be used with care since
its complexity.
We implement a method adding a given general symmetry tensor to the
full symmetry tensor supposing the variables of the general symmetry
tensor are stacked giving only one variable of the full symmetry
tensor. For instance, if $x=[y^T, u^T]^T$, then we can add tensor
$\left[g_{y^2u}\right]$ to tensor $g_{x^3}$. This is done in method
|addSubTensor|. Consult |@<|FGSTensor| class declaration@>| to know
tensor. For instance, if x=[yᵀ, uᵀ]ᵀ, then we can add tensor
[g_y²u] to tensor [g_x³]. This is done in method
addSubTensor(). Consult FGSTensor class declaration to know
what is general symmetry tensor.
Note that the past-the-end index is of the form $(nv, \ldots, nv)$, because
Note that the past-the-end index is of the form (nv,…,nv), because
of the specific implementation of FFSTensor::increment().
*/
@ -53,10 +53,11 @@ public:
}
/* Constructs a tensor by one-dimensional contraction from the higher
dimensional tensor |t|. This is, it constructs a tensor
$$\left[g_{y^n}\right]_{\alpha_1\ldots\alpha_n}=
\left[t_{y^{n+1}}\right]_{\alpha_1\ldots\alpha_n\beta}[x]^\beta$$ See the
implementation for details. */
dimensional tensor t. This is, it constructs a tensor
[g_yⁿ]_α₁…αₙ = [t_yⁿ⁺¹]_α₁…αₙβ [x]^β
See the implementation for details. */
FFSTensor(const FFSTensor &t, const ConstVector &x);
/* Converts from sparse tensor (which is fully symmetric and folded by
@ -95,7 +96,7 @@ public:
};
/* Unfolded fully symmetric tensor is almost the same in structure as
|FFSTensor|, but the method |unfoldData|. It takes columns which also
FFSTensor, but the method unfoldData(). It takes columns which also
exist in folded version and copies them to all their symmetrical
locations. This is useful when constructing unfolded tensor from
folded one. */

4
dynare++/tl/cc/ps_tensor.hh
View File

@ -269,11 +269,11 @@ protected:
};
/* Here we define an abstraction of the permuted symmetry folded
tensor. It is needed in context of the Faa Di Bruno formula for folded
tensor. It is needed in context of the Faà Di Bruno formula for folded
stack container multiplied with container of dense folded tensors, or
multiplied by one full symmetry sparse tensor.
For example, if we perform the Faa Di Bruno for $F=f(z)$, where
For example, if we perform the Faà Di Bruno for $F=f(z)$, where
$z=[g(x,y,u,v), h(x,y,u), x, y]^T$, we get for one concrete
equivalence:
$$

4
dynare++/tl/cc/pyramid_prod.hh
View File

@ -2,8 +2,8 @@
// Multiplying tensor columns.
/* In here, we implement the Faa Di Bruno for folded
tensors. Recall, that one step of the Faa Di Bruno is a formula:
/* In here, we implement the Faà Di Bruno for folded
tensors. Recall, that one step of the Faà Di Bruno is a formula:
$$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_k}=
[h_{y^l}]_{\gamma_1\ldots\gamma_l}
\prod_{m=1}^l\left[g_{s^{\vert c_m\vert}}\right]^{\gamma_m}_{c_m(\alpha)}

2
dynare++/tl/cc/stack_container.cc
View File

@ -28,7 +28,7 @@ FoldedStackContainer::multAndAdd(const FSSparseTensor &t,
}
// |FoldedStackContainer::multAndAdd| dense code
/* Here we perform the Faa Di Bruno step for a given dimension |dim|, and for
/* Here we perform the Faà Di Bruno step for a given dimension |dim|, and for
the dense fully symmetric tensor which is scattered in the container
of general symmetric tensors. The implementation is pretty the same as
|@<|UnfoldedStackContainer::multAndAdd| dense code@>|. */

41
dynare++/tl/cc/stack_container.hh
View File

@ -3,28 +3,35 @@
// Stack of containers.
/* Here we develop abstractions for stacked containers of tensors. For
instance, in perturbation methods for DSGE we need function
$$z(y,u,u',\sigma)=\left[\matrix{G(y,u,u',\sigma)\cr g(y,u,\sigma)\cr y\cr u}\right]$$
and we need to calculate one step of Faa Di Bruno formula
$$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_l}=\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
\sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s^k(c_m)}\right]^{\beta_m}_{c_m(\alpha)}$$
where we have containers for derivatives of $G$ and $g$.
instance, in perturbation methods for DSGE models, we need the function:
⎡ G(y*,u,u′,σ) ⎤
z(y*,u,u′,σ) = ⎢ g(y*,u,σ) ⎥
⎢ y* ⎥
⎣ u ⎦
and we need to calculate one step of Faà Di Bruno formula:
ₗ
[B_sᵏ]_α₁,…,αₗ = [f_zˡ]_β₁,…,βₗ ∑ ∏ [z_(s^|cₘ|)]_cₘ(α)^βₘ
c∈ℳₗ,ₖ ᵐ⁼¹
where we have containers for derivatives of G and g.
The main purpose of this file is to define abstractions for stack of
containers and possibly raw variables, and code |multAndAdd| method
calculating (one step of) the Faa Di Bruno formula for folded and
unfolded tensors. Note also, that tensors $\left[f_{z^l}\right]$ are
sparse.
containers and possibly raw variables, and code multAndAdd() method
calculating (one step of) the Faà Di Bruno formula for folded and
unfolded tensors. Note also, that tensors [f_zˡ] are sparse.
The abstractions are built as follows. At the top, there is an
interface describing stack of columns. It contains pure virtual
methods needed for manipulating the container stack. For technical
reasons it is a template. Both versions (folded, and unfolded) provide
all interface necessary for implementation of |multAndAdd|. The second
all interface necessary for implementation of multAndAdd(). The second
way of inheritance is first general implementation of the interface
|StackContainer|, and then specific (|ZContainer| for our specific
$z$). The only method which is virtual also after |StackContainer| is
|getType|, which is implemented in the specialization and determines
StackContainer, and then specific (ZContainer for our specific z).
The only method which is virtual also after StackContainer is
getType(), which is implemented in the specialization and determines
behaviour of the stack. The complete classes are obtained by
inheriting from the both branches, as it is drawn below:
@ -53,7 +60,7 @@
{|UnfoldedStackContainer|}{|ZContainer<UGSTensor>|}{|UnfoldedZContainer|}
}
We have also two supporting classes |StackProduct| and |KronProdStack|
We have also two supporting classes StackProduct and KronProdStack
and a number of worker classes used as threads. */
#ifndef STACK_CONTAINER_H
@ -85,8 +92,8 @@
Method |createPackedColumn| returns a vector of stack derivatives with
respect to the given symmetry and of the given column, where all zeros
from zero types, or unit matrices are deleted. See {\tt
kron\_prod2.hweb} for explanation. */
from zero types, or unit matrices are deleted. See kron_prod.hh for an
explanation. */
template <class _Ttype>
class StackContainerInterface

2
dynare++/tl/cc/t_container.cc
View File

@ -55,7 +55,7 @@ FGSContainer::FGSContainer(const UGSContainer &c)
}
// |FGSContainer::multAndAdd| folded code
/* Here we perform one step of the Faa Di Bruno operation. We call the
/* Here we perform one step of the Faà Di Bruno operation. We call the
|multAndAdd| for unfolded tensor. */
void
FGSContainer::multAndAdd(const FGSTensor &t, FGSTensor &out) const

6
dynare++/tl/cc/t_container.hh
View File

@ -3,7 +3,7 @@
// Tensor containers.
/* One of primary purposes of the tensor library is to perform one step
of the Faa Di Bruno formula:
of the Faà Di Bruno formula:
$$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_k}=
[h_{y^l}]_{\gamma_1\ldots\gamma_l}\sum_{c\in M_{l,k}}
\prod_{m=1}^l\left[g_{s^{\vert c_m\vert}}\right]^{\gamma_m}_{c_m(\alpha)}
@ -30,14 +30,14 @@
most one tensor.
The class has two purposes: The first is to provide storage (insert
and retrieve). The second is to perform the above step of Faa Di Bruno. This is
and retrieve). The second is to perform the above step of Faà Di Bruno. This is
going through all equivalences with $l$ classes, perform the tensor
product and add to the result.
We define a template class |TensorContainer|. From different
instantiations of the template class we will inherit to create concrete
classes, for example container of unfolded general symmetric
tensors. The one step of the Faa Di Bruno (we call it |multAndAdd|) is
tensors. The one step of the Faà Di Bruno (we call it |multAndAdd|) is
implemented in the concrete subclasses, because the implementation
depends on storage. Note even, that |multAndAdd| has not a template
common declaration. This is because sparse tensor $h$ is multiplied by

11
dynare++/utils/cc/pascal_triangle.cc
View File

@ -11,10 +11,13 @@ PascalRow::setFromPrevious(const PascalRow &prev)
prolong(prev);
}
/** This prolongs the PascalRow. If it is empty, we set the first item
* to k+1, which is noverk(k+1,k) which is the second item in the real
* pascal row, which starts from noverk(k,k)=1. Then we calculate
* other items from the provided row which must be the one with k-1.*/
/* This prolongs the PascalRow. If it is empty, we set the first item
⎛k+1⎞
to k+1, which is ⎝ k ⎠, which is the second item in the real
⎛k⎞
pascal row, which starts from ⎝k⎠=1. Then we calculate
other items from the provided row which must be the one with k-1.
*/
void
PascalRow::prolong(const PascalRow &prev)
{

3
dynare++/utils/cc/pascal_triangle.hh
View File

@ -24,6 +24,9 @@ public:
namespace PascalTriangle
{
void ensure(int n, int k);
/* ⎛n⎞
Computes ⎝k⎠, hence the function name (“n over k”).
*/
int noverk(int n, int k);
void print();
};