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

Dynare++: give more explicit names to some methods on quasi-triangular matrices

time-shift
Sébastien Villemot 2019-03-29 15:27:34 +01:00
parent 39896d8d0a
commit c41216819a
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 56 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -16,7 +16,7 @@ class BlockDiagonal : public QuasiTriangular
public:
BlockDiagonal(ConstVector d, int d_size);
BlockDiagonal(const BlockDiagonal &b) = default;
BlockDiagonal(const QuasiTriangular &t);
explicit BlockDiagonal(const QuasiTriangular &t);
BlockDiagonal &operator=(const QuasiTriangular &t)
{
GeneralMatrix::operator=(t);

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

@ -79,17 +79,13 @@ Diagonal::Diagonal(double *data, int d_size)
if ((jbar < d_size-1) && !isZero(data[ill]))
{
// it is not last column and we have nonzero below diagonal
DiagonalBlock b(jbar, false, &data[id], &data[idd],
&data[iur], &data[ill]);
blocks.push_back(b);
blocks.emplace_back(jbar, false, &data[id], &data[idd],
&data[iur], &data[ill]);
jbar++;
}
else
{
// it is last column or we have zero below diagonal
DiagonalBlock b(jbar, true, &data[id], &data[id], nullptr, nullptr);
blocks.push_back(b);
}
// it is last column or we have zero below diagonal
blocks.emplace_back(jbar, true, &data[id], &data[id], nullptr, nullptr);
jbar++;
j++;
}
@ -100,21 +96,20 @@ Diagonal::Diagonal(double *data, const Diagonal &d)
num_all = d.num_all;
num_real = d.num_real;
int d_size = d.getSize();
for (const auto & dit : d)
for (const auto &block : d)
{
double *beta1 = nullptr;
double *beta2 = nullptr;
int id = dit.getIndex()*(d_size+1);
int id = block.getIndex()*(d_size+1);
int idd = id;
if (!dit.isReal())
if (!block.isReal())
{
beta1 = &data[id+d_size];
beta2 = &data[id+1];
idd = id + d_size + 1;
}
DiagonalBlock b(dit.getIndex(), dit.isReal(),
&data[id], &data[idd], beta1, beta2);
blocks.push_back(b);
blocks.emplace_back(block.getIndex(), block.isReal(),
&data[id], &data[idd], beta1, beta2);
}
}
@ -231,8 +226,7 @@ Diagonal::checkConsistency(diag_iter it)
it->real = true;
it->beta1 = nullptr;
it->beta2 = nullptr;
DiagonalBlock b(jbar+1, d2);
blocks.insert(++it, b);
blocks.emplace(++it, jbar+1, d2);
num_real += 2;
num_all++;
}
@ -397,11 +391,11 @@ QuasiTriangular::QuasiTriangular(double r, const QuasiTriangular &t)
}
QuasiTriangular::QuasiTriangular(double r, const QuasiTriangular &t,
double rr, const QuasiTriangular &tt)
double r2, const QuasiTriangular &t2)
: SqSylvMatrix(t.numRows()), diagonal(getData().base(), t.diagonal)
{
setMatrix(r, t);
addMatrix(rr, tt);
addMatrix(r2, t2);
}
QuasiTriangular::QuasiTriangular(const QuasiTriangular &t)
@ -414,7 +408,7 @@ QuasiTriangular::QuasiTriangular(const ConstVector &d, int d_size)
{
}
QuasiTriangular::QuasiTriangular(int p, const QuasiTriangular &t)
QuasiTriangular::QuasiTriangular(const std::string &dummy, const QuasiTriangular &t)
: SqSylvMatrix(t.numRows()), diagonal(getData().base(), t.diagonal)
{
Vector aux(t.getData());

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

@ -332,12 +332,15 @@ protected:
Diagonal diagonal;
public:
QuasiTriangular(const ConstVector &d, int d_size);
// Initializes with r·t
QuasiTriangular(double r, const QuasiTriangular &t);
// Initializes with r·t+r₂·t₂
QuasiTriangular(double r, const QuasiTriangular &t,
double rr, const QuasiTriangular &tt);
QuasiTriangular(int p, const QuasiTriangular &t);
QuasiTriangular(const SchurDecomp &decomp);
QuasiTriangular(const SchurDecompZero &decomp);
double r2, const QuasiTriangular &t2);
// Initializes with t²
QuasiTriangular(const std::string &dummy, const QuasiTriangular &t);
explicit QuasiTriangular(const SchurDecomp &decomp);
explicit QuasiTriangular(const SchurDecompZero &decomp);
QuasiTriangular(const QuasiTriangular &t);
~QuasiTriangular() override = default;
@ -414,20 +417,23 @@ public:
{
return std::make_unique<QuasiTriangular>(*this);
}
// Returns this²
virtual std::unique_ptr<QuasiTriangular>
clone(int p, const QuasiTriangular &t) const
square() const
{
return std::make_unique<QuasiTriangular>(p, t);
return std::make_unique<QuasiTriangular>("square", *this);
}
// Returns r·this
virtual std::unique_ptr<QuasiTriangular>
clone(double r) const
scale(double r) const
{
return std::make_unique<QuasiTriangular>(r, *this);
}
// Returns r·this + r₂·t₂
virtual std::unique_ptr<QuasiTriangular>
clone(double r, double rr, const QuasiTriangular &tt) const
linearlyCombine(double r, double r2, const QuasiTriangular &t2) const
{
return std::make_unique<QuasiTriangular>(r, *this, rr, tt);
return std::make_unique<QuasiTriangular>(r, *this, r2, t2);
}
protected:
// this = r·t

