diff --git a/dynare++/sylv/cc/BlockDiagonal.cc b/dynare++/sylv/cc/BlockDiagonal.cc
index 83d173c24a8e7ad12e25905a52621b7b88ba34dc..88db350916c5a9cc2f7ca0b71a2ed72cb5899a91 100644
--- a/dynare++/sylv/cc/BlockDiagonal.cc
+++ b/dynare++/sylv/cc/BlockDiagonal.cc
@@ -4,7 +4,7 @@
#include "BlockDiagonal.hh"
-#include <cstdio>
+#include <iostream>
#include <cstring>
BlockDiagonal::BlockDiagonal(const double *d, int d_size)
@@ -315,7 +315,7 @@ BlockDiagonal::multKronTrans(KronVector &x) const
void
BlockDiagonal::printInfo() const
{
- printf("Block sizes:");
+ std::cout << "Block sizes:";
int num_blocks = 0;
const_diag_iter start = diag_begin();
const_diag_iter end = findBlockStart(start);
@@ -325,14 +325,14 @@ BlockDiagonal::printInfo() const
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
- printf(" %d", ei-si);
+ std::cout << ' ' << ei-si;
num_blocks++;
start = end;
end = findBlockStart(start);
}
- printf("\nNum blocks: %d\n", num_blocks);
- printf("There are %d zeros out of %d\n",
- getNumZeros(), getNumOffdiagonal());
+ std::cout << std::endl
+ << "Num blocks: " << num_blocks << std::endl
+ << "There are " << getNumZeros() << " zeros out of " << getNumOffdiagonal() << std::endl;
}
int
diff --git a/dynare++/sylv/cc/BlockDiagonal.hh b/dynare++/sylv/cc/BlockDiagonal.hh
index 3b108bbd328c7a1fbd6ebc9a5503c159d1c7b304..3b5c18546ec37efd7dd4f695383032509e20269d 100644
--- a/dynare++/sylv/cc/BlockDiagonal.hh
+++ b/dynare++/sylv/cc/BlockDiagonal.hh
@@ -18,12 +18,12 @@ public:
BlockDiagonal(int p, const BlockDiagonal &b);
BlockDiagonal(const BlockDiagonal &b);
BlockDiagonal(const QuasiTriangular &t);
- const BlockDiagonal &
+ BlockDiagonal &
operator=(const QuasiTriangular &t)
{
GeneralMatrix::operator=(t); return *this;
}
- const BlockDiagonal &operator=(const BlockDiagonal &b);
+ BlockDiagonal &operator=(const BlockDiagonal &b);
~BlockDiagonal() override
{
delete [] row_len; delete [] col_len;
diff --git a/dynare++/sylv/cc/GeneralMatrix.cc b/dynare++/sylv/cc/GeneralMatrix.cc
index 2c40db05384f8322e2286e562984e2705829abee..912a4657f64d3058c19ba40b77dfe73200691788 100644
--- a/dynare++/sylv/cc/GeneralMatrix.cc
+++ b/dynare++/sylv/cc/GeneralMatrix.cc
@@ -8,11 +8,13 @@
#include <dynblas.h>
#include <dynlapack.h>
-#include <cstdio>
+#include <iostream>
+#include <iomanip>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <limits>
+#include <vector>
int GeneralMatrix::md_length = 32;
@@ -80,9 +82,6 @@ GeneralMatrix::GeneralMatrix(const GeneralMatrix &a, const char *dum1,
gemm("T", a, "T", b, 1.0, 0.0);
}
-GeneralMatrix::~GeneralMatrix()
-= default;
-
void
GeneralMatrix::place(const ConstGeneralMatrix &m, int i, int j)
{
@@ -178,11 +177,9 @@ GeneralMatrix::zeros()
if (ld == rows)
data.zeros();
else
- {
- for (int i = 0; i < rows; i++)
- for (int j = 0; j < cols; j++)
- get(i, j) = 0;
- }
+ for (int i = 0; i < rows; i++)
+ for (int j = 0; j < cols; j++)
+ get(i, j) = 0;
}
void
@@ -219,11 +216,9 @@ GeneralMatrix::mult(double a)
if (ld == rows)
data.mult(a);
else
- {
- for (int i = 0; i < rows; i++)
- for (int j = 0; j < cols; j++)
- get(i, j) *= a;
- }
+ for (int i = 0; i < rows; i++)
+ for (int j = 0; j < cols; j++)
+ get(i, j) *= a;
}
// here we must be careful for ld
@@ -236,11 +231,9 @@ GeneralMatrix::add(double a, const ConstGeneralMatrix &m)
if (ld == rows && m.ld == m.rows)
data.add(a, m.data);
else
- {
- for (int i = 0; i < rows; i++)
- for (int j = 0; j < cols; j++)
- get(i, j) += a*m.get(i, j);
- }
+ for (int i = 0; i < rows; i++)
+ for (int j = 0; j < cols; j++)
+ get(i, j) += a*m.get(i, j);
}
void
@@ -396,9 +389,7 @@ void
ConstGeneralMatrix::multVec(double a, Vector &x, double b, const ConstVector &d) const
{
if (x.length() != rows || cols != d.length())
- {
- throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
if (rows > 0)
{
blas_int mm = rows;
@@ -419,9 +410,7 @@ ConstGeneralMatrix::multVecTrans(double a, Vector &x, double b,
const ConstVector &d) const
{
if (x.length() != cols || rows != d.length())
- {
- throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
if (rows > 0)
{
blas_int mm = rows;
@@ -441,24 +430,19 @@ void
ConstGeneralMatrix::multInvLeft(const char *trans, int mrows, int mcols, int mld, double *d) const
{
if (rows != cols)
- {
- throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
- }
+ throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
if (cols != mrows)
- {
- throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
if (rows > 0)
{
GeneralMatrix inv(*this);
- auto *ipiv = new lapack_int[rows];
+ std::vector<lapack_int> ipiv(rows);
lapack_int info;
lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld;
- dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv, &info);
- dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv, d,
+ dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv.data(), &info);
+ dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv.data(), d,
&mld2, &info);
- delete [] ipiv;
}
}
@@ -481,9 +465,7 @@ void
ConstGeneralMatrix::multInvLeft(Vector &d) const
{
if (d.skip() != 1)
- {
- throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
- }
+ throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
multInvLeft("N", d.length(), 1, d.length(), d.base());
}
@@ -493,9 +475,7 @@ void
ConstGeneralMatrix::multInvLeftTrans(Vector &d) const
{
if (d.skip() != 1)
- {
- throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
- }
+ throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
multInvLeft("T", d.length(), 1, d.length(), d.base());
}
@@ -523,16 +503,17 @@ ConstGeneralMatrix::isZero() const
void
ConstGeneralMatrix::print() const
{
- printf("rows=%d, cols=%d\n", rows, cols);
+ auto ff = std::cout.flags();
+ std::cout << "rows=" << rows << ", cols=" << cols << std::endl;
for (int i = 0; i < rows; i++)
{
- printf("row %d:\n", i);
+ std::cout << "row " << i << ':' << std::endl
+ << std::setprecision(3);
for (int j = 0; j < cols; j++)
- {
- printf("%6.3g ", get(i, j));
- }
- printf("\n");
+ std::cout << std::setw(6) << get(i, j) << ' ';
+ std::cout << std::endl;
}
+ std::cout.flags(ff);
}
void
@@ -563,16 +544,15 @@ SVDDecomp::construct(const GeneralMatrix &A)
lapack_int lwork = -1;
lapack_int info;
- auto *iwork = new lapack_int[8*minmn];
+ std::vector<lapack_int> iwork(8*minmn);
// query for optimal lwork
dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &tmpwork,
- &lwork, iwork, &info);
+ &lwork, iwork.data(), &info);
lwork = (lapack_int) tmpwork;
Vector work(lwork);
// do the decomposition
dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work.base(),
- &lwork, iwork, &info);
- delete [] iwork;
+ &lwork, iwork.data(), &info);
if (info < 0)
throw SYLV_MES_EXCEPTION("Internal error in SVDDecomp constructor");
if (info == 0)
diff --git a/dynare++/sylv/cc/GeneralMatrix.hh b/dynare++/sylv/cc/GeneralMatrix.hh
index e7ff5611db2d0d3cbe62827ad598ec8ee2f00d2d..968061e4dea5e42a614c05413830753d9067180e 100644
--- a/dynare++/sylv/cc/GeneralMatrix.hh
+++ b/dynare++/sylv/cc/GeneralMatrix.hh
@@ -145,7 +145,7 @@ public:
GeneralMatrix(const GeneralMatrix &a, const char *dum1,
const GeneralMatrix &b, const char *dum2);
- virtual ~GeneralMatrix();
+ virtual ~GeneralMatrix() = default;
GeneralMatrix &operator=(const GeneralMatrix &m) = default;
const double &
diff --git a/dynare++/sylv/cc/KronVector.cc b/dynare++/sylv/cc/KronVector.cc
index 494a8f13b13948d8f6956835eded9084f7c0e2a7..db27fa13cd505a8dd0bd9874c722521c841dee4c 100644
--- a/dynare++/sylv/cc/KronVector.cc
+++ b/dynare++/sylv/cc/KronVector.cc
@@ -40,7 +40,7 @@ KronVector::KronVector(const ConstKronVector &v)
Vector::operator=(v);
}
-const KronVector &
+KronVector &
KronVector::operator=(const ConstKronVector &v)
{
Vector::operator=(v);
@@ -50,7 +50,7 @@ KronVector::operator=(const ConstKronVector &v)
return *this;
}
-const KronVector &
+KronVector &
KronVector::operator=(const Vector &v)
{
if (length() != v.length())
diff --git a/dynare++/sylv/cc/KronVector.hh b/dynare++/sylv/cc/KronVector.hh
index 8bddd29f79e7b06dece37de623232ed98a10ca5f..f58d0f67df75cbe6a785e4b78f6945cf5fb8b66a 100644
--- a/dynare++/sylv/cc/KronVector.hh
+++ b/dynare++/sylv/cc/KronVector.hh
@@ -24,8 +24,8 @@ public:
KronVector(KronVector &, int i); // picks i-th subvector
KronVector(const ConstKronVector &v); // new instance and copy
KronVector &operator=(const KronVector &v) = default;
- const KronVector &operator=(const ConstKronVector &v);
- const KronVector &operator=(const Vector &v);
+ KronVector &operator=(const ConstKronVector &v);
+ KronVector &operator=(const Vector &v);
int
getM() const
{
diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc
index 348403a72f183fbef7dd7971435bc28b1d86382c..618c81db5db537f1f1dbfe10853102029c26754b 100644
--- a/dynare++/sylv/cc/QuasiTriangular.cc
+++ b/dynare++/sylv/cc/QuasiTriangular.cc
@@ -8,8 +8,9 @@
#include <dynblas.h>
-#include <cstdio>
#include <cmath>
+#include <iostream>
+#include <sstream>
using namespace std;
@@ -137,9 +138,7 @@ Diagonal::getNumComplex(const double *data, int d_size)
{
num_complex++;
if (in < d_size - 2 && !isZero(data[in + d_size +1]))
- {
- throw SYLV_MES_EXCEPTION("Matrix is not quasi-triangular");
- }
+ throw SYLV_MES_EXCEPTION("Matrix is not quasi-triangular");
}
}
return num_complex;
@@ -175,19 +174,17 @@ Diagonal::getEigenValues(Vector &eig) const
int d_size = getSize();
if (eig.length() != 2*d_size)
{
- char mes[500];
- sprintf(mes, "Wrong length of vector for eigenvalues len=%d, should be=%d.\n",
- eig.length(), 2*d_size);
- throw SYLV_MES_EXCEPTION(mes);
+ ostringstream mes;
+ mes << "Wrong length of vector for eigenvalues len=" << eig.length()
+ << ", should be=" << 2*d_size << '.' << std::endl;
+ throw SYLV_MES_EXCEPTION(mes.str());
}
for (const auto & b : *this)
{
int ind = b.getIndex();
eig[2*ind] = *(b.getAlpha());
if (b.isReal())
- {
- eig[2*ind+1] = 0.0;
- }
+ eig[2*ind+1] = 0.0;
else
{
double beta = sqrt(b.getSBeta());
@@ -304,23 +301,20 @@ Diagonal::findNextLargerBlock(diag_iter start, diag_iter end, double a)
void
Diagonal::print() const
{
- printf("Num real: %d, num complex: %d\n", getNumReal(), getNumComplex());
+ auto ff = std::cout.flags();
+ std::cout << "Num real: " << getNumReal() << ", num complex: " << getNumComplex() << std::endl
+ << std::fixed;
for (const auto & it : *this)
- {
- if (it.isReal())
- {
- printf("real: jbar=%d, alpha=%f\n", it.getIndex(), *(it.getAlpha()));
- }
- else
- {
- printf("complex: jbar=%d, alpha=%f, beta1=%f, beta2=%f\n",
- it.getIndex(), *(it.getAlpha()), it.getBeta1(), it.getBeta2());
- }
- }
+ if (it.isReal())
+ std::cout << "real: jbar=" << it.getIndex() << ", alpha=" << *(it.getAlpha()) << std::endl;
+ else
+ std::cout << "complex: jbar=" << it.getIndex()
+ << ", alpha=" << *(it.getAlpha())
+ << ", beta1=" << it.getBeta1()
+ << ", beta2=" << it.getBeta2() << std::endl;
+ std::cout.flags(ff);
}
-double Diagonal::EPS = 1.0e-300;
-
bool
Diagonal::isZero(double p)
{
diff --git a/dynare++/sylv/cc/QuasiTriangular.hh b/dynare++/sylv/cc/QuasiTriangular.hh
index d4f39847a3eea3e116257a32a17ec352553bc795..773508020064583405dce9cf0856ec159fbb86aa 100644
--- a/dynare++/sylv/cc/QuasiTriangular.hh
+++ b/dynare++/sylv/cc/QuasiTriangular.hh
@@ -100,7 +100,7 @@ public:
{
copy(b);
}
- const DiagonalBlock &
+ DiagonalBlock &
operator=(const DiagonalBlock &b)
{
copy(b); return *this;
@@ -179,7 +179,7 @@ public:
{
return x.it != it;
}
- const _Self &
+ _Self &
operator=(const _Self &x)
{
it = x.it; return *this;
@@ -209,7 +209,7 @@ public:
{
copy(d);
}
- const Diagonal &
+ Diagonal &
operator=(const Diagonal &d)
{
copy(d); return *this;
@@ -268,7 +268,7 @@ public:
/* redefine pointers as data start at p */
void changeBase(double *p);
private:
- static double EPS;
+ constexpr static double EPS = 1.0e-300;
static int getNumComplex(const double *data, int d_size);
static bool isZero(double p);
};
@@ -286,7 +286,7 @@ public:
ptr = base; d_size = ds; real = r;
}
virtual ~_matrix_iter() = default;
- const _Self &
+ _Self &
operator=(const _Self &it)
{
ptr = it.ptr; d_size = it.d_size; real = it.real; return *this;
diff --git a/dynare++/sylv/cc/QuasiTriangularZero.cc b/dynare++/sylv/cc/QuasiTriangularZero.cc
index 37bd9ff24e20ce0f722c73216d05592c1494ff85..665b6cd902f365d2b80c611de19d1a6ddf44b174 100644
--- a/dynare++/sylv/cc/QuasiTriangularZero.cc
+++ b/dynare++/sylv/cc/QuasiTriangularZero.cc
@@ -7,7 +7,7 @@
#include "SylvMatrix.hh"
#include "SylvException.hh"
-#include <cstdio>
+#include <iostream>
QuasiTriangularZero::QuasiTriangularZero(int num_zeros, const double *d,
int d_size)
@@ -62,9 +62,6 @@ QuasiTriangularZero::QuasiTriangularZero(const QuasiTriangular &t)
{
}
-QuasiTriangularZero::~QuasiTriangularZero()
-= default;
-
void
QuasiTriangularZero::solvePre(Vector &x, double &eig_min)
{
@@ -136,10 +133,10 @@ QuasiTriangularZero::multLeftOther(GeneralMatrix &a) const
void
QuasiTriangularZero::print() const
{
- printf("super=\n");
+ std::cout << "super=" << std::endl;
QuasiTriangular::print();
- printf("nz=%d\n", nz);
- printf("ru=\n");
+ std::cout << "nz=" << nz << std::endl
+ << "ru=" << std::endl;
ru.print();
}
diff --git a/dynare++/sylv/cc/QuasiTriangularZero.hh b/dynare++/sylv/cc/QuasiTriangularZero.hh
index f4bb4a9e5a535dd725c80bdf58ff4c5d98607151..154bc9ff3feba006b3cf65e2394b1f48deb35bd8 100644
--- a/dynare++/sylv/cc/QuasiTriangularZero.hh
+++ b/dynare++/sylv/cc/QuasiTriangularZero.hh
@@ -22,7 +22,7 @@ public:
QuasiTriangularZero(int p, const QuasiTriangularZero &t);
QuasiTriangularZero(const QuasiTriangular &t);
QuasiTriangularZero(const SchurDecompZero &decomp);
- ~QuasiTriangularZero() override;
+ ~QuasiTriangularZero() override = default;
void solvePre(Vector &x, double &eig_min) override;
void solvePreTrans(Vector &x, double &eig_min) override;
void multVec(Vector &x, const ConstVector &b) const override;
diff --git a/dynare++/sylv/cc/SchurDecomp.cc b/dynare++/sylv/cc/SchurDecomp.cc
index 1989dabdcb57ae112d58bea22de26b6d6c1eb63e..f4f85ba3d0bddbc5eca2648fba3c6bc63cf3d9ac 100644
--- a/dynare++/sylv/cc/SchurDecomp.cc
+++ b/dynare++/sylv/cc/SchurDecomp.cc
@@ -4,6 +4,8 @@
#include "SchurDecomp.hh"
+#include <vector>
+
#include <dynlapack.h>
SchurDecomp::SchurDecomp(const SqSylvMatrix &m)
@@ -12,17 +14,13 @@ SchurDecomp::SchurDecomp(const SqSylvMatrix &m)
lapack_int rows = m.numRows();
SqSylvMatrix auxt(m);
lapack_int sdim;
- auto *const wr = new double[rows];
- auto *const wi = new double[rows];
+ std::vector<double> wr(rows), wi(rows);
lapack_int lwork = 6*rows;
- auto *const work = new double[lwork];
+ std::vector<double> work(lwork);
lapack_int info;
dgees("V", "N", nullptr, &rows, auxt.base(), &rows, &sdim,
- wr, wi, q.base(), &rows,
- work, &lwork, nullptr, &info);
- delete [] work;
- delete [] wi;
- delete [] wr;
+ wr.data(), wi.data(), q.base(), &rows,
+ work.data(), &lwork, nullptr, &info);
t_storage = std::make_unique<QuasiTriangular>(auxt.base(), rows);
t = t_storage.get();
}
diff --git a/dynare++/sylv/cc/SchurDecompEig.cc b/dynare++/sylv/cc/SchurDecompEig.cc
index 394be8f37387dd6519a9673adac16d3e1b075845..3c5955d32625d01005abbfe41cd971414c6c7b6b 100644
--- a/dynare++/sylv/cc/SchurDecompEig.cc
+++ b/dynare++/sylv/cc/SchurDecompEig.cc
@@ -3,6 +3,8 @@
#include <dynlapack.h>
+#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
@@ -55,15 +57,12 @@ SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
lapack_int n = getDim();
lapack_int ifst = (*it).getIndex() + 1;
lapack_int ilst = (*itadd).getIndex() + 1;
- auto *work = new double[n];
+ std::vector<double> work(n);
lapack_int info;
- dtrexc("V", &n, getT().base(), &n, getQ().base(), &n, &ifst, &ilst, work,
+ dtrexc("V", &n, getT().base(), &n, getQ().base(), &n, &ifst, &ilst, work.data(),
&info);
- delete [] work;
if (info < 0)
- {
- throw SYLV_MES_EXCEPTION("Wrong argument to dtrexc.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong argument to dtrexc.");
if (info == 0)
{
diff --git a/dynare++/sylv/cc/SylvException.cc b/dynare++/sylv/cc/SylvException.cc
index 1065c283fc515c6a0d4b2187179cc8434fea5b62..599ea9bcc372a56c318aeb757ab59f2ddb5e74fc 100644
--- a/dynare++/sylv/cc/SylvException.cc
+++ b/dynare++/sylv/cc/SylvException.cc
@@ -4,58 +4,34 @@
#include "SylvException.hh"
-#include <cstring>
-#include <cstdio>
+#include <iostream>
+#include <utility>
-SylvException::SylvException(const char *f, int l)
+SylvException::SylvException(std::string f, int l)
+ : file{std::move(f)}, line{l}
{
- strcpy(file, f);
- line = l;
}
void
SylvException::printMessage() const
{
- char mes[1500];
- mes[0] = '\0';
- printMessage(mes, 1499);
- puts(mes);
+ std::cout << getMessage();
}
-int
-SylvException::printMessage(char *str, int maxlen) const
+std::string
+SylvException::getMessage() const
{
- int remain = maxlen;
- char aux[100];
- sprintf(aux, "From %s:%d\n", file, line);
- int newremain = remain - strlen(aux);
- if (newremain < 0)
- {
- aux[remain] = '\0';
- newremain = 0;
- }
- strcat(str, aux);
- return newremain;
+ return "From " + file + ':' + std::to_string(line) + '\n';
}
-SylvExceptionMessage::SylvExceptionMessage(const char *f, int i,
- const char *mes)
- : SylvException(f, i)
+SylvExceptionMessage::SylvExceptionMessage(std::string f, int i,
+ std::string mes)
+ : SylvException{std::move(f), i}, message{std::move(mes)}
{
- strcpy(message, mes);
}
-int
-SylvExceptionMessage::printMessage(char *str, int maxlen) const
+std::string
+SylvExceptionMessage::getMessage() const
{
- char aux[600];
- sprintf(aux, "At %s:%d:%s\n", file, line, message);
- int newremain = maxlen - strlen(aux);
- if (newremain < 0)
- {
- aux[maxlen] = '\0';
- newremain = 0;
- }
- strcat(str, aux);
- return newremain;
+ return "At " + file + ':' + std::to_string(line) + ':' + message + '\n';
}
diff --git a/dynare++/sylv/cc/SylvException.hh b/dynare++/sylv/cc/SylvException.hh
index 786a2a69a0ad9a947636ab625f09b62339c72beb..598c2e1e6c79fda7872fdff325fb28d01a6915f3 100644
--- a/dynare++/sylv/cc/SylvException.hh
+++ b/dynare++/sylv/cc/SylvException.hh
@@ -5,26 +5,28 @@
#ifndef SYLV_EXCEPTION_H
#define SYLV_EXCEPTION_H
+#include <string>
+
#include "SylvMemory.hh"
class SylvException : public MallocAllocator
{
protected:
- char file[50];
+ std::string file;
int line;
public:
- SylvException(const char *f, int l);
+ SylvException(std::string f, int l);
virtual ~SylvException() = default;
- virtual int printMessage(char *str, int maxlen) const;
void printMessage() const;
+ virtual std::string getMessage() const;
};
class SylvExceptionMessage : public SylvException
{
- char message[500];
+ std::string message;
public:
- SylvExceptionMessage(const char *f, int l, const char *mes);
- int printMessage(char *str, int maxlen) const override;
+ SylvExceptionMessage(std::string f, int l, std::string mes);
+ std::string getMessage() const override;
};
// define macros:
diff --git a/dynare++/sylv/cc/SylvMatrix.cc b/dynare++/sylv/cc/SylvMatrix.cc
index 74e777e96b72e1898f54fc81342bd263ce44d024..7744e895c395f782b835a2300559449322101e28 100644
--- a/dynare++/sylv/cc/SylvMatrix.cc
+++ b/dynare++/sylv/cc/SylvMatrix.cc
@@ -8,18 +8,17 @@
#include <dynblas.h>
#include <dynlapack.h>
-#include <cstdio>
#include <cstring>
#include <cmath>
+#include <vector>
void
SylvMatrix::multLeftI(const SqSylvMatrix &m)
{
int off = rows - m.numRows();
if (off < 0)
- {
- throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeftI.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeftI.");
+
GeneralMatrix subtmp(*this, off, 0, m.numRows(), cols);
subtmp.multLeft(m);
}
@@ -29,9 +28,8 @@ SylvMatrix::multLeftITrans(const SqSylvMatrix &m)
{
int off = rows - m.numRows();
if (off < 0)
- {
- throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeftITrans.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeftITrans.");
+
GeneralMatrix subtmp(*this, off, 0, m.numRows(), cols);
subtmp.multLeftTrans(m);
}
@@ -42,9 +40,8 @@ SylvMatrix::multLeft(int zero_cols, const GeneralMatrix &a, const GeneralMatrix
int off = a.numRows() - a.numCols();
if (off < 0 || a.numRows() != rows || off != zero_cols
|| rows != b.numRows() || cols != b.numCols())
- {
- throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeft.");
- }
+ 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)
@@ -67,9 +64,8 @@ void
SylvMatrix::multRightKron(const SqSylvMatrix &m, int order)
{
if (power(m.numRows(), order) != cols)
- {
- throw SYLV_MES_EXCEPTION("Wrong number of cols for right kron multiply.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong number of cols for right kron multiply.");
+
KronVector auxrow(m.numRows(), m.numRows(), order-1);
for (int i = 0; i < rows; i++)
{
@@ -84,9 +80,7 @@ void
SylvMatrix::multRightKronTrans(const SqSylvMatrix &m, int order)
{
if (power(m.numRows(), order) != cols)
- {
- throw SYLV_MES_EXCEPTION("Wrong number of cols for right kron multiply.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong number of cols for right kron multiply.");
KronVector auxrow(m.numRows(), m.numRows(), order-1);
for (int i = 0; i < rows; i++)
@@ -179,9 +173,7 @@ SqSylvMatrix::multVecKron(KronVector &x, const KronVector &d) const
{
x.zeros();
if (d.getDepth() == 0)
- {
- multaVec(x, d);
- }
+ multaVec(x, d);
else
{
KronVector aux(x.getM(), x.getN(), x.getDepth());
@@ -208,9 +200,7 @@ SqSylvMatrix::multVecKronTrans(KronVector &x, const KronVector &d) const
{
x.zeros();
if (d.getDepth() == 0)
- {
- multaVecTrans(x, d);
- }
+ multaVecTrans(x, d);
else
{
KronVector aux(x.getM(), x.getN(), x.getDepth());
@@ -237,51 +227,43 @@ SqSylvMatrix::multInvLeft2(GeneralMatrix &a, GeneralMatrix &b,
double &rcond1, double &rcondinf) const
{
if (rows != a.numRows() || rows != b.numRows())
- {
- throw SYLV_MES_EXCEPTION("Wrong dimensions for multInvLeft2.");
- }
+ throw SYLV_MES_EXCEPTION("Wrong dimensions for multInvLeft2.");
+
// PLU factorization
Vector inv(data);
- auto *const ipiv = new lapack_int[rows];
+ std::vector<lapack_int> ipiv(rows);
lapack_int info;
lapack_int rows2 = rows;
- dgetrf(&rows2, &rows2, inv.base(), &rows2, ipiv, &info);
+ dgetrf(&rows2, &rows2, inv.base(), &rows2, ipiv.data(), &info);
// solve a
lapack_int acols = a.numCols();
double *abase = a.base();
- dgetrs("N", &rows2, &acols, inv.base(), &rows2, ipiv,
+ dgetrs("N", &rows2, &acols, inv.base(), &rows2, ipiv.data(),
abase, &rows2, &info);
// solve b
lapack_int bcols = b.numCols();
double *bbase = b.base();
- dgetrs("N", &rows2, &bcols, inv.base(), &rows2, ipiv,
+ dgetrs("N", &rows2, &bcols, inv.base(), &rows2, ipiv.data(),
bbase, &rows2, &info);
- delete [] ipiv;
// condition numbers
- auto *const work = new double[4*rows];
- auto *const iwork = new lapack_int[rows];
+ std::vector<double> work(4*rows);
+ std::vector<lapack_int> iwork(rows);
double norm1 = getNorm1();
dgecon("1", &rows2, inv.base(), &rows2, &norm1, &rcond1,
- work, iwork, &info);
+ work.data(), iwork.data(), &info);
double norminf = getNormInf();
dgecon("I", &rows2, inv.base(), &rows2, &norminf, &rcondinf,
- work, iwork, &info);
- delete [] iwork;
- delete [] work;
+ work.data(), iwork.data(), &info);
}
void
SqSylvMatrix::setUnit()
{
for (int i = 0; i < rows; i++)
- {
- for (int j = 0; j < cols; j++)
- {
- if (i == j)
- get(i, j) = 1.0;
- else
- get(i, j) = 0.0;
- }
- }
+ for (int j = 0; j < cols; j++)
+ if (i == j)
+ get(i, j) = 1.0;
+ else
+ get(i, j) = 0.0;
}
diff --git a/dynare++/sylv/cc/SylvMatrix.hh b/dynare++/sylv/cc/SylvMatrix.hh
index 9c97913956b4952c4a8f780c2ff3fff254b20903..54a59509b8ea4ad6810c935d0c96af17147d5caa 100644
--- a/dynare++/sylv/cc/SylvMatrix.hh
+++ b/dynare++/sylv/cc/SylvMatrix.hh
@@ -84,7 +84,7 @@ public:
{
}
SqSylvMatrix(const GeneralMatrix &a, const GeneralMatrix &b);
- const SqSylvMatrix &
+ SqSylvMatrix &
operator=(const SqSylvMatrix &m)
{
GeneralMatrix::operator=(m);
diff --git a/dynare++/sylv/cc/SylvMemory.cc b/dynare++/sylv/cc/SylvMemory.cc
index a156637704f9b0594bd42539c1ad6c7aaccc2bc0..fdbaeed395390bbb64efaef0a6d80bc98540cfd6 100644
--- a/dynare++/sylv/cc/SylvMemory.cc
+++ b/dynare++/sylv/cc/SylvMemory.cc
@@ -11,7 +11,6 @@
#endif
#include <cmath>
-#include <cstdio>
#include <cstdlib>
/**********************************************************/
diff --git a/dynare++/sylv/cc/SylvMemory.hh b/dynare++/sylv/cc/SylvMemory.hh
index cd7e7f4a7e876efc6073b962fb85499ce069fb91..da1a8f1811bf6e5948f08b3fb75f904f56846f25 100644
--- a/dynare++/sylv/cc/SylvMemory.hh
+++ b/dynare++/sylv/cc/SylvMemory.hh
@@ -36,7 +36,7 @@ class SylvMemoryPool
public:
SylvMemoryPool() = default;
SylvMemoryPool(const SylvMemoryPool &) = delete;
- const SylvMemoryPool &operator=(const SylvMemoryPool &) = delete;
+ SylvMemoryPool &operator=(const SylvMemoryPool &) = delete;
~SylvMemoryPool();
void init(size_t size);
void *allocate(size_t size);
@@ -51,7 +51,7 @@ public:
SylvMemoryDriver(int num_d, int m, int n, int order);
SylvMemoryDriver(const SylvParams &pars, int num_d, int m, int n, int order);
SylvMemoryDriver(const SylvMemoryDriver &) = delete;
- const SylvMemoryDriver &operator=(const SylvMemoryDriver &) = delete;
+ SylvMemoryDriver &operator=(const SylvMemoryDriver &) = delete;
static void setStackMode(bool);
~SylvMemoryDriver();
protected:
diff --git a/dynare++/sylv/cc/SylvParams.cc b/dynare++/sylv/cc/SylvParams.cc
index 4eb35bddc38c56217669eb6aed5ec5969beedc95..00da72e7b86ba54b98064326b117ed141d4614ec 100644
--- a/dynare++/sylv/cc/SylvParams.cc
+++ b/dynare++/sylv/cc/SylvParams.cc
@@ -4,46 +4,47 @@
#include "SylvParams.hh"
+#include <iostream>
+
void
-SylvParams::print(const char *prefix) const
+SylvParams::print(const std::string &prefix) const
{
- print(stdout, prefix);
+ print(std::cout, prefix);
}
void
-SylvParams::print(FILE *fdesc, const char *prefix) const
+SylvParams::print(std::ostream &fdesc, const std::string &prefix) const
{
- rcondA1.print(fdesc, prefix, "reci. cond1 A ", "%8.4g");
- rcondAI.print(fdesc, prefix, "reci. condInf A ", "%8.4g");
- bs_norm.print(fdesc, prefix, "log10 diag norm ", "%8.4g");
- f_err1.print(fdesc, prefix, "abs. err 1 F diag ", "%8.4g");
- f_errI.print(fdesc, prefix, "abs. err I F diag ", "%8.4g");
- viv_err1.print(fdesc, prefix, "abs. err 1 V*invV ", "%8.4g");
- viv_errI.print(fdesc, prefix, "abs. err I V*invV ", "%8.4g");
- ivv_err1.print(fdesc, prefix, "abs. err 1 invV*V ", "%8.4g");
- ivv_errI.print(fdesc, prefix, "abs. err I invV*V ", "%8.4g");
- f_blocks.print(fdesc, prefix, "num blocks in F ", "%d");
- f_largest.print(fdesc, prefix, "largest block in F ", "%d");
- f_zeros.print(fdesc, prefix, "num zeros in F ", "%d");
- f_offdiag.print(fdesc, prefix, "num offdiag in F ", "%d");
+ rcondA1.print(fdesc, prefix, "reci. cond1 A ");
+ rcondAI.print(fdesc, prefix, "reci. condInf A ");
+ bs_norm.print(fdesc, prefix, "log10 diag norm ");
+ f_err1.print(fdesc, prefix, "abs. err 1 F diag ");
+ f_errI.print(fdesc, prefix, "abs. err I F diag ");
+ viv_err1.print(fdesc, prefix, "abs. err 1 V*invV ");
+ viv_errI.print(fdesc, prefix, "abs. err I V*invV ");
+ ivv_err1.print(fdesc, prefix, "abs. err 1 invV*V ");
+ ivv_errI.print(fdesc, prefix, "abs. err I invV*V ");
+ f_blocks.print(fdesc, prefix, "num blocks in F ");
+ f_largest.print(fdesc, prefix, "largest block in F ");
+ f_zeros.print(fdesc, prefix, "num zeros in F ");
+ f_offdiag.print(fdesc, prefix, "num offdiag in F ");
if (*method == iter)
{
- converged.print(fdesc, prefix, "converged ", "%d");
- convergence_tol.print(fdesc, prefix, "convergence tol. ", "%8.4g");
- iter_last_norm.print(fdesc, prefix, "last norm ", "%8.4g");
- max_num_iter.print(fdesc, prefix, "max num iter ", "%d");
- num_iter.print(fdesc, prefix, "num iter ", "%d");
+ converged.print(fdesc, prefix, "converged ");
+ convergence_tol.print(fdesc, prefix, "convergence tol. ");
+ iter_last_norm.print(fdesc, prefix, "last norm ");
+ max_num_iter.print(fdesc, prefix, "max num iter ");
+ num_iter.print(fdesc, prefix, "num iter ");
}
else
- {
- eig_min.print(fdesc, prefix, "minimum eigenvalue ", "%8.4g");
- }
- mat_err1.print(fdesc, prefix, "rel. matrix norm1 ", "%8.4g");
- mat_errI.print(fdesc, prefix, "rel. matrix normInf", "%8.4g");
- mat_errF.print(fdesc, prefix, "rel. matrix normFro", "%8.4g");
- vec_err1.print(fdesc, prefix, "rel. vector norm1 ", "%8.4g");
- vec_errI.print(fdesc, prefix, "rel. vector normInf", "%8.4g");
- cpu_time.print(fdesc, prefix, "time (CPU secs) ", "%8.4g");
+ eig_min.print(fdesc, prefix, "minimum eigenvalue ");
+
+ mat_err1.print(fdesc, prefix, "rel. matrix norm1 ");
+ mat_errI.print(fdesc, prefix, "rel. matrix normInf");
+ mat_errF.print(fdesc, prefix, "rel. matrix normFro");
+ vec_err1.print(fdesc, prefix, "rel. vector norm1 ");
+ vec_errI.print(fdesc, prefix, "rel. vector normInf");
+ cpu_time.print(fdesc, prefix, "time (CPU secs) ");
}
void
diff --git a/dynare++/sylv/cc/SylvParams.hh b/dynare++/sylv/cc/SylvParams.hh
index 34b836ac60ce6f778dc222acfdc26e284051ec88..ad0a6a88fc6ba9ee76ffc527916bbab1a516d7be 100644
--- a/dynare++/sylv/cc/SylvParams.hh
+++ b/dynare++/sylv/cc/SylvParams.hh
@@ -5,8 +5,8 @@
#ifndef SYLV_PARAMS_H
#define SYLV_PARAMS_H
-#include <cstdio>
-#include <cstring>
+#include <string>
+#include <ostream>
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
# include <dynmex.h>
@@ -34,12 +34,12 @@ public:
{
value = item.value; s = item.s;
}
- const _Self &
+ _Self &
operator=(const _Self &item)
{
value = item.value; s = item.s; return *this;
}
- const _Self &
+ _Self &
operator=(const _Type &val)
{
value = val; s = changed; return *this;
@@ -55,19 +55,14 @@ public:
return s;
}
void
- print(FILE *f, const char *prefix, const char *str, const char *fmt) const
+ print(std::ostream &out, const std::string &prefix, const std::string &str) const
{
if (s == undef)
return;
- char out[1000];
- strcpy(out, prefix);
- strcat(out, str);
- strcat(out, "= ");
- strcat(out, fmt);
+ out << prefix << str << "= " << value;
if (s == def)
- strcat(out, " <default>");
- strcat(out, "\n");
- fprintf(f, out, value);
+ out << " <default>";
+ out << std::endl;
}
};
@@ -86,9 +81,8 @@ protected:
DoubleParamItem(double val) : ParamItem<double>(val)
{
}
- DoubleParamItem(const DoubleParamItem &item)
- = default;
- const DoubleParamItem &
+ DoubleParamItem(const DoubleParamItem &item) = default;
+ DoubleParamItem &
operator=(const double &val)
{
ParamItem<double>::operator=(val); return *this;
@@ -107,9 +101,8 @@ protected:
IntParamItem(int val) : ParamItem<int>(val)
{
}
- IntParamItem(const IntParamItem &item)
- = default;
- const IntParamItem &
+ IntParamItem(const IntParamItem &item) = default;
+ IntParamItem &
operator=(const int &val)
{
ParamItem<int>::operator=(val); return *this;
@@ -128,9 +121,8 @@ protected:
BoolParamItem(bool val) : ParamItem<bool>(val)
{
}
- BoolParamItem(const BoolParamItem &item)
- = default;
- const BoolParamItem &
+ BoolParamItem(const BoolParamItem &item) = default;
+ BoolParamItem &
operator=(const bool &val)
{
ParamItem<bool>::operator=(val); return *this;
@@ -149,9 +141,8 @@ protected:
MethodParamItem(solve_method val) : ParamItem<solve_method>(val)
{
}
- MethodParamItem(const MethodParamItem &item)
- = default;
- const MethodParamItem
+ MethodParamItem(const MethodParamItem &item) = default;
+ MethodParamItem
operator=(const solve_method &val)
{
ParamItem<solve_method>::operator=(val); return *this;
@@ -202,15 +193,14 @@ public:
{
copy(p);
}
- const SylvParams &
+ SylvParams &
operator=(const SylvParams &p)
{
copy(p); return *this;
}
- ~SylvParams()
- = default;
- void print(const char *prefix) const;
- void print(FILE *fdesc, const char *prefix) const;
+ ~SylvParams() = default;
+ void print(const std::string &prefix) const;
+ void print(std::ostream &fdesc, const std::string &prefix) const;
void setArrayNames(int &num, const char **names) const;
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mxArray *createStructArray() const;
diff --git a/dynare++/sylv/cc/SylvesterSolver.hh b/dynare++/sylv/cc/SylvesterSolver.hh
index 9855c6d65369b279eba7d73f04a5d3de4c765622..622dc266b1207840a499ed008bf134b6ebeb2a3a 100644
--- a/dynare++/sylv/cc/SylvesterSolver.hh
+++ b/dynare++/sylv/cc/SylvesterSolver.hh
@@ -12,11 +12,13 @@
#include "SylvParams.hh"
#include "SchurDecomp.hh"
+#include <memory>
+
class SylvesterSolver
{
protected:
- const QuasiTriangular *const matrixK;
- const QuasiTriangular *const matrixF;
+ 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 */
@@ -28,26 +30,25 @@ private:
}
public:
SylvesterSolver(const QuasiTriangular &k, const QuasiTriangular &f)
- : matrixK(new QuasiTriangular(k)),
- matrixF(new QuasiTriangular(f))
+ : matrixK(std::make_unique<QuasiTriangular>(k)),
+ matrixF(std::make_unique<QuasiTriangular>(f))
{
}
SylvesterSolver(const SchurDecompZero &kdecomp, const SchurDecomp &fdecomp)
: matrixK((zeroPad(kdecomp)) ?
- new QuasiTriangular(kdecomp) : new QuasiTriangularZero(kdecomp)),
- matrixF(new QuasiTriangular(fdecomp))
+ std::make_unique<QuasiTriangular>(kdecomp) :
+ std::make_unique<QuasiTriangularZero>(kdecomp)),
+ matrixF(std::make_unique<QuasiTriangular>(fdecomp))
{
}
SylvesterSolver(const SchurDecompZero &kdecomp, const SimilarityDecomp &fdecomp)
: matrixK((zeroPad(kdecomp)) ?
- new QuasiTriangular(kdecomp) : new QuasiTriangularZero(kdecomp)),
- matrixF(new BlockDiagonal(fdecomp.getB()))
- {
- }
- virtual ~SylvesterSolver()
+ std::make_unique<QuasiTriangular>(kdecomp) :
+ std::make_unique<QuasiTriangularZero>(kdecomp)),
+ matrixF(std::make_unique<BlockDiagonal>(fdecomp.getB()))
{
- delete matrixK; delete matrixF;
}
+ virtual ~SylvesterSolver() = default;
virtual void solve(SylvParams &pars, KronVector &x) const = 0;
};
diff --git a/dynare++/sylv/cc/SymSchurDecomp.cc b/dynare++/sylv/cc/SymSchurDecomp.cc
index 746c43465d7c294ee72eb7d616639b5817f10cb0..b11ee363cbe45dc8ec74697cf7a245c97f0e1ef5 100644
--- a/dynare++/sylv/cc/SymSchurDecomp.cc
+++ b/dynare++/sylv/cc/SymSchurDecomp.cc
@@ -9,6 +9,7 @@
#include <algorithm>
#include <cmath>
+#include <vector>
SymSchurDecomp::SymSchurDecomp(const GeneralMatrix &mata)
: lambda(mata.numRows()), q(mata.numRows())
@@ -36,7 +37,7 @@ SymSchurDecomp::SymSchurDecomp(const GeneralMatrix &mata)
double *w = lambda.base();
double *z = q.base();
lapack_int ldz = q.getLD();
- lapack_int *isuppz = new lapack_int[2*std::max(1, (int) m)];
+ std::vector<lapack_int> isuppz(2*std::max(1, (int) m));
double tmpwork;
lapack_int lwork = -1;
lapack_int tmpiwork;
@@ -45,25 +46,21 @@ SymSchurDecomp::SymSchurDecomp(const GeneralMatrix &mata)
// query for lwork and liwork
dsyevr(jobz, range, uplo, &n, a, &lda, vl, vu, il, iu, &abstol,
- &m, w, z, &ldz, isuppz, &tmpwork, &lwork, &tmpiwork, &liwork, &info);
+ &m, w, z, &ldz, isuppz.data(), &tmpwork, &lwork, &tmpiwork, &liwork, &info);
lwork = (int) tmpwork;
liwork = tmpiwork;
// allocate work arrays
- auto *work = new double[lwork];
- auto *iwork = new lapack_int[liwork];
+ std::vector<double> work(lwork);
+ std::vector<lapack_int> iwork(liwork);
// do the calculation
dsyevr(jobz, range, uplo, &n, a, &lda, vl, vu, il, iu, &abstol,
- &m, w, z, &ldz, isuppz, work, &lwork, iwork, &liwork, &info);
+ &m, w, z, &ldz, isuppz.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
if (info < 0)
throw SYLV_MES_EXCEPTION("Internal error in SymSchurDecomp constructor");
if (info > 0)
throw SYLV_MES_EXCEPTION("Internal LAPACK error in DSYEVR");
-
- delete [] work;
- delete [] iwork;
- delete [] isuppz;
}
void
diff --git a/dynare++/sylv/cc/SymSchurDecomp.hh b/dynare++/sylv/cc/SymSchurDecomp.hh
index 2324f08c431ee3fd359085f24a9a760b0682de82..085649d24f84bd1f89f86590275b7d2452bbce62 100644
--- a/dynare++/sylv/cc/SymSchurDecomp.hh
+++ b/dynare++/sylv/cc/SymSchurDecomp.hh
@@ -16,11 +16,8 @@ public:
/** Calculates A = Q*Lambda*Q^T, where A is assummed to be
* symmetric and Lambda real diagonal, hence a vector. */
SymSchurDecomp(const GeneralMatrix &a);
- SymSchurDecomp(const SymSchurDecomp &ssd)
-
- = default;
- virtual ~SymSchurDecomp()
- = default;
+ SymSchurDecomp(const SymSchurDecomp &ssd) = default;
+ virtual ~SymSchurDecomp() = default;
const Vector &
getLambda() const
{
diff --git a/dynare++/sylv/cc/TriangularSylvester.cc b/dynare++/sylv/cc/TriangularSylvester.cc
index e6373395799f6738a542b5d62c68937a624a996e..8ba2e8bc2a0583e252e4e440da5b884b2ca79e04 100644
--- a/dynare++/sylv/cc/TriangularSylvester.cc
+++ b/dynare++/sylv/cc/TriangularSylvester.cc
@@ -7,12 +7,9 @@
#include "KronUtils.hh"
#include "BlockDiagonal.hh"
-#include <cstdio>
+#include <iostream>
#include <cmath>
-double TriangularSylvester::diag_zero = 1.e-15;
-double TriangularSylvester::diag_zero_sq = 1.e-30;
-
TriangularSylvester::TriangularSylvester(const QuasiTriangular &k,
const QuasiTriangular &f)
: SylvesterSolver(k, f),
@@ -40,9 +37,9 @@ TriangularSylvester::TriangularSylvester(const SchurDecompZero &kdecomp,
void
TriangularSylvester::print() const
{
- printf("matrix K (%d):\n", matrixK->getDiagonal().getSize());
+ std::cout << "matrix K (" << matrixK->getDiagonal().getSize() << "):" << std::endl;
matrixK->print();
- printf("matrix F (%d):\n", matrixF->getDiagonal().getSize());
+ std::cout << "matrix F (" << matrixF->getDiagonal().getSize() << "):" << std::endl;
matrixF->print();
}
@@ -67,16 +64,10 @@ TriangularSylvester::solvi(double r, KronVector &d, double &eig_min) const
for (const_diag_iter di = matrixF->diag_begin();
di != matrixF->diag_end();
++di)
- {
- if ((*di).isReal())
- {
- solviRealAndEliminate(r, di, d, eig_min);
- }
- else
- {
- solviComplexAndEliminate(r, di, d, eig_min);
- }
- }
+ if ((*di).isReal())
+ solviRealAndEliminate(r, di, d, eig_min);
+ else
+ solviComplexAndEliminate(r, di, d, eig_min);
}
}
@@ -115,16 +106,10 @@ TriangularSylvester::solviip(double alpha, double betas,
const_diag_iter di = matrixF->diag_begin();
const_diag_iter dsi = matrixFF->diag_begin();
for (; di != matrixF->diag_end(); ++di, ++dsi)
- {
- if ((*di).isReal())
- {
- solviipRealAndEliminate(alpha, betas, di, dsi, d, eig_min);
- }
- else
- {
- solviipComplexAndEliminate(alpha, betas, di, dsi, d, eig_min);
- }
- }
+ if ((*di).isReal())
+ solviipRealAndEliminate(alpha, betas, di, dsi, d, eig_min);
+ else
+ solviipComplexAndEliminate(alpha, betas, di, dsi, d, eig_min);
}
}
@@ -138,9 +123,7 @@ TriangularSylvester::solviRealAndEliminate(double r, const_diag_iter di,
KronVector dj(d, jbar);
// solve system
if (abs(r*f) > diag_zero)
- {
- solvi(r*f, dj, eig_min);
- }
+ solvi(r*f, dj, eig_min);
// calculate y
KronVector y((const KronVector &)dj);
KronUtils::multKron(*matrixF, *matrixK, y);
@@ -177,9 +160,7 @@ TriangularSylvester::solviComplexAndEliminate(double r, const_diag_iter di,
KronVector djj(d, jbar+1);
// solve
if (r*r*aspbs > diag_zero_sq)
- {
- solvii(r*alpha, r*beta1, r*beta2, dj, djj, eig_min);
- }
+ solvii(r*alpha, r*beta1, r*beta2, dj, djj, eig_min);
KronVector y1(dj);
KronVector y2(djj);
KronUtils::multKron(*matrixF, *matrixK, y1);
@@ -219,9 +200,7 @@ TriangularSylvester::solviipRealAndEliminate(double alpha, double betas,
KronVector dj(d, jbar);
// solve
if (fs*aspbs > diag_zero_sq)
- {
- solviip(f*alpha, fs*betas, dj, eig_min);
- }
+ solviip(f*alpha, fs*betas, dj, eig_min);
KronVector y1((const KronVector &)dj);
KronVector y2((const KronVector &)dj);
KronUtils::multKron(*matrixF, *matrixK, y1);
@@ -265,9 +244,7 @@ TriangularSylvester::solviipComplexAndEliminate(double alpha, double betas,
KronVector dj(d, jbar);
KronVector djj(d, jbar+1);
if (gspds*aspbs > diag_zero_sq)
- {
- solviipComplex(alpha, betas, gamma, delta1, delta2, dj, djj, eig_min);
- }
+ solviipComplex(alpha, betas, gamma, delta1, delta2, dj, djj, eig_min);
// here dj, djj is solution, set y1, y2, y11, y22
// y1
KronVector y1((const KronVector &)dj);
@@ -404,9 +381,7 @@ TriangularSylvester::multEigVector(KronVector &eig, const Vector &feig,
int n = eig.getN();
if (depth == 0)
- {
- eig = keig;
- }
+ eig = keig;
else
{
KronVector aux(m, n, depth-1);
diff --git a/dynare++/sylv/cc/TriangularSylvester.hh b/dynare++/sylv/cc/TriangularSylvester.hh
index e804e3343ac8ad157e151899400155626f402bbf..58e0ddde58a4f5de8ff5bb41d1ff6110fa537811 100644
--- a/dynare++/sylv/cc/TriangularSylvester.hh
+++ b/dynare++/sylv/cc/TriangularSylvester.hh
@@ -113,8 +113,8 @@ private:
KronVector &d1, KronVector &d2,
double &eig_min) const;
/* norms for what we consider zero on diagonal of F */
- static double diag_zero;
- static double diag_zero_sq; // square of diag_zero
+ static constexpr double diag_zero = 1.e-15;
+ static constexpr double diag_zero_sq = 1.e-30; // square of diag_zero
};
#endif /* TRIANGULAR_SYLVESTER_H */
diff --git a/dynare++/sylv/cc/Vector.hh b/dynare++/sylv/cc/Vector.hh
index ce08790fd5261409f55c1194afddbdbd994c1ff3..34e2ef42241eb957e074cd75c38de3a3acee2933 100644
--- a/dynare++/sylv/cc/Vector.hh
+++ b/dynare++/sylv/cc/Vector.hh
@@ -133,8 +133,8 @@ public:
}
private:
void copy(const double *d, int inc);
- const Vector &operator=(int); // must not be used (not implemented)
- const Vector &operator=(double); // must not be used (not implemented)
+ Vector &operator=(int); // must not be used (not implemented)
+ Vector &operator=(double); // must not be used (not implemented)
};
class BaseConstVector
diff --git a/dynare++/sylv/matlab/gensylv.cc b/dynare++/sylv/matlab/gensylv.cc
index 113a776dc15f2b19045136dcd0ef5c3f67f87309..2cd7691b20c6e0139198c86de2aac8ba2635823e 100644
--- a/dynare++/sylv/matlab/gensylv.cc
+++ b/dynare++/sylv/matlab/gensylv.cc
@@ -70,23 +70,17 @@ extern "C" {
try
{
if (nlhs == 2)
- {
- gen_sylv_solve(order, n, m, zero_cols,
- mxGetPr(A), mxGetPr(B), mxGetPr(C),
- mxGetPr(X));
- }
+ gen_sylv_solve(order, n, m, zero_cols,
+ mxGetPr(A), mxGetPr(B), mxGetPr(C),
+ mxGetPr(X));
else if (nlhs == 3)
- {
- gen_sylv_solve_and_check(order, n, m, zero_cols,
- mxGetPr(A), mxGetPr(B), mxGetPr(C),
- mxGetPr(D), mxGetPr(X), plhs[2]);
- }
+ gen_sylv_solve_and_check(order, n, m, zero_cols,
+ mxGetPr(A), mxGetPr(B), mxGetPr(C),
+ mxGetPr(D), mxGetPr(X), plhs[2]);
}
catch (const SylvException &e)
{
- char mes[1000];
- e.printMessage(mes, 999);
- DYN_MEX_FUNC_ERR_MSG_TXT(mes);
+ DYN_MEX_FUNC_ERR_MSG_TXT(e.getMessage().c_str());
}
plhs[1] = X;
plhs[0] = mxCreateDoubleScalar(0);