diff --git a/dynare++/sylv/cc/BlockDiagonal.cc b/dynare++/sylv/cc/BlockDiagonal.cc index ecaca4b5a..a2182e792 100644 --- a/dynare++/sylv/cc/BlockDiagonal.cc +++ b/dynare++/sylv/cc/BlockDiagonal.cc @@ -5,12 +5,11 @@ #include "BlockDiagonal.hh" #include -#include #include BlockDiagonal::BlockDiagonal(ConstVector d, int d_size) : QuasiTriangular(std::move(d), d_size), - row_len(new int[d_size]), col_len(new int[d_size]) + row_len(d_size), col_len(d_size) { for (int i = 0; i < d_size; i++) { @@ -21,7 +20,7 @@ BlockDiagonal::BlockDiagonal(ConstVector d, int d_size) BlockDiagonal::BlockDiagonal(const QuasiTriangular &t) : QuasiTriangular(t), - row_len(new int[t.numRows()]), col_len(new int[t.numRows()]) + row_len(t.numRows()), col_len(t.numRows()) { for (int i = 0; i < t.numRows(); i++) { @@ -32,18 +31,8 @@ BlockDiagonal::BlockDiagonal(const QuasiTriangular &t) BlockDiagonal::BlockDiagonal(int p, const BlockDiagonal &b) : QuasiTriangular(p, b), - row_len(new int[b.numRows()]), col_len(new int[b.numRows()]) + row_len(b.row_len), col_len(b.col_len) { - memcpy(row_len, b.row_len, b.numRows()*sizeof(int)); - memcpy(col_len, b.col_len, b.numRows()*sizeof(int)); -} - -BlockDiagonal::BlockDiagonal(const BlockDiagonal &b) - : QuasiTriangular(b), - row_len(new int[b.numRows()]), col_len(new int[b.numRows()]) -{ - memcpy(row_len, b.row_len, b.numRows()*sizeof(int)); - memcpy(col_len, b.col_len, b.numRows()*sizeof(int)); } /* put zeroes to right upper submatrix whose first column is defined @@ -135,9 +124,7 @@ BlockDiagonal::getNumZeros() const { int sum = 0; for (int i = 0; i < diagonal.getSize(); i++) - { - sum += diagonal.getSize() - row_len[i]; - } + sum += diagonal.getSize() - row_len[i]; return sum; } diff --git a/dynare++/sylv/cc/BlockDiagonal.hh b/dynare++/sylv/cc/BlockDiagonal.hh index a2ab6233a..c0cbcc18e 100644 --- a/dynare++/sylv/cc/BlockDiagonal.hh +++ b/dynare++/sylv/cc/BlockDiagonal.hh @@ -6,28 +6,25 @@ #define BLOCK_DIAGONAL_H #include +#include #include "QuasiTriangular.hh" class BlockDiagonal : public QuasiTriangular { - int *const row_len; - int *const col_len; + std::vector row_len, col_len; public: BlockDiagonal(ConstVector d, int d_size); BlockDiagonal(int p, const BlockDiagonal &b); - BlockDiagonal(const BlockDiagonal &b); + BlockDiagonal(const BlockDiagonal &b) = default; BlockDiagonal(const QuasiTriangular &t); - BlockDiagonal & - operator=(const QuasiTriangular &t) + BlockDiagonal &operator=(const QuasiTriangular &t) { - GeneralMatrix::operator=(t); return *this; - } - BlockDiagonal &operator=(const BlockDiagonal &b); - ~BlockDiagonal() override - { - delete [] row_len; delete [] col_len; + GeneralMatrix::operator=(t); + return *this; } + BlockDiagonal &operator=(const BlockDiagonal &b) = default; + ~BlockDiagonal() override = default; void setZeroBlockEdge(diag_iter edge); int getNumZeros() const; int getNumBlocks() const; diff --git a/dynare++/sylv/cc/GeneralMatrix.cc b/dynare++/sylv/cc/GeneralMatrix.cc index 20041cfe4..7b7e43584 100644 --- a/dynare++/sylv/cc/GeneralMatrix.cc +++ b/dynare++/sylv/cc/GeneralMatrix.cc @@ -10,7 +10,6 @@ #include #include -#include #include #include #include @@ -28,7 +27,7 @@ GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &m) copy(m); } -GeneralMatrix::GeneralMatrix(const GeneralMatrix &m, const char *dummy) +GeneralMatrix::GeneralMatrix(const GeneralMatrix &m, const std::string &dummy) : data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols) { for (int i = 0; i < m.rows; i++) @@ -36,7 +35,7 @@ GeneralMatrix::GeneralMatrix(const GeneralMatrix &m, const char *dummy) get(j, i) = m.get(i, j); } -GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &m, const char *dummy) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &m, const std::string &dummy) : data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols) { for (int i = 0; i < m.rows; i++) @@ -61,20 +60,20 @@ GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatr gemm("N", a, "N", b, 1.0, 0.0); } -GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, const char *dum) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, const std::string &dum) : data(a.rows*b.rows), rows(a.rows), cols(b.rows), ld(a.rows) { gemm("N", a, "T", b, 1.0, 0.0); } -GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const char *dum, const ConstGeneralMatrix &b) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const std::string &dum, const ConstGeneralMatrix &b) : data(a.cols*b.cols), rows(a.cols), cols(b.cols), ld(a.cols) { gemm("T", a, "N", b, 1.0, 0.0); } -GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const char *dum1, - const ConstGeneralMatrix &b, const char *dum2) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const std::string &dum1, + const ConstGeneralMatrix &b, const std::string &dum2) : data(a.cols*b.rows), rows(a.cols), cols(b.rows), ld(a.cols) { gemm("T", a, "T", b, 1.0, 0.0); @@ -140,21 +139,21 @@ GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix void GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, - const char *dum, double mult) + const std::string &dum, double mult) { gemm("N", a, "T", b, mult, 1.0); } void -GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const char *dum, +GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const std::string &dum, const ConstGeneralMatrix &b, double mult) { gemm("T", a, "N", b, mult, 1.0); } void -GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const char *dum1, - const ConstGeneralMatrix &b, const char *dum2, double mult) +GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const std::string &dum1, + const ConstGeneralMatrix &b, const std::string &dum2, double mult) { gemm("T", a, "T", b, mult, 1.0); } @@ -269,7 +268,7 @@ GeneralMatrix::add(double a, const ConstGeneralMatrix &m) } void -GeneralMatrix::add(double a, const ConstGeneralMatrix &m, const char *dum) +GeneralMatrix::add(double a, const ConstGeneralMatrix &m, const std::string &dum) { if (m.numRows() != cols || m.numCols() != rows) throw SYLV_MES_EXCEPTION("Matrix has different size in GeneralMatrix::add."); @@ -288,20 +287,20 @@ GeneralMatrix::copy(const ConstGeneralMatrix &m, int ioff, int joff) } void -GeneralMatrix::gemm(const char *transa, const ConstGeneralMatrix &a, - const char *transb, const ConstGeneralMatrix &b, +GeneralMatrix::gemm(const std::string &transa, const ConstGeneralMatrix &a, + const std::string &transb, const ConstGeneralMatrix &b, double alpha, double beta) { int opa_rows = a.numRows(); int opa_cols = a.numCols(); - if (!strcmp(transa, "T")) + if (transa == "T") { opa_rows = a.numCols(); opa_cols = a.numRows(); } int opb_rows = b.numRows(); int opb_cols = b.numCols(); - if (!strcmp(transb, "T")) + if (transb == "T") { opb_rows = b.numCols(); opb_cols = b.numRows(); @@ -322,7 +321,7 @@ GeneralMatrix::gemm(const char *transa, const ConstGeneralMatrix &a, blas_int ldc = ld; if (lda > 0 && ldb > 0 && ldc > 0) { - dgemm(transa, transb, &m, &n, &k, &alpha, a.data.base(), &lda, + dgemm(transa.c_str(), transb.c_str(), &m, &n, &k, &alpha, a.data.base(), &lda, b.data.base(), &ldb, &beta, data.base(), &ldc); } else if (numRows()*numCols() > 0) @@ -335,7 +334,7 @@ GeneralMatrix::gemm(const char *transa, const ConstGeneralMatrix &a, } void -GeneralMatrix::gemm_partial_left(const char *trans, const ConstGeneralMatrix &m, +GeneralMatrix::gemm_partial_left(const std::string &trans, const ConstGeneralMatrix &m, double alpha, double beta) { int icol; @@ -354,7 +353,7 @@ GeneralMatrix::gemm_partial_left(const char *trans, const ConstGeneralMatrix &m, } void -GeneralMatrix::gemm_partial_right(const char *trans, const ConstGeneralMatrix &m, +GeneralMatrix::gemm_partial_right(const std::string &trans, const ConstGeneralMatrix &m, double alpha, double beta) { int irow; @@ -469,7 +468,7 @@ ConstGeneralMatrix::multVecTrans(double a, Vector &x, double b, /* m = inv(this)*m */ void -ConstGeneralMatrix::multInvLeft(const char *trans, int mrows, int mcols, int mld, double *d) const +ConstGeneralMatrix::multInvLeft(const std::string &trans, int mrows, int mcols, int mld, double *d) const { if (rows != cols) throw SYLV_MES_EXCEPTION("The matrix is not square for inversion."); @@ -483,7 +482,7 @@ ConstGeneralMatrix::multInvLeft(const char *trans, int mrows, int mcols, int mld lapack_int info; lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld, lda = inv.ld; dgetrf(&rows2, &rows2, inv.getData().base(), &lda, ipiv.data(), &info); - dgetrs(trans, &rows2, &mcols2, inv.base(), &lda, ipiv.data(), d, + dgetrs(trans.c_str(), &rows2, &mcols2, inv.base(), &lda, ipiv.data(), d, &mld2, &info); } } diff --git a/dynare++/sylv/cc/GeneralMatrix.hh b/dynare++/sylv/cc/GeneralMatrix.hh index b88ac6413..37947dded 100644 --- a/dynare++/sylv/cc/GeneralMatrix.hh +++ b/dynare++/sylv/cc/GeneralMatrix.hh @@ -11,6 +11,7 @@ #include #include #include +#include class GeneralMatrix; @@ -120,7 +121,7 @@ public: virtual void print() const; protected: - void multInvLeft(const char *trans, int mrows, int mcols, int mld, double *d) const; + void multInvLeft(const std::string &trans, int mrows, int mcols, int mld, double *d) const; }; class GeneralMatrix @@ -152,19 +153,19 @@ public: explicit GeneralMatrix(const ConstGeneralMatrix &m); GeneralMatrix(GeneralMatrix &&m) = default; - GeneralMatrix(const GeneralMatrix &m, const char *dummy); // transpose - GeneralMatrix(const ConstGeneralMatrix &m, const char *dummy); // transpose + GeneralMatrix(const GeneralMatrix &m, const std::string &dummy); // transpose + GeneralMatrix(const ConstGeneralMatrix &m, const std::string &dummy); // transpose GeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols); GeneralMatrix(GeneralMatrix &m, int i, int j, int nrows, int ncols); /* this = a*b */ GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b); /* this = a*b' */ - GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, const char *dum); + GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, const std::string &dum); /* this = a'*b */ - GeneralMatrix(const ConstGeneralMatrix &a, const char *dum, const ConstGeneralMatrix &b); + GeneralMatrix(const ConstGeneralMatrix &a, const std::string &dum, const ConstGeneralMatrix &b); /* this = a'*b */ - GeneralMatrix(const ConstGeneralMatrix &a, const char *dum1, - const ConstGeneralMatrix &b, const char *dum2); + GeneralMatrix(const ConstGeneralMatrix &a, const std::string &dum1, + const ConstGeneralMatrix &b, const std::string &dum2); virtual ~GeneralMatrix() = default; GeneralMatrix &operator=(const GeneralMatrix &m) = default; @@ -260,30 +261,30 @@ public: /* this = this + scalar*a*b' */ void multAndAdd(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, - const char *dum, double mult = 1.0); + const std::string &dum, double mult = 1.0); void multAndAdd(const GeneralMatrix &a, const GeneralMatrix &b, - const char *dum, double mult = 1.0) + const std::string &dum, double mult = 1.0) { multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult); } /* this = this + scalar*a'*b */ - void multAndAdd(const ConstGeneralMatrix &a, const char *dum, const ConstGeneralMatrix &b, + void multAndAdd(const ConstGeneralMatrix &a, const std::string &dum, const ConstGeneralMatrix &b, double mult = 1.0); void - multAndAdd(const GeneralMatrix &a, const char *dum, const GeneralMatrix &b, + multAndAdd(const GeneralMatrix &a, const std::string &dum, const GeneralMatrix &b, double mult = 1.0) { multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult); } /* this = this + scalar*a'*b' */ - void multAndAdd(const ConstGeneralMatrix &a, const char *dum1, - const ConstGeneralMatrix &b, const char *dum2, double mult = 1.0); + void multAndAdd(const ConstGeneralMatrix &a, const std::string &dum1, + const ConstGeneralMatrix &b, const std::string &dum2, double mult = 1.0); void - multAndAdd(const GeneralMatrix &a, const char *dum1, - const GeneralMatrix &b, const char *dum2, double mult = 1.0) + multAndAdd(const GeneralMatrix &a, const std::string &dum1, + const GeneralMatrix &b, const std::string &dum2, double mult = 1.0) { multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b), dum2, mult); } @@ -394,9 +395,9 @@ public: } /* this = this + scalar*m' */ - void add(double a, const ConstGeneralMatrix &m, const char *dum); + void add(double a, const ConstGeneralMatrix &m, const std::string &dum); void - add(double a, const GeneralMatrix &m, const char *dum) + add(double a, const GeneralMatrix &m, const std::string &dum) { add(a, ConstGeneralMatrix(m), dum); } @@ -426,12 +427,12 @@ private: copy(ConstGeneralMatrix(m), ioff, joff); } - void gemm(const char *transa, const ConstGeneralMatrix &a, - const char *transb, const ConstGeneralMatrix &b, + void gemm(const std::string &transa, const ConstGeneralMatrix &a, + const std::string &transb, const ConstGeneralMatrix &b, double alpha, double beta); void - gemm(const char *transa, const GeneralMatrix &a, - const char *transb, const GeneralMatrix &b, + gemm(const std::string &transa, const GeneralMatrix &a, + const std::string &transb, const GeneralMatrix &b, double alpha, double beta) { gemm(transa, ConstGeneralMatrix(a), transb, ConstGeneralMatrix(b), @@ -439,20 +440,20 @@ private: } /* this = this * op(m) (without whole copy of this) */ - void gemm_partial_right(const char *trans, const ConstGeneralMatrix &m, + void gemm_partial_right(const std::string &trans, const ConstGeneralMatrix &m, double alpha, double beta); void - gemm_partial_right(const char *trans, const GeneralMatrix &m, + gemm_partial_right(const std::string &trans, const GeneralMatrix &m, double alpha, double beta) { gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta); } /* this = op(m) *this (without whole copy of this) */ - void gemm_partial_left(const char *trans, const ConstGeneralMatrix &m, + void gemm_partial_left(const std::string &trans, const ConstGeneralMatrix &m, double alpha, double beta); void - gemm_partial_left(const char *trans, const GeneralMatrix &m, + gemm_partial_left(const std::string &trans, const GeneralMatrix &m, double alpha, double beta) { gemm_partial_left(trans, ConstGeneralMatrix(m), alpha, beta); diff --git a/dynare++/sylv/cc/GeneralSylvester.cc b/dynare++/sylv/cc/GeneralSylvester.cc index 7cce950be..53ed0953a 100644 --- a/dynare++/sylv/cc/GeneralSylvester.cc +++ b/dynare++/sylv/cc/GeneralSylvester.cc @@ -71,7 +71,7 @@ GeneralSylvester::init() cdecomp = std::make_unique(c.getData(), c.numRows(), *(pars.bs_norm)); cdecomp->check(pars, c); cdecomp->infoToPars(pars); - if (*(pars.method) == SylvParams::recurse) + if (*(pars.method) == SylvParams::solve_method::recurse) sylv = std::make_unique(*bdecomp, *cdecomp); else sylv = std::make_unique(*bdecomp, *cdecomp); diff --git a/dynare++/sylv/cc/KronUtils.cc b/dynare++/sylv/cc/KronUtils.cc index 46ab2dd55..11d542bc1 100644 --- a/dynare++/sylv/cc/KronUtils.cc +++ b/dynare++/sylv/cc/KronUtils.cc @@ -27,9 +27,7 @@ KronUtils::multAtLevel(int level, const QuasiTriangular &t, t.multVec(x, b); } else // 0 < level == depth - { - t.multKron(x); - } + t.multKron(x); } void @@ -55,9 +53,7 @@ KronUtils::multAtLevelTrans(int level, const QuasiTriangular &t, t.multVecTrans(x, b); } else // 0 < level == depth - { - t.multKronTrans(x); - } + t.multKronTrans(x); } void diff --git a/dynare++/sylv/cc/KronVector.cc b/dynare++/sylv/cc/KronVector.cc index 8fc940f2c..5e0f2cd18 100644 --- a/dynare++/sylv/cc/KronVector.cc +++ b/dynare++/sylv/cc/KronVector.cc @@ -37,18 +37,17 @@ KronVector::KronVector(KronVector &v, int i) } KronVector::KronVector(const ConstKronVector &v) - : Vector(v.length()), m(v.getM()), n(v.getN()), depth(v.getDepth()) + : Vector(v), m(v.m), n(v.n), depth(v.depth) { - Vector::operator=(v); } KronVector & KronVector::operator=(const ConstKronVector &v) { Vector::operator=(v); - m = v.getM(); - n = v.getN(); - depth = v.getDepth(); + m = v.m; + n = v.n; + depth = v.depth; return *this; } @@ -62,13 +61,7 @@ KronVector::operator=(const Vector &v) } ConstKronVector::ConstKronVector(const KronVector &v) - : ConstVector(v), m(v.getM()), n(v.getN()), depth(v.getDepth()) -{ -} - -ConstKronVector::ConstKronVector(const ConstKronVector &v) - : ConstVector(power(v.getM(), v.getDepth())*v.getN()), m(v.getM()), n(v.getN()), - depth(v.getDepth()) + : ConstVector(v), m(v.m), n(v.n), depth(v.depth) { } diff --git a/dynare++/sylv/cc/KronVector.hh b/dynare++/sylv/cc/KronVector.hh index 88e005d44..150021e14 100644 --- a/dynare++/sylv/cc/KronVector.hh +++ b/dynare++/sylv/cc/KronVector.hh @@ -11,18 +11,22 @@ class ConstKronVector; class KronVector : public Vector { + friend class ConstKronVector; protected: int m{0}; int n{0}; int depth{0}; public: KronVector() = default; + KronVector(const KronVector &v) = default; + KronVector(KronVector &&v) = default; KronVector(int mm, int nn, int dp); // new instance KronVector(Vector &v, int mm, int nn, int dp); // conversion KronVector(KronVector &, int i); // picks i-th subvector // We don't want implict conversion from ConstKronVector, since it's expensive explicit KronVector(const ConstKronVector &v); // new instance and copy KronVector &operator=(const KronVector &v) = default; + KronVector &operator=(KronVector &&v) = default; KronVector &operator=(const ConstKronVector &v); KronVector &operator=(const Vector &v); int @@ -44,6 +48,7 @@ public: class ConstKronVector : public ConstVector { + friend class KronVector; protected: int m; int n; @@ -51,11 +56,14 @@ protected: public: // Implicit conversion from KronVector is ok, since it's cheap ConstKronVector(const KronVector &v); - ConstKronVector(const ConstKronVector &v); + ConstKronVector(const ConstKronVector &v) = default; + ConstKronVector(ConstKronVector &&v) = default; ConstKronVector(const Vector &v, int mm, int nn, int dp); ConstKronVector(ConstVector v, int mm, int nn, int dp); ConstKronVector(const KronVector &v, int i); ConstKronVector(const ConstKronVector &v, int i); + ConstKronVector &operator=(const ConstKronVector &v) = delete; + ConstKronVector &operator=(ConstKronVector &&v) = delete; int getM() const { diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc index 5f14e2d1b..caa340931 100644 --- a/dynare++/sylv/cc/QuasiTriangular.cc +++ b/dynare++/sylv/cc/QuasiTriangular.cc @@ -12,8 +12,6 @@ #include #include -using namespace std; - double DiagonalBlock::getDeterminant() const { @@ -30,9 +28,9 @@ double DiagonalBlock::getSize() const { if (real) - return abs(*alpha); + return std::abs(*alpha); else - return sqrt(getDeterminant()); + return std::sqrt(getDeterminant()); } // this function makes Diagonal inconsistent, it should only be used @@ -119,28 +117,18 @@ Diagonal::Diagonal(double *data, const Diagonal &d) } } -void -Diagonal::copy(const Diagonal &d) -{ - num_all = d.num_all; - num_real = d.num_real; - blocks = d.blocks; -} - int Diagonal::getNumComplex(const double *data, int d_size) { int num_complex = 0; int in = 1; for (int i = 0; i < d_size-1; i++, in = in + d_size + 1) - { - if (!isZero(data[in])) - { - num_complex++; - if (in < d_size - 2 && !isZero(data[in + d_size +1])) - throw SYLV_MES_EXCEPTION("Matrix is not quasi-triangular"); - } - } + if (!isZero(data[in])) + { + num_complex++; + if (in < d_size - 2 && !isZero(data[in + d_size +1])) + throw SYLV_MES_EXCEPTION("Matrix is not quasi-triangular"); + } return num_complex; } @@ -174,7 +162,7 @@ Diagonal::getEigenValues(Vector &eig) const int d_size = getSize(); if (eig.length() != 2*d_size) { - ostringstream mes; + std::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()); @@ -187,7 +175,7 @@ Diagonal::getEigenValues(Vector &eig) const eig[2*ind+1] = 0.0; else { - double beta = sqrt(b.getSBeta()); + double beta = std::sqrt(b.getSBeta()); eig[2*ind+1] = beta; eig[2*ind+2] = eig[2*ind]; eig[2*ind+3] = -beta; @@ -204,28 +192,28 @@ Diagonal::swapLogically(diag_iter it) diag_iter itp = it; ++itp; - if ((*it).isReal() && !(*itp).isReal()) + if (it->isReal() && !itp->isReal()) { // first is real, second is complex - double *d1 = (*it).alpha.a1; - double *d2 = (*itp).alpha.a1; - double *d3 = (*itp).alpha.a2; + double *d1 = it->alpha.a1; + double *d2 = itp->alpha.a1; + double *d3 = itp->alpha.a2; // swap - DiagonalBlock new_it((*it).jbar, d1, d2); + DiagonalBlock new_it(it->jbar, d1, d2); *it = new_it; - DiagonalBlock new_itp((*itp).jbar+1, d3); + DiagonalBlock new_itp(itp->jbar+1, d3); *itp = new_itp; } - else if (!(*it).isReal() && (*itp).isReal()) + else if (!it->isReal() && itp->isReal()) { // first is complex, second is real - double *d1 = (*it).alpha.a1; - double *d2 = (*it).alpha.a2; - double *d3 = (*itp).alpha.a1; + double *d1 = it->alpha.a1; + double *d2 = it->alpha.a2; + double *d3 = itp->alpha.a1; // swap - DiagonalBlock new_it((*it).jbar, d1); + DiagonalBlock new_it(it->jbar, d1); *it = new_it; - DiagonalBlock new_itp((*itp).jbar-1, d2, d3); + DiagonalBlock new_itp(itp->jbar-1, d2, d3); *itp = new_itp; } } @@ -233,17 +221,17 @@ Diagonal::swapLogically(diag_iter it) void Diagonal::checkConsistency(diag_iter it) { - if (!(*it).isReal() && isZero((*it).getBeta2())) + if (!it->isReal() && isZero(it->getBeta2())) { - (*it).getBeta2() = 0.0; // put exact zero - int jbar = (*it).getIndex(); - double *d2 = (*it).alpha.a2; - (*it).alpha.a2 = (*it).alpha.a1; - (*it).real = true; - (*it).beta1 = nullptr; - (*it).beta2 = nullptr; + it->getBeta2() = 0.0; // put exact zero + int jbar = it->getIndex(); + double *d2 = it->alpha.a2; + it->alpha.a2 = it->alpha.a1; + it->real = true; + it->beta1 = nullptr; + it->beta2 = nullptr; DiagonalBlock b(jbar+1, d2); - blocks.insert((++it).iter(), b); + blocks.insert(++it, b); num_real += 2; num_all++; } @@ -257,7 +245,7 @@ Diagonal::getAverageSize(diag_iter start, diag_iter end) for (diag_iter run = start; run != end; ++run) { num++; - res += (*run).getSize(); + res += run->getSize(); } if (num > 0) res = res/num; @@ -271,7 +259,7 @@ Diagonal::findClosestBlock(diag_iter start, diag_iter end, double a) double minim = 1.0e100; for (diag_iter run = start; run != end; ++run) { - double dist = abs(a - (*run).getSize()); + double dist = std::abs(a - run->getSize()); if (dist < minim) { minim = dist; @@ -288,7 +276,7 @@ Diagonal::findNextLargerBlock(diag_iter start, diag_iter end, double a) double minim = 1.0e100; for (diag_iter run = start; run != end; ++run) { - double dist = (*run).getSize() - a; + double dist = run->getSize() - a; if ((0 <= dist) && (dist < minim)) { minim = dist; @@ -318,7 +306,7 @@ Diagonal::print() const bool Diagonal::isZero(double p) { - return (abs(p) < EPS); + return (std::abs(p) < EPS); } QuasiTriangular::const_col_iter @@ -453,21 +441,15 @@ QuasiTriangular::QuasiTriangular(const SchurDecompZero &decomp) zv.zeros(); // fill right upper part with decomp.getRU() for (int i = 0; i < decomp.getRU().numRows(); i++) - { - for (int j = 0; j < decomp.getRU().numCols(); j++) - { - getData()[(j+decomp.getZeroCols())*decomp.getDim()+i] = decomp.getRU().get(i, j); - } - } + for (int j = 0; j < decomp.getRU().numCols(); j++) + getData()[(j+decomp.getZeroCols())*decomp.getDim()+i] = decomp.getRU().get(i, j); + // fill right lower part with decomp.getT() for (int i = 0; i < decomp.getT().numRows(); i++) - { - for (int j = 0; j < decomp.getT().numCols(); j++) - { - getData()[(j+decomp.getZeroCols())*decomp.getDim()+decomp.getZeroCols()+i] - = decomp.getT().get(i, j); - } - } + for (int j = 0; j < decomp.getT().numCols(); j++) + getData()[(j+decomp.getZeroCols())*decomp.getDim()+decomp.getZeroCols()+i] + = decomp.getT().get(i, j); + // construct diagonal diagonal = Diagonal{getData().base(), decomp.getDim()}; } @@ -487,26 +469,22 @@ QuasiTriangular::setMatrixViaIter(double r, const QuasiTriangular &t) const_diag_iter dir = t.diag_begin(); for (; dil != diag_end(); ++dil, ++dir) { - (*dil).getAlpha() = rr*(*(*dir).getAlpha()); - if (!(*dil).isReal()) + dil->getAlpha() = rr*(*dir->getAlpha()); + if (!dil->isReal()) { - (*dil).getBeta1() = rr*(*dir).getBeta1(); - (*dil).getBeta2() = rr*(*dir).getBeta2(); + dil->getBeta1() = rr*dir->getBeta1(); + dil->getBeta2() = rr*dir->getBeta2(); } col_iter cil = col_begin(*dil); const_col_iter cir = t.col_begin(*dir); for (; cil != col_end(*dil); ++cil, ++cir) - { - if ((*dil).isReal()) - { - *cil = rr*(*cir); - } - else - { - cil.a() = rr*cir.a(); - cil.b() = rr*cir.b(); - } - } + if (dil->isReal()) + *cil = rr*(*cir); + else + { + cil.a() = rr*cir.a(); + cil.b() = rr*cir.b(); + } } } @@ -524,26 +502,22 @@ QuasiTriangular::addMatrixViaIter(double r, const QuasiTriangular &t) const_diag_iter dir = t.diag_begin(); for (; dil != diag_end(); ++dil, ++dir) { - (*dil).getAlpha() = (*(*dil).getAlpha()) + rr*(*(*dir).getAlpha()); - if (!(*dil).isReal()) + dil->getAlpha() = (*dil->getAlpha()) + rr*(*dir->getAlpha()); + if (!dil->isReal()) { - (*dil).getBeta1() += rr*(*dir).getBeta1(); - (*dil).getBeta2() += rr*(*dir).getBeta2(); + dil->getBeta1() += rr*dir->getBeta1(); + dil->getBeta2() += rr*dir->getBeta2(); } col_iter cil = col_begin(*dil); const_col_iter cir = t.col_begin(*dir); for (; cil != col_end(*dil); ++cil, ++cir) - { - if ((*dil).isReal()) - { - *cil += rr*(*cir); - } - else - { - cil.a() += rr*cir.a(); - cil.b() += rr*cir.b(); - } - } + if (dil->isReal()) + *cil += rr*(*cir); + else + { + cil.a() += rr*cir.a(); + cil.b() += rr*cir.b(); + } } } @@ -551,9 +525,7 @@ void QuasiTriangular::addUnit() { for (diag_iter di = diag_begin(); di != diag_end(); ++di) - { - (*di).getAlpha() = *((*di).getAlpha()) + 1.0; - } + di->getAlpha() = *(di->getAlpha()) + 1.0; } void @@ -577,15 +549,13 @@ QuasiTriangular::solvePre(Vector &x, double &eig_min) for (diag_iter di = diag_begin(); di != diag_end(); ++di) { double eig_size; - if (!(*di).isReal()) + if (!di->isReal()) { - eig_size = (*di).getDeterminant(); - eliminateLeft((*di).getIndex()+1, (*di).getIndex(), x); + eig_size = di->getDeterminant(); + eliminateLeft(di->getIndex()+1, di->getIndex(), x); } else - { - eig_size = *(*di).getAlpha()*(*(*di).getAlpha()); - } + eig_size = *di->getAlpha()*(*di->getAlpha()); if (eig_size < eig_min) eig_min = eig_size; } @@ -603,15 +573,13 @@ QuasiTriangular::solvePreTrans(Vector &x, double &eig_min) for (diag_iter di = diag_begin(); di != diag_end(); ++di) { double eig_size; - if (!(*di).isReal()) + if (!di->isReal()) { - eig_size = (*di).getDeterminant(); - eliminateRight((*di).getIndex()+1, (*di).getIndex(), x); + eig_size = di->getDeterminant(); + eliminateRight(di->getIndex()+1, di->getIndex(), x); } else - { - eig_size = *(*di).getAlpha()*(*(*di).getAlpha()); - } + eig_size = *di->getAlpha()*(*di->getAlpha()); if (eig_size < eig_min) eig_min = eig_size; } @@ -632,13 +600,11 @@ QuasiTriangular::multVec(Vector &x, const ConstVector &b) const blas_int incx = x.skip(); dtrmv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx); for (const_diag_iter di = diag_begin(); di != diag_end(); ++di) - { - if (!(*di).isReal()) - { - int jbar = (*di).getIndex(); - x[jbar+1] += (*di).getBeta2()*(b[jbar]); - } - } + if (!di->isReal()) + { + int jbar = di->getIndex(); + x[jbar+1] += di->getBeta2()*(b[jbar]); + } } void @@ -650,13 +616,11 @@ QuasiTriangular::multVecTrans(Vector &x, const ConstVector &b) const blas_int incx = x.skip(); dtrmv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx); for (const_diag_iter di = diag_begin(); di != diag_end(); ++di) - { - if (!(*di).isReal()) - { - int jbar = (*di).getIndex(); - x[jbar] += (*di).getBeta2()*b[jbar+1]; - } - } + if (!di->isReal()) + { + int jbar = di->getIndex(); + x[jbar] += di->getBeta2()*b[jbar+1]; + } } void diff --git a/dynare++/sylv/cc/QuasiTriangular.hh b/dynare++/sylv/cc/QuasiTriangular.hh index f8d3c8e7e..1ba51081f 100644 --- a/dynare++/sylv/cc/QuasiTriangular.hh +++ b/dynare++/sylv/cc/QuasiTriangular.hh @@ -12,8 +12,6 @@ #include #include -using namespace std; - class DiagonalBlock; class Diagonal; class DiagPair @@ -23,19 +21,17 @@ private: double *a2; public: DiagPair() = default; - DiagPair(double *aa1, double *aa2) + DiagPair(double *aa1, double *aa2) : a1{aa1}, a2{aa2} { - a1 = aa1; a2 = aa2; - } - DiagPair(const DiagPair &p) - { - a1 = p.a1; a2 = p.a2; } + DiagPair(const DiagPair &p) = default; DiagPair &operator=(const DiagPair &p) = default; DiagPair & operator=(double v) { - *a1 = v; *a2 = v; return *this; + *a1 = v; + *a2 = v; + return *this; } const double & operator*() const @@ -48,6 +44,9 @@ public: friend class DiagonalBlock; }; +// Stores a diagonal block: either a scalar, or a 2x2 block +/* alpha points to the diagonal element(s); beta1 and beta2 point to the + off-diagonal elements of the 2x2 block */ class DiagonalBlock { private: @@ -57,54 +56,25 @@ private: double *beta1; double *beta2; - void - copy(const DiagonalBlock &b) - { - jbar = b.jbar; - real = b.real; - alpha = b.alpha; - beta1 = b.beta1; - beta2 = b.beta2; - } - public: DiagonalBlock() = default; DiagonalBlock(int jb, bool r, double *a1, double *a2, double *b1, double *b2) - : alpha(a1, a2) + : jbar{jb}, real{r}, alpha{a1, a2}, beta1{b1}, beta2{b2} { - jbar = jb; - real = r; - beta1 = b1; - beta2 = b2; } // construct complex block DiagonalBlock(int jb, double *a1, double *a2) - : alpha(a1, a2) + : jbar{jb}, real{false}, alpha{a1, a2}, beta1{a2-1}, beta2{a1+1} { - jbar = jb; - real = false; - beta1 = a2 - 1; - beta2 = a1 + 1; } // construct real block DiagonalBlock(int jb, double *a1) - : alpha(a1, a1) + : jbar{jb}, real{true}, alpha{a1, a1}, beta1{nullptr}, beta2{nullptr} { - jbar = jb; - real = true; - beta1 = nullptr; - beta2 = nullptr; - } - DiagonalBlock(const DiagonalBlock &b) - { - copy(b); - } - DiagonalBlock & - operator=(const DiagonalBlock &b) - { - copy(b); return *this; } + DiagonalBlock(const DiagonalBlock &b) = default; + DiagonalBlock &operator=(const DiagonalBlock &b) = default; int getIndex() const { @@ -144,76 +114,21 @@ public: friend class Diagonal; }; -template -struct _diag_iter -{ - using _Self = _diag_iter<_Tdiag, _Tblock, _Titer>; - _Tdiag diag; - _Titer it; -public: - _diag_iter(_Tdiag d, _Titer iter) : diag(d), it(iter) - { - } - _Tblock - operator*() const - { - return *it; - } - _Self & - operator++() - { - ++it; return *this; - } - _Self & - operator--() - { - --it; return *this; - } - bool - operator==(const _Self &x) const - { - return x.it == it; - } - bool - operator!=(const _Self &x) const - { - return x.it != it; - } - _Self & - operator=(const _Self &x) - { - it = x.it; return *this; - } - _Titer - iter() const - { - return it; - } -}; - class Diagonal { public: - using const_diag_iter = _diag_iter::const_iterator>; - using diag_iter = _diag_iter::iterator>; + using const_diag_iter = std::list::const_iterator; + using diag_iter = std::list::iterator; private: int num_all{0}; - list blocks; + std::list blocks; int num_real{0}; - void copy(const Diagonal &); public: Diagonal() = default; Diagonal(double *data, int d_size); Diagonal(double *data, const Diagonal &d); - Diagonal(const Diagonal &d) - { - copy(d); - } - Diagonal & - operator=(const Diagonal &d) - { - copy(d); return *this; - } + Diagonal(const Diagonal &d) = default; + Diagonal &operator=(const Diagonal &d) = default; virtual ~Diagonal() = default; int @@ -247,22 +162,22 @@ public: diag_iter begin() { - return diag_iter(*this, blocks.begin()); + return blocks.begin(); } const_diag_iter begin() const { - return const_diag_iter(*this, blocks.begin()); + return blocks.begin(); } diag_iter end() { - return diag_iter(*this, blocks.end()); + return blocks.end(); } const_diag_iter end() const { - return const_diag_iter(*this, blocks.end()); + return blocks.end(); } /* redefine pointers as data start at p */ @@ -283,14 +198,11 @@ struct _matrix_iter public: _matrix_iter(_TPtr base, int ds, bool r) { - ptr = base; d_size = ds; real = r; + ptr = base; + d_size = ds; + real = r; } virtual ~_matrix_iter() = default; - _Self & - operator=(const _Self &it) - { - ptr = it.ptr; d_size = it.d_size; real = it.real; return *this; - } bool operator==(const _Self &it) const { @@ -328,19 +240,17 @@ public: _Self & operator++() override { - _Tparent::ptr++; row++; return *this; + _Tparent::ptr++; + row++; + return *this; } _TRef b() const { if (_Tparent::real) - { - return *(_Tparent::ptr); - } + return *(_Tparent::ptr); else - { - return *(_Tparent::ptr+_Tparent::d_size); - } + return *(_Tparent::ptr+_Tparent::d_size); } int getRow() const @@ -363,19 +273,17 @@ public: _Self & operator++() override { - _Tparent::ptr += _Tparent::d_size; col++; return *this; + _Tparent::ptr += _Tparent::d_size; + col++; + return *this; } virtual _TRef b() const { if (_Tparent::real) - { - return *(_Tparent::ptr); - } + return *(_Tparent::ptr); else - { - return *(_Tparent::ptr+1); - } + return *(_Tparent::ptr+1); } int getCol() const @@ -387,6 +295,12 @@ public: class SchurDecomp; class SchurDecompZero; +/* Represents an upper quasi-triangular matrix. + All the elements are stored in the SqSylvMatrix super-class. + Additionally, a list of the diagonal blocks (1x1 or 2x2), is stored in the + "diagonal" member, in order to optimize some operations (where the matrix is + seen as an upper-triangular matrix, plus sub-diagonal elements of the 2x2 + diagonal blocks) */ class QuasiTriangular : public SqSylvMatrix { public: diff --git a/dynare++/sylv/cc/SchurDecompEig.cc b/dynare++/sylv/cc/SchurDecompEig.cc index 5193ed7bf..962a3d3c0 100644 --- a/dynare++/sylv/cc/SchurDecompEig.cc +++ b/dynare++/sylv/cc/SchurDecompEig.cc @@ -19,9 +19,7 @@ SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to) { diag_iter runm = run; if (!tryToSwap(run, runm) && runm == to) - { - ++to; - } + ++to; else { // bubble all eigenvalues from runm(incl.) to run(excl.), diff --git a/dynare++/sylv/cc/SimilarityDecomp.cc b/dynare++/sylv/cc/SimilarityDecomp.cc index ec16e76d7..8857f85e3 100644 --- a/dynare++/sylv/cc/SimilarityDecomp.cc +++ b/dynare++/sylv/cc/SimilarityDecomp.cc @@ -120,9 +120,7 @@ SimilarityDecomp::diagonalize(double norm) ++end; } else - { - bringGuiltyBlock(start, end); // moves with end - } + bringGuiltyBlock(start, end); // moves with end } } diff --git a/dynare++/sylv/cc/SylvMatrix.cc b/dynare++/sylv/cc/SylvMatrix.cc index 1e0b3241b..770916e8f 100644 --- a/dynare++/sylv/cc/SylvMatrix.cc +++ b/dynare++/sylv/cc/SylvMatrix.cc @@ -8,7 +8,6 @@ #include #include -#include #include #include @@ -102,9 +101,7 @@ SylvMatrix::eliminateLeft(int row, int col, Vector &x) get(row, col) = 0.0; double mult = e/d; for (int i = col + 1; i < numCols(); i++) - { - get(row, i) = get(row, i) - mult*get(col, i); - } + get(row, i) = get(row, i) - mult*get(col, i); x[row] = x[row] - mult*x[col]; } else if (std::abs(e) > std::abs(d)) @@ -137,9 +134,7 @@ SylvMatrix::eliminateRight(int row, int col, Vector &x) get(row, col) = 0.0; double mult = e/d; for (int i = 0; i < row; i++) - { - get(i, col) = get(i, col) - mult*get(i, row); - } + get(i, col) = get(i, col) - mult*get(i, row); x[col] = x[col] - mult*x[row]; } else if (std::abs(e) > std::abs(d)) diff --git a/dynare++/sylv/cc/SylvMatrix.hh b/dynare++/sylv/cc/SylvMatrix.hh index 5e96a034a..55993ef2e 100644 --- a/dynare++/sylv/cc/SylvMatrix.hh +++ b/dynare++/sylv/cc/SylvMatrix.hh @@ -39,6 +39,10 @@ public: : GeneralMatrix(a, b) { } + SylvMatrix(const SylvMatrix &m) = default; + SylvMatrix(SylvMatrix &&m) = default; + SylvMatrix &operator=(const SylvMatrix &m) = default; + SylvMatrix &operator=(SylvMatrix &&m) = default; /* this = |I 0|* this |0 m| */ @@ -70,6 +74,7 @@ public: { } SqSylvMatrix(const SqSylvMatrix &m) = default; + SqSylvMatrix(SqSylvMatrix &&m) = default; SqSylvMatrix(const GeneralMatrix &m, int i, int j, int nrows) : SylvMatrix(m, i, j, nrows, nrows) { @@ -79,12 +84,8 @@ public: { } SqSylvMatrix(const GeneralMatrix &a, const GeneralMatrix &b); - SqSylvMatrix & - operator=(const SqSylvMatrix &m) - { - GeneralMatrix::operator=(m); - return *this; - } + SqSylvMatrix &operator=(const SqSylvMatrix &m) = default; + SqSylvMatrix &operator=(SqSylvMatrix &&m) = default; /* x = (this \otimes this..\otimes this)*d */ void multVecKron(KronVector &x, const ConstKronVector &d) const; /* x = (this' \otimes this'..\otimes this')*d */ diff --git a/dynare++/sylv/cc/SylvParams.cc b/dynare++/sylv/cc/SylvParams.cc index 623845005..0053177ba 100644 --- a/dynare++/sylv/cc/SylvParams.cc +++ b/dynare++/sylv/cc/SylvParams.cc @@ -29,7 +29,7 @@ SylvParams::print(std::ostream &fdesc, const std::string &prefix) const 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) + if (*method == solve_method::iter) { converged.print(fdesc, prefix, "converged "); convergence_tol.print(fdesc, prefix, "convergence tol. "); @@ -52,57 +52,57 @@ void SylvParams::setArrayNames(int &num, const char **names) const { num = 0; - if (method.getStatus() != undef) + if (method.getStatus() != status::undef) names[num++] = "method"; - if (convergence_tol.getStatus() != undef) + if (convergence_tol.getStatus() != status::undef) names[num++] = "convergence_tol"; - if (max_num_iter.getStatus() != undef) + if (max_num_iter.getStatus() != status::undef) names[num++] = "max_num_iter"; - if (bs_norm.getStatus() != undef) + if (bs_norm.getStatus() != status::undef) names[num++] = "bs_norm"; - if (converged.getStatus() != undef) + if (converged.getStatus() != status::undef) names[num++] = "converged"; - if (iter_last_norm.getStatus() != undef) + if (iter_last_norm.getStatus() != status::undef) names[num++] = "iter_last_norm"; - if (num_iter.getStatus() != undef) + if (num_iter.getStatus() != status::undef) names[num++] = "num_iter"; - if (f_err1.getStatus() != undef) + if (f_err1.getStatus() != status::undef) names[num++] = "f_err1"; - if (f_errI.getStatus() != undef) + if (f_errI.getStatus() != status::undef) names[num++] = "f_errI"; - if (viv_err1.getStatus() != undef) + if (viv_err1.getStatus() != status::undef) names[num++] = "viv_err1"; - if (viv_errI.getStatus() != undef) + if (viv_errI.getStatus() != status::undef) names[num++] = "viv_errI"; - if (ivv_err1.getStatus() != undef) + if (ivv_err1.getStatus() != status::undef) names[num++] = "ivv_err1"; - if (ivv_errI.getStatus() != undef) + if (ivv_errI.getStatus() != status::undef) names[num++] = "ivv_errI"; - if (f_blocks.getStatus() != undef) + if (f_blocks.getStatus() != status::undef) names[num++] = "f_blocks"; - if (f_largest.getStatus() != undef) + if (f_largest.getStatus() != status::undef) names[num++] = "f_largest"; - if (f_zeros.getStatus() != undef) + if (f_zeros.getStatus() != status::undef) names[num++] = "f_zeros"; - if (f_offdiag.getStatus() != undef) + if (f_offdiag.getStatus() != status::undef) names[num++] = "f_offdiag"; - if (rcondA1.getStatus() != undef) + if (rcondA1.getStatus() != status::undef) names[num++] = "rcondA1"; - if (rcondAI.getStatus() != undef) + if (rcondAI.getStatus() != status::undef) names[num++] = "rcondAI"; - if (eig_min.getStatus() != undef) + if (eig_min.getStatus() != status::undef) names[num++] = "eig_min"; - if (mat_err1.getStatus() != undef) + if (mat_err1.getStatus() != status::undef) names[num++] = "mat_err1"; - if (mat_errI.getStatus() != undef) + if (mat_errI.getStatus() != status::undef) names[num++] = "mat_errI"; - if (mat_errF.getStatus() != undef) + if (mat_errF.getStatus() != status::undef) names[num++] = "mat_errF"; - if (vec_err1.getStatus() != undef) + if (vec_err1.getStatus() != status::undef) names[num++] = "vec_err1"; - if (vec_errI.getStatus() != undef) + if (vec_errI.getStatus() != status::undef) names[num++] = "vec_errI"; - if (cpu_time.getStatus() != undef) + if (cpu_time.getStatus() != status::undef) names[num++] = "cpu_time"; } @@ -133,7 +133,7 @@ SylvParams::BoolParamItem::createMatlabArray() const mxArray * SylvParams::MethodParamItem::createMatlabArray() const { - if (value == iter) + if (value == solve_method::iter) return mxCreateString("iterative"); else return mxCreateString("recursive"); @@ -149,57 +149,57 @@ SylvParams::createStructArray() const mxArray *const res = mxCreateStructArray(2, dims, num, names); int i = 0; - if (method.getStatus() != undef) + if (method.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, method.createMatlabArray()); - if (convergence_tol.getStatus() != undef) + if (convergence_tol.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, convergence_tol.createMatlabArray()); - if (max_num_iter.getStatus() != undef) + if (max_num_iter.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, max_num_iter.createMatlabArray()); - if (bs_norm.getStatus() != undef) + if (bs_norm.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, bs_norm.createMatlabArray()); - if (converged.getStatus() != undef) + if (converged.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, converged.createMatlabArray()); - if (iter_last_norm.getStatus() != undef) + if (iter_last_norm.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, iter_last_norm.createMatlabArray()); - if (num_iter.getStatus() != undef) + if (num_iter.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, num_iter.createMatlabArray()); - if (f_err1.getStatus() != undef) + if (f_err1.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, f_err1.createMatlabArray()); - if (f_errI.getStatus() != undef) + if (f_errI.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, f_errI.createMatlabArray()); - if (viv_err1.getStatus() != undef) + if (viv_err1.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, viv_err1.createMatlabArray()); - if (viv_errI.getStatus() != undef) + if (viv_errI.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, viv_errI.createMatlabArray()); - if (ivv_err1.getStatus() != undef) + if (ivv_err1.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, ivv_err1.createMatlabArray()); - if (ivv_errI.getStatus() != undef) + if (ivv_errI.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, ivv_errI.createMatlabArray()); - if (f_blocks.getStatus() != undef) + if (f_blocks.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, f_blocks.createMatlabArray()); - if (f_largest.getStatus() != undef) + if (f_largest.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, f_largest.createMatlabArray()); - if (f_zeros.getStatus() != undef) + if (f_zeros.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, f_zeros.createMatlabArray()); - if (f_offdiag.getStatus() != undef) + if (f_offdiag.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, f_offdiag.createMatlabArray()); - if (rcondA1.getStatus() != undef) + if (rcondA1.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, rcondA1.createMatlabArray()); - if (rcondAI.getStatus() != undef) + if (rcondAI.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, rcondAI.createMatlabArray()); - if (eig_min.getStatus() != undef) + if (eig_min.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, eig_min.createMatlabArray()); - if (mat_err1.getStatus() != undef) + if (mat_err1.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, mat_err1.createMatlabArray()); - if (mat_errI.getStatus() != undef) + if (mat_errI.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, mat_errI.createMatlabArray()); - if (mat_errF.getStatus() != undef) + if (mat_errF.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, mat_errF.createMatlabArray()); - if (vec_err1.getStatus() != undef) + if (vec_err1.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, vec_err1.createMatlabArray()); - if (vec_errI.getStatus() != undef) + if (vec_errI.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, vec_errI.createMatlabArray()); - if (cpu_time.getStatus() != undef) + if (cpu_time.getStatus() != status::undef) mxSetFieldByNumber(res, 0, i++, cpu_time.createMatlabArray()); return res; diff --git a/dynare++/sylv/cc/SylvParams.hh b/dynare++/sylv/cc/SylvParams.hh index a183697a8..58520b1e2 100644 --- a/dynare++/sylv/cc/SylvParams.hh +++ b/dynare++/sylv/cc/SylvParams.hh @@ -12,7 +12,7 @@ # include #endif -typedef enum {def, changed, undef} status; +enum class status {def, changed, undef}; template struct ParamItem @@ -24,25 +24,21 @@ protected: public: ParamItem() { - s = undef; + s = status::undef; } ParamItem(_Type val) { - value = val; s = def; - } - ParamItem(const _Self &item) - { - value = item.value; s = item.s; - } - _Self & - operator=(const _Self &item) - { - value = item.value; s = item.s; return *this; + value = val; + s = status::def; } + ParamItem(const _Self &item) = default; + _Self &operator=(const _Self &item) = default; _Self & operator=(const _Type &val) { - value = val; s = changed; return *this; + value = val; + s = status::changed; + return *this; } _Type operator*() const @@ -57,10 +53,10 @@ public: void print(std::ostream &out, const std::string &prefix, const std::string &str) const { - if (s == undef) + if (s == status::undef) return; out << prefix << str << "= " << value; - if (s == def) + if (s == status::def) out << " "; out << std::endl; } @@ -69,7 +65,7 @@ public: class SylvParams { public: - using solve_method = enum {iter, recurse}; + enum class solve_method {iter, recurse}; protected: class DoubleParamItem : public ParamItem @@ -184,7 +180,7 @@ public: DoubleParamItem cpu_time; // time of the job in CPU seconds SylvParams(bool wc = false) - : method(recurse), convergence_tol(1.e-30), max_num_iter(15), + : method(solve_method::recurse), convergence_tol(1.e-30), max_num_iter(15), bs_norm(1.3), want_check(wc) { } @@ -206,10 +202,10 @@ operator<<(std::ostream &out, SylvParams::solve_method m) { switch (m) { - case SylvParams::iter: + case SylvParams::solve_method::iter: out << "iterative"; break; - case SylvParams::recurse: + case SylvParams::solve_method::recurse: out << "recurse (a.k.a. triangular)"; break; } diff --git a/dynare++/sylv/cc/SymSchurDecomp.cc b/dynare++/sylv/cc/SymSchurDecomp.cc index c0f281a08..236c9e15f 100644 --- a/dynare++/sylv/cc/SymSchurDecomp.cc +++ b/dynare++/sylv/cc/SymSchurDecomp.cc @@ -19,9 +19,6 @@ SymSchurDecomp::SymSchurDecomp(const ConstGeneralMatrix &mata) throw SYLV_MES_EXCEPTION("Matrix is not square in SymSchurDecomp constructor"); // prepare for dsyevr - const char *jobz = "V"; - const char *range = "A"; - const char *uplo = "U"; lapack_int n = mata.numRows(); GeneralMatrix tmpa(mata); double *a = tmpa.base(); @@ -45,7 +42,7 @@ SymSchurDecomp::SymSchurDecomp(const ConstGeneralMatrix &mata) lapack_int info; // query for lwork and liwork - dsyevr(jobz, range, uplo, &n, a, &lda, vl, vu, il, iu, &abstol, + dsyevr("V", "A", "U", &n, a, &lda, vl, vu, il, iu, &abstol, &m, w, z, &ldz, isuppz.data(), &tmpwork, &lwork, &tmpiwork, &liwork, &info); lwork = (int) tmpwork; liwork = tmpiwork; @@ -54,7 +51,7 @@ SymSchurDecomp::SymSchurDecomp(const ConstGeneralMatrix &mata) std::vector iwork(liwork); // do the calculation - dsyevr(jobz, range, uplo, &n, a, &lda, vl, vu, il, iu, &abstol, + dsyevr("V", "A", "U", &n, a, &lda, vl, vu, il, iu, &abstol, &m, w, z, &ldz, isuppz.data(), work.data(), &lwork, iwork.data(), &liwork, &info); if (info < 0) diff --git a/dynare++/sylv/cc/TriangularSylvester.cc b/dynare++/sylv/cc/TriangularSylvester.cc index 3092b9d1f..bc899d778 100644 --- a/dynare++/sylv/cc/TriangularSylvester.cc +++ b/dynare++/sylv/cc/TriangularSylvester.cc @@ -48,7 +48,7 @@ TriangularSylvester::solve(SylvParams &pars, KronVector &d) const { double eig_min = 1e30; solvi(1., d, eig_min); - pars.eig_min = sqrt(eig_min); + pars.eig_min = std::sqrt(eig_min); } void @@ -122,7 +122,7 @@ TriangularSylvester::solviRealAndEliminate(double r, const_diag_iter di, double f = *((*di).getAlpha()); KronVector dj(d, jbar); // solve system - if (abs(r*f) > diag_zero) + if (std::abs(r*f) > diag_zero) solvi(r*f, dj, eig_min); // calculate y KronVector y((const KronVector &)dj); @@ -277,8 +277,8 @@ TriangularSylvester::solviipComplex(double alpha, double betas, double gamma, KronVector d2tmp(d2); quaEval(alpha, betas, gamma, delta1, delta2, d1, d2, d1tmp, d2tmp); - double delta = sqrt(delta1*delta2); - double beta = sqrt(betas); + double delta = std::sqrt(delta1*delta2); + double beta = std::sqrt(betas); double a1 = alpha*gamma - beta*delta; double b1 = alpha*delta + gamma*beta; double a2 = alpha*gamma + beta*delta; diff --git a/dynare++/sylv/cc/Vector.cc b/dynare++/sylv/cc/Vector.cc index c902c52f6..f826fc960 100644 --- a/dynare++/sylv/cc/Vector.cc +++ b/dynare++/sylv/cc/Vector.cc @@ -14,8 +14,6 @@ #include #include -using namespace std; - ZeroPad zero_pad; Vector::Vector(const Vector &v) @@ -25,9 +23,9 @@ Vector::Vector(const Vector &v) } Vector::Vector(const ConstVector &v) - : len(v.length()), data{new double[len], [](double *arr) { delete[] arr; }} + : len(v.len), data{new double[len], [](double *arr) { delete[] arr; }} { - copy(v.base(), v.skip()); + copy(v.base(), v.s); } Vector & @@ -38,17 +36,13 @@ Vector::operator=(const Vector &v) if (v.len != len) throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths."); - /* + if (s == v.s - && (data <= v.data && v.data < data+len*s - || v.data <= data && data < v.data+v.len*v.s) - && (data-v.data) % s == 0) - { - std::cout << "this destroy=" << destroy << ", v destroy=" << v.destroy - << ", data-v.data=" << (unsigned long) (data-v.data) - << ", len=" << len << std::endl; - throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors."); - } */ + && (base() <= v.base() && v.base() < base()+len*s + || v.base() <= base() && base() < v.base()+v.len*v.s) + && (base()-v.base()) % s == 0) + throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors."); + copy(v.base(), v.s); return *this; } @@ -65,16 +59,15 @@ Vector::operator=(Vector &&v) Vector & Vector::operator=(const ConstVector &v) { - if (v.length() != len) + if (v.len != len) throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths."); - /* - if (v.skip() == 1 && skip() == 1 && ( - (base() < v.base() + v.length() && base() >= v.base()) - || (base() + length() < v.base() + v.length() - && base() + length() > v.base()))) + if (s == v.s + && (base() <= v.base() && v.base() < base()+len*s + || v.base() <= base() && base() < v.base()+v.len*v.s) + && (base()-v.base()) % s == 0) throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors."); - */ - copy(v.base(), v.skip()); + + copy(v.base(), v.s); return *this; } @@ -210,7 +203,7 @@ void Vector::add(double r, const ConstVector &v) { blas_int n = len; - blas_int incx = v.skip(); + blas_int incx = v.s; blas_int incy = s; daxpy(&n, &r, v.base(), &incx, base(), &incy); } @@ -225,7 +218,7 @@ void Vector::addComplex(const std::complex &z, const ConstVector &v) { blas_int n = len/2; - blas_int incx = v.skip(); + blas_int incx = v.s; blas_int incy = s; zaxpy(&n, reinterpret_cast(z), v.base(), &incx, base(), &incy); } @@ -374,8 +367,8 @@ ConstVector::getMax() const { double r = 0; for (int i = 0; i < len; i++) - if (abs(operator[](i)) > r) - r = abs(operator[](i)); + if (std::abs(operator[](i)) > r) + r = std::abs(operator[](i)); return r; } @@ -384,7 +377,7 @@ ConstVector::getNorm1() const { double norm = 0.0; for (int i = 0; i < len; i++) - norm += abs(operator[](i)); + norm += std::abs(operator[](i)); return norm; } @@ -403,7 +396,7 @@ bool ConstVector::isFinite() const { int i = 0; - while (i < len && isfinite(operator[](i))) + while (i < len && std::isfinite(operator[](i))) i++; return i == len; } diff --git a/dynare++/sylv/cc/Vector.hh b/dynare++/sylv/cc/Vector.hh index 025dc4c56..eca0f4bc9 100644 --- a/dynare++/sylv/cc/Vector.hh +++ b/dynare++/sylv/cc/Vector.hh @@ -140,6 +140,7 @@ class ConstGeneralMatrix; class ConstVector { + friend class Vector; protected: int len; int off{0}; // offset to double* pointer diff --git a/dynare++/sylv/testing/MMMatrix.cc b/dynare++/sylv/testing/MMMatrix.cc index 9d3493d99..e9efa49e5 100644 --- a/dynare++/sylv/testing/MMMatrix.cc +++ b/dynare++/sylv/testing/MMMatrix.cc @@ -4,70 +4,52 @@ #include "MMMatrix.hh" -#include -#include +#include +#include -MMMatrixIn::MMMatrixIn(const char *fname) +MMMatrixIn::MMMatrixIn(const std::string &fname) { - FILE *fd; - if (nullptr == (fd = fopen(fname, "r"))) - throw MMException(string("Cannot open file ")+fname+" for reading\n"); + std::ifstream fd{fname}; + if (fd.fail()) + throw MMException("Cannot open file "+fname+" for reading\n"); - char buffer[1000]; // jump over initial comments - while (fgets(buffer, 1000, fd) && strncmp(buffer, "%%", 2)) - { - } + while (fd.peek() == '%') + fd.ignore(std::numeric_limits::max(), '\n'); + // read in number of rows and cols - if (!fgets(buffer, 1000, fd)) - throw MMException(string("Cannot read rows and cols while reading ")+fname+"\n"); - if (2 != sscanf(buffer, "%d %d", &rows, &cols)) + fd >> rows >> cols; + if (fd.fail()) throw MMException("Couldn't parse rows and cols\n"); // read in data data = std::shared_ptr(static_cast(operator new[](rows*cols*sizeof(double))), [](double *arr) { operator delete[](static_cast(arr)); }); int len = rows*cols; int i = 0; - while (fgets(buffer, 1000, fd) && i < len) + while (!fd.eof() && i < len) { - if (1 != sscanf(buffer, "%lf", const_cast(data.get())+i)) - throw MMException(string("Couldn't parse float number ")+buffer+"\n"); + fd >> data.get()[i]; + if (fd.fail()) + throw MMException("Couldn't parse float number\n"); i++; } if (i < len) - { - char mes[1000]; - sprintf(mes, "Couldn't read all %d lines, read %d so far\n", len, i); - throw MMException(mes); - } - fclose(fd); + throw MMException("Couldn't read all " + std::to_string(len) + " elements, read " + + std::to_string(i) + " so far\n"); + fd.close(); } void -MMMatrixOut::write(const char *fname, int rows, int cols, const double *data) +MMMatrixOut::write(const std::string &fname, const GeneralMatrix &m) { - FILE *fd; - if (nullptr == (fd = fopen(fname, "w"))) - throw MMException(string("Cannot open file ")+fname+" for writing\n"); + std::ofstream fd{fname, std::ios::out | std::ios::trunc}; + if (fd.fail()) + throw MMException("Cannot open file "+fname+" for writing\n"); - if (0 > fprintf(fd, "%%%%MatrixMarket matrix array real general\n")) - throw MMException(string("Output error when writing file ")+fname); - if (0 > fprintf(fd, "%d %d\n", rows, cols)) - throw MMException(string("Output error when writing file ")+fname); - int running = 0; - for (int i = 0; i < cols; i++) - { - for (int j = 0; j < rows; j++) - { - if (0 > fprintf(fd, "%40.35g\n", data[running])) - throw MMException(string("Output error when writing file ")+fname); - running++; - } - } - fclose(fd); -} - -void -MMMatrixOut::write(const char *fname, const GeneralMatrix &m) -{ - write(fname, m.numRows(), m.numCols(), m.base()); + fd << "%%%%MatrixMarket matrix array real general" << std::endl + << m.numRows() << ' ' << m.numCols() << std::endl + << std::setprecision(35); + for (int i = 0; i < m.numCols(); i++) + for (int j = 0; j < m.numRows(); j++) + fd << std::setw(40) << m.get(i, j); + fd.close(); } diff --git a/dynare++/sylv/testing/MMMatrix.hh b/dynare++/sylv/testing/MMMatrix.hh index 3c6bb0e14..1b9c7243b 100644 --- a/dynare++/sylv/testing/MMMatrix.hh +++ b/dynare++/sylv/testing/MMMatrix.hh @@ -12,22 +12,17 @@ #include #include -using namespace std; - class MMException : public MallocAllocator { - string message; + std::string message; public: - MMException(string mes) : message(std::move(mes)) + MMException(std::string mes) : message(std::move(mes)) { } - MMException(const char *mes) : message(mes) - { - } - const char * + std::string getMessage() const { - return message.data(); + return message; } }; @@ -37,7 +32,7 @@ class MMMatrixIn : public MallocAllocator int rows; int cols; public: - MMMatrixIn(const char *fname); + MMMatrixIn(const std::string &fname); ~MMMatrixIn() = default; Vector getData() const @@ -64,8 +59,7 @@ public: class MMMatrixOut : public MallocAllocator { public: - static void write(const char *fname, int rows, int cols, const double *data); - static void write(const char *fname, const GeneralMatrix &m); + static void write(const std::string &fname, const GeneralMatrix &m); }; #endif /* MM_MATRIX_H */ diff --git a/dynare++/sylv/testing/tests.cc b/dynare++/sylv/testing/tests.cc index 2e768a864..54ac2c3dc 100644 --- a/dynare++/sylv/testing/tests.cc +++ b/dynare++/sylv/testing/tests.cc @@ -18,76 +18,65 @@ #include "MMMatrix.hh" -#include -#include #include - #include +#include +#include +#include +#include +#include class TestRunnable : public MallocAllocator { - char name[100]; - static double eps_norm; public: - TestRunnable(const char *n) + const std::string name; + static constexpr double eps_norm = 1.0e-10; + TestRunnable(std::string n) : name(std::move(n)) { - strncpy(name, n, 100); } virtual ~TestRunnable() = default; bool test() const; virtual bool run() const = 0; - const char * - getName() const - { - return name; - } protected: // declaration of auxiliary static methods - static bool quasi_solve(bool trans, const char *mname, const char *vname); - static bool mult_kron(bool trans, const char *mname, const char *vname, - const char *cname, int m, int n, int depth); - static bool level_kron(bool trans, const char *mname, const char *vname, - const char *cname, int level, int m, int n, int depth); - static bool kron_power(const char *m1name, const char *m2name, const char *vname, - const char *cname, int m, int n, int depth); - static bool lin_eval(const char *m1name, const char *m2name, const char *vname, - const char *cname, int m, int n, int depth, + static bool quasi_solve(bool trans, const std::string &mname, const std::string &vname); + static bool mult_kron(bool trans, const std::string &mname, const std::string &vname, + const std::string &cname, int m, int n, int depth); + static bool level_kron(bool trans, const std::string &mname, const std::string &vname, + const std::string &cname, int level, int m, int n, int depth); + static bool kron_power(const std::string &m1name, const std::string &m2name, const std::string &vname, + const std::string &cname, int m, int n, int depth); + static bool lin_eval(const std::string &m1name, const std::string &m2name, const std::string &vname, + const std::string &cname, int m, int n, int depth, double alpha, double beta1, double beta2); - static bool qua_eval(const char *m1name, const char *m2name, const char *vname, - const char *cname, int m, int n, int depth, + static bool qua_eval(const std::string &m1name, const std::string &m2name, const std::string &vname, + const std::string &cname, int m, int n, int depth, double alpha, double betas, double gamma, double delta1, double delta2); - static bool tri_sylv(const char *m1name, const char *m2name, const char *vname, + static bool tri_sylv(const std::string &m1name, const std::string &m2name, const std::string &vname, int m, int n, int depth); - static bool gen_sylv(const char *aname, const char *bname, const char *cname, - const char *dname, int m, int n, int order); - static bool eig_bubble(const char *aname, int from, int to); - static bool block_diag(const char *aname, double log10norm = 3.0); - static bool iter_sylv(const char *m1name, const char *m2name, const char *vname, + static bool gen_sylv(const std::string &aname, const std::string &bname, const std::string &cname, + const std::string &dname, int m, int n, int order); + static bool eig_bubble(const std::string &aname, int from, int to); + static bool block_diag(const std::string &aname, double log10norm = 3.0); + static bool iter_sylv(const std::string &m1name, const std::string &m2name, const std::string &vname, int m, int n, int depth); }; -double TestRunnable::eps_norm = 1.0e-10; - bool TestRunnable::test() const { - printf("Running test <%s>\n", name); + std::cout << "Running test <" << name << '>' << std::endl; clock_t start = clock(); bool passed = run(); clock_t end = clock(); - printf("CPU time %8.4g (CPU seconds)..................", - ((double) (end-start))/CLOCKS_PER_SEC); + std::cout << "CPU time " << ((double) (end-start))/CLOCKS_PER_SEC << " (CPU seconds).................."; if (passed) - { - printf("passed\n\n"); - return passed; - } + std::cout << "passed"; else - { - printf("FAILED\n\n"); - return passed; - } + std::cout << "FAILED"; + std::cout << std::endl << std::endl; + return passed; } /**********************************************************/ @@ -95,27 +84,27 @@ TestRunnable::test() const /**********************************************************/ bool -TestRunnable::quasi_solve(bool trans, const char *mname, const char *vname) +TestRunnable::quasi_solve(bool trans, const std::string &mname, const std::string &vname) { MMMatrixIn mmt(mname); MMMatrixIn mmv(vname); SylvMemoryDriver memdriver(1, mmt.row(), mmt.row(), 1); - QuasiTriangular *t; - QuasiTriangular *tsave; + std::unique_ptr t; + std::unique_ptr tsave; if (mmt.row() == mmt.col()) { - t = new QuasiTriangular(mmt.getData(), mmt.row()); - tsave = new QuasiTriangular(*t); + t = std::make_unique(mmt.getData(), mmt.row()); + tsave = std::make_unique(*t); } else if (mmt.row() > mmt.col()) { - t = new QuasiTriangularZero(mmt.row()-mmt.col(), mmt.getData(), mmt.col()); - tsave = new QuasiTriangularZero((const QuasiTriangularZero &) *t); + t = std::make_unique(mmt.row()-mmt.col(), mmt.getData(), mmt.col()); + tsave = std::make_unique((const QuasiTriangularZero &) *t); } else { - printf(" Wrong quasi triangular dimensions, rows must be >= cols.\n"); + std::cout << " Wrong quasi triangular dimensions, rows must be >= cols.\n"; return false; } ConstVector v{mmv.getData()}; @@ -125,24 +114,22 @@ TestRunnable::quasi_solve(bool trans, const char *mname, const char *vname) t->solveTrans(x, v, eig_min); else t->solve(x, v, eig_min); - printf("eig_min = %8.4g\n", eig_min); + std::cout << "eig_min = " << eig_min << std::endl; Vector xx(v.length()); if (trans) tsave->multVecTrans(xx, ConstVector(x)); else tsave->multVec(xx, ConstVector(x)); - delete tsave; - delete t; xx.add(-1.0, v); xx.add(1.0, x); double norm = xx.getNorm(); - printf("\terror norm = %8.4g\n", norm); + std::cout << "\terror norm = " << norm << std::endl; return (norm < eps_norm); } bool -TestRunnable::mult_kron(bool trans, const char *mname, const char *vname, - const char *cname, int m, int n, int depth) +TestRunnable::mult_kron(bool trans, const std::string &mname, const std::string &vname, + const std::string &cname, int m, int n, int depth) { MMMatrixIn mmt(mname); MMMatrixIn mmv(vname); @@ -153,7 +140,10 @@ TestRunnable::mult_kron(bool trans, const char *mname, const char *vname, || mmv.row() != length || mmc.row() != length) { - printf(" Incompatible sizes for krom mult action, len=%d, matrow=%d, m=%d, vrow=%d, crow=%d \n", length, mmt.row(), m, mmv.row(), mmc.row()); + std::cout << " Incompatible sizes for kron mult action, len=" << length + << ", matrow=" << mmt.row() << ", m=" << m + << ", vrow=" << mmv.row() << ", crow=" << mmc.row() + << std::endl; return false; } @@ -169,13 +159,13 @@ TestRunnable::mult_kron(bool trans, const char *mname, const char *vname, t.multKron(v); c.add(-1.0, v); double norm = c.getNorm(); - printf("\terror norm = %8.4g\n", norm); + std::cout << "\terror norm = " << norm << std::endl; return (norm < eps_norm); } bool -TestRunnable::level_kron(bool trans, const char *mname, const char *vname, - const char *cname, int level, int m, int n, int depth) +TestRunnable::level_kron(bool trans, const std::string &mname, const std::string &vname, + const std::string &cname, int level, int m, int n, int depth) { MMMatrixIn mmt(mname); MMMatrixIn mmv(vname); @@ -187,7 +177,10 @@ TestRunnable::level_kron(bool trans, const char *mname, const char *vname, || mmv.row() != length || mmc.row() != length) { - printf(" Incompatible sizes for krom mult action, len=%d, matrow=%d, m=%d, n=%d, vrow=%d, crow=%d \n", length, mmt.row(), m, n, mmv.row(), mmc.row()); + std::cout << " Incompatible sizes for kron mult action, len=" << length + << ", matrow=" << mmt.row() << ", m=" << m << ", n=" << n + << ", vrow=" << mmv.row() << ", crow=" << mmc.row() + << std::endl; return false; } @@ -204,13 +197,13 @@ TestRunnable::level_kron(bool trans, const char *mname, const char *vname, KronUtils::multAtLevel(level, t, x); x.add(-1, c); double norm = x.getNorm(); - printf("\terror norm = %8.4g\n", norm); + std::cout << "\terror norm = " << norm << std::endl; return (norm < eps_norm); } bool -TestRunnable::kron_power(const char *m1name, const char *m2name, const char *vname, - const char *cname, int m, int n, int depth) +TestRunnable::kron_power(const std::string &m1name, const std::string &m2name, const std::string &vname, + const std::string &cname, int m, int n, int depth) { MMMatrixIn mmt1(m1name); MMMatrixIn mmt2(m2name); @@ -223,7 +216,11 @@ TestRunnable::kron_power(const char *m1name, const char *m2name, const char *vna || mmv.row() != length || mmc.row() != length) { - printf(" Incompatible sizes for krom power mult action, len=%d, row1=%d, row2=%d, m=%d, n=%d, vrow=%d, crow=%d \n", length, mmt1.row(), mmt2.row(), m, n, mmv.row(), mmc.row()); + std::cout << " Incompatible sizes for kron power mult action, len=" << length + << ", row1=" << mmt1.row() << ", row2=" << mmt2.row() + << ", m=" << m << ", n=" << n + << ", vrow=" << mmv.row() << ", crow=" << mmc.row() + << std::endl; return false; } @@ -240,13 +237,13 @@ TestRunnable::kron_power(const char *m1name, const char *m2name, const char *vna memdriver.setStackMode(false); x.add(-1, c); double norm = x.getNorm(); - printf("\terror norm = %8.4g\n", norm); + std::cout << "\terror norm = " << norm << std::endl; return (norm < eps_norm); } bool -TestRunnable::lin_eval(const char *m1name, const char *m2name, const char *vname, - const char *cname, int m, int n, int depth, +TestRunnable::lin_eval(const std::string &m1name, const std::string &m2name, const std::string &vname, + const std::string &cname, int m, int n, int depth, double alpha, double beta1, double beta2) { MMMatrixIn mmt1(m1name); @@ -260,7 +257,11 @@ TestRunnable::lin_eval(const char *m1name, const char *m2name, const char *vname || mmv.row() != 2*length || mmc.row() != 2*length) { - printf(" Incompatible sizes for lin eval action, len=%d, row1=%d, row2=%d, m=%d, n=%d, vrow=%d, crow=%d \n", length, mmt1.row(), mmt2.row(), m, n, mmv.row(), mmc.row()); + std::cout << " Incompatible sizes for lin eval action, len=" << length + << ", row1=" << mmt1.row() << ", row2=" << mmt2.row() + << ", m=" << m << ", n=" << n + << ", vrow=" << mmv.row() << ", crow=" << mmc.row() + << std::endl; return false; } @@ -285,13 +286,13 @@ TestRunnable::lin_eval(const char *m1name, const char *m2name, const char *vname x2.add(-1, c2); double norm1 = x1.getNorm(); double norm2 = x2.getNorm(); - printf("\terror norm1 = %8.4g\n\terror norm2 = %8.4g\n", norm1, norm2); + std::cout << "\terror norm1 = " << norm1 << "\n\terror norm2 = " << norm2 << '\n'; return (norm1*norm1+norm2*norm2 < eps_norm*eps_norm); } bool -TestRunnable::qua_eval(const char *m1name, const char *m2name, const char *vname, - const char *cname, int m, int n, int depth, +TestRunnable::qua_eval(const std::string &m1name, const std::string &m2name, const std::string &vname, + const std::string &cname, int m, int n, int depth, double alpha, double betas, double gamma, double delta1, double delta2) { @@ -306,7 +307,11 @@ TestRunnable::qua_eval(const char *m1name, const char *m2name, const char *vname || mmv.row() != 2*length || mmc.row() != 2*length) { - printf(" Incompatible sizes for qua eval action, len=%d, row1=%d, row2=%d, m=%d, n=%d, vrow=%d, crow=%d \n", length, mmt1.row(), mmt2.row(), m, n, mmv.row(), mmc.row()); + std::cout << " Incompatible sizes for qua eval action, len=" << length + << ", row1=" << mmt1.row() << ", row2=" << mmt2.row() + << ", m=" << m << ", n=" << n + << ", vrow=" << mmv.row() << ", crow=" << mmc.row() + << std::endl; return false; } @@ -331,12 +336,12 @@ TestRunnable::qua_eval(const char *m1name, const char *m2name, const char *vname x2.add(-1, c2); double norm1 = x1.getNorm(); double norm2 = x2.getNorm(); - printf("\terror norm1 = %8.4g\n\terror norm2 = %8.4g\n", norm1, norm2); + std::cout << "\terror norm1 = " << norm1 << "\n\terror norm2 = " << norm2 << std::endl; return (norm1*norm1+norm2*norm2 < 100*eps_norm*eps_norm); // relax norm } bool -TestRunnable::tri_sylv(const char *m1name, const char *m2name, const char *vname, +TestRunnable::tri_sylv(const std::string &m1name, const std::string &m2name, const std::string &vname, int m, int n, int depth) { MMMatrixIn mmt1(m1name); @@ -348,7 +353,11 @@ TestRunnable::tri_sylv(const char *m1name, const char *m2name, const char *vname || mmt2.row() != n || mmv.row() != length) { - printf(" Incompatible sizes for triangular sylvester action, len=%d, row1=%d, row2=%d, m=%d, n=%d, vrow=%d\n", length, mmt1.row(), mmt2.row(), m, n, mmv.row()); + std::cout << " Incompatible sizes for triangular sylvester action, len=" << length + << ", row1=" << mmt1.row() << ", row2=" << mmt2.row() + << ", m=" << m << ", n=" << n + << ", vrow=" << mmv.row() + << std::endl; return false; } @@ -369,17 +378,17 @@ TestRunnable::tri_sylv(const char *m1name, const char *m2name, const char *vname dcheck.add(-1.0, v); double norm = dcheck.getNorm(); double xnorm = v.getNorm(); - printf("\trel. error norm = %8.4g\n", norm/xnorm); + std::cout << "\trel. error norm = " << norm/xnorm << std::endl; double max = dcheck.getMax(); double xmax = v.getMax(); - printf("\trel. error max = %8.4g\n", max/xmax); + std::cout << "\trel. error max = " << max/xmax << std::endl; memdriver.setStackMode(false); return (norm < xnorm*eps_norm); } bool -TestRunnable::gen_sylv(const char *aname, const char *bname, const char *cname, - const char *dname, int m, int n, int order) +TestRunnable::gen_sylv(const std::string &aname, const std::string &bname, const std::string &cname, + const std::string &dname, int m, int n, int order) { MMMatrixIn mma(aname); MMMatrixIn mmb(bname); @@ -391,7 +400,7 @@ TestRunnable::gen_sylv(const char *aname, const char *bname, const char *cname, || n != mmb.row() || n < mmb.col() || n != mmd.row() || power(m, order) != mmd.col()) { - printf(" Incompatible sizes for gen_sylv.\n"); + std::cout << " Incompatible sizes for gen_sylv.\n"; return false; } @@ -410,13 +419,13 @@ TestRunnable::gen_sylv(const char *aname, const char *bname, const char *cname, } bool -TestRunnable::eig_bubble(const char *aname, int from, int to) +TestRunnable::eig_bubble(const std::string &aname, int from, int to) { MMMatrixIn mma(aname); if (mma.row() != mma.col()) { - printf(" Matrix is not square\n"); + std::cout << " Matrix is not square\n"; return false; } @@ -438,21 +447,21 @@ TestRunnable::eig_bubble(const char *aname, int from, int to) double normInf = check.getNormInf(); double onorm1 = orig.getNorm1(); double onormInf = orig.getNormInf(); - printf("\tabs. error1 = %8.4g\n", norm1); - printf("\tabs. errorI = %8.4g\n", normInf); - printf("\trel. error1 = %8.4g\n", norm1/onorm1); - printf("\trel. errorI = %8.4g\n", normInf/onormInf); + std:: cout << "\tabs. error1 = " << norm1 << std::endl + << "\tabs. errorI = " << normInf << std::endl + << "\trel. error1 = " << norm1/onorm1 << std::endl + << "\trel. errorI = " << normInf/onormInf << std::endl; return (norm1 < eps_norm*onorm1 && normInf < eps_norm*onormInf); } bool -TestRunnable::block_diag(const char *aname, double log10norm) +TestRunnable::block_diag(const std::string &aname, double log10norm) { MMMatrixIn mma(aname); if (mma.row() != mma.col()) { - printf(" Matrix is not square\n"); + std::cout << " Matrix is not square\n"; return false; } @@ -468,25 +477,25 @@ TestRunnable::block_diag(const char *aname, double log10norm) double normInf = check.getNormInf(); double onorm1 = orig.getNorm1(); double onormInf = orig.getNormInf(); - printf("\terror Q*B*invQ:\n"); - printf("\tabs. error1 = %8.4g\n", norm1); - printf("\tabs. errorI = %8.4g\n", normInf); - printf("\trel. error1 = %8.4g\n", norm1/onorm1); - printf("\trel. errorI = %8.4g\n", normInf/onormInf); + std::cout << "\terror Q*B*invQ:" << std::endl + << "\tabs. error1 = " << norm1 << std::endl + << "\tabs. errorI = " << normInf << std::endl + << "\trel. error1 = " << norm1/onorm1 << std::endl + << "\trel. errorI = " << normInf/onormInf << std::endl; SqSylvMatrix check2(dec.getQ(), dec.getInvQ()); SqSylvMatrix in(n); in.setUnit(); check2.add(-1, in); double nor1 = check2.getNorm1(); double norInf = check2.getNormInf(); - printf("\terror Q*invQ:\n"); - printf("\tabs. error1 = %8.4g\n", nor1); - printf("\tabs. errorI = %8.4g\n", norInf); + std::cout << "\terror Q*invQ:" << std::endl + << "\tabs. error1 = " << nor1 << std::endl + << "\tabs. errorI = " << norInf << std::endl; return (norm1 < eps_norm*pow(10, log10norm)*onorm1); } bool -TestRunnable::iter_sylv(const char *m1name, const char *m2name, const char *vname, +TestRunnable::iter_sylv(const std::string &m1name, const std::string &m2name, const std::string &vname, int m, int n, int depth) { MMMatrixIn mmt1(m1name); @@ -498,7 +507,11 @@ TestRunnable::iter_sylv(const char *m1name, const char *m2name, const char *vnam || mmt2.row() != n || mmv.row() != length) { - printf(" Incompatible sizes for triangular sylvester iteration, len=%d, row1=%d, row2=%d, m=%d, n=%d, vrow=%d\n", length, mmt1.row(), mmt2.row(), m, n, mmv.row()); + std::cout << " Incompatible sizes for triangular sylvester iteration, len=" << length + << ", row1=" << mmt1.row() << ", row2=" << mmt2.row() + << ", m=" << m << ", n=" << n + << ", vrow=" << mmv.row() + << std::endl; return false; } @@ -511,7 +524,7 @@ TestRunnable::iter_sylv(const char *m1name, const char *m2name, const char *vnam ConstKronVector v(vraw, m, n, depth); KronVector d(v); // copy of v SylvParams pars; - pars.method = SylvParams::iter; + pars.method = SylvParams::solve_method::iter; is.solve(pars, d); pars.print("\t"); KronVector dcheck((const KronVector &)d); @@ -520,10 +533,10 @@ TestRunnable::iter_sylv(const char *m1name, const char *m2name, const char *vnam dcheck.add(-1.0, v); double cnorm = dcheck.getNorm(); double xnorm = v.getNorm(); - printf("\trel. error norm = %8.4g\n", cnorm/xnorm); + std::cout << "\trel. error norm = " << cnorm/xnorm << std::endl; double max = dcheck.getMax(); double xmax = v.getMax(); - printf("\trel. error max = %8.4g\n", max/xmax); + std::cout << "\trel. error max = " << max/xmax << std::endl; memdriver.setStackMode(false); return (cnorm < xnorm*eps_norm); } @@ -1149,79 +1162,76 @@ BlockDiagBigTest::run() const int main() { - TestRunnable *all_tests[50]; + std::vector> all_tests; // fill in vector of all tests - int num_tests = 0; - all_tests[num_tests++] = new PureTriangTest(); - all_tests[num_tests++] = new PureTriangTransTest(); - all_tests[num_tests++] = new PureTrLargeTest(); - all_tests[num_tests++] = new PureTrLargeTransTest(); - all_tests[num_tests++] = new QuasiTriangTest(); - all_tests[num_tests++] = new QuasiTriangTransTest(); - all_tests[num_tests++] = new QuasiTrLargeTest(); - all_tests[num_tests++] = new QuasiTrLargeTransTest(); - all_tests[num_tests++] = new QuasiZeroSmallTest(); - all_tests[num_tests++] = new MultKronSmallTest(); - all_tests[num_tests++] = new MultKronTest(); - all_tests[num_tests++] = new MultKronSmallTransTest(); - all_tests[num_tests++] = new MultKronTransTest(); - all_tests[num_tests++] = new LevelKronTest(); - all_tests[num_tests++] = new LevelKronTransTest(); - all_tests[num_tests++] = new LevelZeroKronTest(); - all_tests[num_tests++] = new LevelZeroKronTransTest(); - all_tests[num_tests++] = new KronPowerTest(); - all_tests[num_tests++] = new SmallLinEvalTest(); - all_tests[num_tests++] = new LinEvalTest(); - all_tests[num_tests++] = new SmallQuaEvalTest(); - all_tests[num_tests++] = new QuaEvalTest(); - all_tests[num_tests++] = new EigBubFrankTest(); - all_tests[num_tests++] = new EigBubSplitTest(); - all_tests[num_tests++] = new EigBubSameTest(); - all_tests[num_tests++] = new BlockDiagSmallTest(); - all_tests[num_tests++] = new BlockDiagFrankTest(); - all_tests[num_tests++] = new BlockDiagIllCondTest(); - all_tests[num_tests++] = new BlockDiagBigTest(); - all_tests[num_tests++] = new TriSylvSmallRealTest(); - all_tests[num_tests++] = new TriSylvSmallComplexTest(); - all_tests[num_tests++] = new TriSylvTest(); - all_tests[num_tests++] = new TriSylvBigTest(); - all_tests[num_tests++] = new TriSylvLargeTest(); - all_tests[num_tests++] = new IterSylvTest(); - all_tests[num_tests++] = new IterSylvLargeTest(); - all_tests[num_tests++] = new GenSylvSmallTest(); - all_tests[num_tests++] = new GenSylvTest(); - all_tests[num_tests++] = new GenSylvSingTest(); - all_tests[num_tests++] = new GenSylvLargeTest(); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); // launch the tests + std::cout << std::setprecision(4); int success = 0; - for (int i = 0; i < num_tests; i++) + for (const auto &test : all_tests) { try { - if (all_tests[i]->test()) + if (test->test()) success++; } catch (const MMException &e) { - printf("Caugth MM exception in <%s>:\n%s", all_tests[i]->getName(), - e.getMessage()); + std::cout << "Caught MM exception in <" << test->name << ">:\n" << e.getMessage(); } catch (SylvException &e) { - printf("Caught Sylv exception in %s:\n", all_tests[i]->getName()); + std::cout << "Caught Sylv exception in " << test->name << ":\n"; e.printMessage(); } } - printf("There were %d tests that failed out of %d tests run.\n", - num_tests - success, num_tests); + int nfailed = all_tests.size() - success; + std::cout << "There were " << nfailed << " tests that failed out of " + << all_tests.size() << " tests run." << std::endl; - // destroy - for (int i = 0; i < num_tests; i++) - { - delete all_tests[i]; - } - - return 0; + if (nfailed) + return EXIT_FAILURE; + else + return EXIT_SUCCESS; }