12
dynare++/sylv/cc/QuasiTriangularZero.cc
View File

@ -30,18 +30,18 @@ QuasiTriangularZero::QuasiTriangularZero(double r,
QuasiTriangularZero::QuasiTriangularZero(double r,
const QuasiTriangularZero &t,
double rr,
const QuasiTriangularZero &tt)
: QuasiTriangular(r, t, rr, tt),
double r2,
const QuasiTriangularZero &t2)
: QuasiTriangular(r, t, r2, t2),
nz(t.nz),
ru(t.ru)
{
ru.mult(r);
ru.add(rr, tt.ru);
ru.add(r2, t2.ru);
}
QuasiTriangularZero::QuasiTriangularZero(int p, const QuasiTriangularZero &t)
: QuasiTriangular(p, t),
QuasiTriangularZero::QuasiTriangularZero(const std::string &dummy, const QuasiTriangularZero &t)
: QuasiTriangular(dummy, t),
nz(t.nz),
ru(t.ru)
{

23
dynare++/sylv/cc/QuasiTriangularZero.hh
View File

@ -16,12 +16,15 @@ class QuasiTriangularZero : public QuasiTriangular
GeneralMatrix ru; // data in right upper part (nz,d_size)
public:
QuasiTriangularZero(int num_zeros, const ConstVector &d, int d_size);
// Initializes with r·t
QuasiTriangularZero(double r, const QuasiTriangularZero &t);
// Initializes with r·t+r₂·t₂
QuasiTriangularZero(double r, const QuasiTriangularZero &t,
double rr, const QuasiTriangularZero &tt);
QuasiTriangularZero(int p, const QuasiTriangularZero &t);
QuasiTriangularZero(const QuasiTriangular &t);
QuasiTriangularZero(const SchurDecompZero &decomp);
double r2, const QuasiTriangularZero &t2);
// Initializes with t²
QuasiTriangularZero(const std::string &dummy, const QuasiTriangularZero &t);
explicit QuasiTriangularZero(const QuasiTriangular &t);
explicit QuasiTriangularZero(const SchurDecompZero &decomp);
~QuasiTriangularZero() override = default;
void solvePre(Vector &x, double &eig_min) override;
void solvePreTrans(Vector &x, double &eig_min) override;
@ -32,26 +35,26 @@ public:
void multKron(KronVector &x) const override;
void multKronTrans(KronVector &x) const override;
void multLeftOther(GeneralMatrix &a) const override;
/* clone */
std::unique_ptr<QuasiTriangular>
clone() const override
{
return std::make_unique<QuasiTriangularZero>(*this);
}
std::unique_ptr<QuasiTriangular>
clone(int p, const QuasiTriangular &t) const override
square() const override
{
return std::make_unique<QuasiTriangularZero>(p, dynamic_cast<const QuasiTriangularZero &>(t));
return std::make_unique<QuasiTriangularZero>("square", *this);
}
std::unique_ptr<QuasiTriangular>
clone(double r) const override
scale(double r) const override
{
return std::make_unique<QuasiTriangularZero>(r, *this);
}
std::unique_ptr<QuasiTriangular>
clone(double r, double rr, const QuasiTriangular &tt) const override
linearlyCombine(double r, double r2, const QuasiTriangular &t2) const override
{
return std::make_unique<QuasiTriangularZero>(r, *this, rr, dynamic_cast<const QuasiTriangularZero &>(tt));
return std::make_unique<QuasiTriangularZero>(r, *this, r2, dynamic_cast<const QuasiTriangularZero &>(t2));
}
void print() const override;
};

16
dynare++/sylv/cc/TriangularSylvester.cc
View File

@ -13,24 +13,24 @@
TriangularSylvester::TriangularSylvester(const QuasiTriangular &k,
const QuasiTriangular &f)
: SylvesterSolver(k, f),
matrixKK(matrixK->clone(2, *matrixK)),
matrixFF(std::make_unique<QuasiTriangular>(2, *matrixF))
matrixKK{matrixK->square()},
matrixFF{matrixF->square()}
{
}
TriangularSylvester::TriangularSylvester(const SchurDecompZero &kdecomp,
const SchurDecomp &fdecomp)
: SylvesterSolver(kdecomp, fdecomp),
matrixKK(matrixK->clone(2, *matrixK)),
matrixFF(std::make_unique<QuasiTriangular>(2, *matrixF))
matrixKK{matrixK->square()},
matrixFF{matrixF->square()}
{
}
TriangularSylvester::TriangularSylvester(const SchurDecompZero &kdecomp,
const SimilarityDecomp &fdecomp)
: SylvesterSolver(kdecomp, fdecomp),
matrixKK(matrixK->clone(2, *matrixK)),
matrixFF(std::make_unique<BlockDiagonal>(2, *matrixF))
matrixKK{matrixK->square()},
matrixFF{matrixF->square()}
{
}
@ -56,7 +56,7 @@ TriangularSylvester::solvi(double r, KronVector &d, double &eig_min) const
{
if (d.getDepth() == 0)
{
auto t = matrixK->clone(r);
auto t = matrixK->scale(r);
t->solvePre(d, eig_min);
}
else
@ -98,7 +98,7 @@ TriangularSylvester::solviip(double alpha, double betas,
if (d.getDepth() == 0)
{
double aspbs = alpha*alpha+betas;
auto t = matrixK->clone(2*alpha, aspbs, *matrixKK);
auto t = matrixK->linearlyCombine(2*alpha, aspbs, *matrixKK);
t->solvePre(d, eig_min);
}
else