diff --git a/dynare++/extern/matlab/dynare_simul.cc b/dynare++/extern/matlab/dynare_simul.cc index 843c38816..e81b70982 100644 --- a/dynare++/extern/matlab/dynare_simul.cc +++ b/dynare++/extern/matlab/dynare_simul.cc @@ -100,7 +100,7 @@ extern "C" { DYN_MEX_FUNC_ERR_MSG_TXT(buf); } ft.zeros(); - ConstTwoDMatrix gk_mat(ft.nrows(), ft.ncols(), mxGetPr(gk)); + ConstTwoDMatrix gk_mat(ft.nrows(), ft.ncols(), ConstVector{gk}); ft.add(1.0, gk_mat); auto *ut = new UFSTensor(ft); pol.insert(ut); @@ -108,17 +108,16 @@ extern "C" { // form the decision rule UnfoldDecisionRule dr(pol, PartitionY(nstat, npred, nboth, nforw), - nexog, ConstVector(mxGetPr(ysteady), ny)); + nexog, ConstVector{ysteady}); // form the shock realization - TwoDMatrix shocks_mat(nexog, nper, (const double *)mxGetPr(shocks)); - TwoDMatrix vcov_mat(nexog, nexog, (const double *)mxGetPr(vcov)); + TwoDMatrix shocks_mat(nexog, nper, ConstVector{shocks}); + TwoDMatrix vcov_mat(nexog, nexog, ConstVector{vcov}); GenShockRealization sr(vcov_mat, shocks_mat, seed); // simulate and copy the results - Vector ystart_vec((const double *)mxGetPr(ystart), ny); TwoDMatrix *res_mat = dr.simulate(DecisionRule::horner, nper, - ystart_vec, sr); - TwoDMatrix res_tmp_mat(ny, nper, mxGetPr(res)); + ConstVector{ystart}, sr); + TwoDMatrix res_tmp_mat{ny, nper, Vector{res}}; res_tmp_mat = (const TwoDMatrix &) (*res_mat); delete res_mat; plhs[1] = res; diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc index 41085d8e6..ddf9a5058 100644 --- a/dynare++/kord/approximation.cc +++ b/dynare++/kord/approximation.cc @@ -144,7 +144,7 @@ Approximation::walkStochSteady() to |ss|. */ model.solveDeterministicSteady(); approxAtSteady(); - Vector steady0(ss, 0); + Vector steady0{ss.getCol(0)}; steady0 = (const Vector &) model.getSteady(); double sigma_so_far = 0.0; @@ -172,7 +172,7 @@ Approximation::walkStochSteady() rec << " Not converged!!" << endrec; KORD_RAISE_X("Fix point calculation not converged", KORD_FP_NOT_CONV); } - Vector steadyi(ss, i); + Vector steadyi{ss.getCol(i)}; steadyi = (const Vector &) model.getSteady(); // calculate |hh| as expectations of the last $g^{**}$ @@ -378,7 +378,7 @@ Approximation::calcYCov() const TwoDMatrix *X = new TwoDMatrix(guSigma, guTrans); GeneralSylvester gs(1, model.numeq(), model.numeq(), 0, - A.base(), B.base(), C.base(), X->base()); + A.getData(), B.getData(), C.getData(), X->getData()); gs.solve(); return X; diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc index de2f54e1b..e1998d0d7 100644 --- a/dynare++/kord/decision_rule.cc +++ b/dynare++/kord/decision_rule.cc @@ -181,7 +181,7 @@ SimResultsStats::writeMat(mat_t *fd, const char *lname) const { char tmp[100]; sprintf(tmp, "%s_mean", lname); - ConstTwoDMatrix m(num_y, 1, mean.base()); + ConstTwoDMatrix m(num_y, 1, mean); m.writeMat(fd, tmp); sprintf(tmp, "%s_vcov", lname); ConstTwoDMatrix(vcov).writeMat(fd, tmp); @@ -198,7 +198,7 @@ SimResultsStats::calcMean() { for (int j = 0; j < num_per; j++) { - ConstVector col(*i, j); + ConstVector col{i->getCol(j)}; mean.add(mult, col); } } @@ -273,10 +273,10 @@ SimResultsDynamicStats::calcMean() double mult = 1.0/data.size(); for (int j = 0; j < num_per; j++) { - Vector meanj(mean, j); + Vector meanj{mean.getCol(j)}; for (auto & i : data) { - ConstVector col(*i, j); + ConstVector col{i->getCol(j)}; meanj.add(mult, col); } } @@ -292,11 +292,11 @@ SimResultsDynamicStats::calcVariance() double mult = 1.0/(data.size()-1); for (int j = 0; j < num_per; j++) { - ConstVector meanj(mean, j); - Vector varj(variance, j); + ConstVector meanj{mean.getCol(j)}; + Vector varj{variance.getCol(j)}; for (auto & i : data) { - Vector col(ConstVector((*i), j)); + Vector col{i->getCol(j)}; col.add(-1.0, meanj); for (int k = 0; k < col.length(); k++) col[k] = col[k]*col[k]; @@ -426,7 +426,7 @@ RTSimResultsStats::writeMat(mat_t *fd, const char *lname) { char tmp[100]; sprintf(tmp, "%s_rt_mean", lname); - ConstTwoDMatrix m(nc.getDim(), 1, mean.base()); + ConstTwoDMatrix m(nc.getDim(), 1, mean); m.writeMat(fd, tmp); sprintf(tmp, "%s_rt_vcov", lname); ConstTwoDMatrix(vcov).writeMat(fd, tmp); @@ -500,7 +500,7 @@ SimulationIRFWorker::operator()() = new ExplicitShockRealization(res.control.getShocks(idata)); esr->addToShock(ishock, 0, imp); const TwoDMatrix &data = res.control.getData(idata); - ConstVector st(data, res.control.getNumBurn()); + Vector st{data.getCol(res.control.getNumBurn())}; TwoDMatrix *m = dr.simulate(em, np, st, *esr); m->add(-1.0, res.control.getData(idata)); { @@ -604,7 +604,7 @@ ExplicitShockRealization::ExplicitShockRealization(ShockRealization &sr, { for (int j = 0; j < num_per; j++) { - Vector jcol(shocks, j); + Vector jcol{shocks.getCol(j)}; sr.get(j, jcol); } } @@ -615,7 +615,7 @@ ExplicitShockRealization::get(int n, Vector &out) KORD_RAISE_IF(out.length() != numShocks(), "Wrong length of out vector in ExplicitShockRealization::get"); int i = n % shocks.ncols(); - ConstVector icol(shocks, i); + ConstVector icol{shocks.getCol(i)}; out = icol; } diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh index a7259dc08..5b9fa80da 100644 --- a/dynare++/kord/decision_rule.hh +++ b/dynare++/kord/decision_rule.hh @@ -57,7 +57,7 @@ public: enum emethod { horner, trad }; virtual ~DecisionRule() = default; - virtual TwoDMatrix *simulate(emethod em, int np, const Vector &ystart, + virtual TwoDMatrix *simulate(emethod em, int np, const ConstVector &ystart, ShockRealization &sr) const = 0; virtual void eval(emethod em, Vector &out, const ConstVector &v) const = 0; virtual void evaluate(emethod em, Vector &out, const ConstVector &ys, @@ -103,18 +103,18 @@ protected: const int nu; public: DecisionRuleImpl(const _Tparent &pol, const PartitionY &yp, int nuu, - const Vector &ys) + const ConstVector &ys) : ctraits::Tpol(pol), ysteady(ys), ypart(yp), nu(nuu) { } DecisionRuleImpl(_Tparent &pol, const PartitionY &yp, int nuu, - const Vector &ys) + const ConstVector &ys) : ctraits::Tpol(0, yp.ny(), pol), ysteady(ys), ypart(yp), nu(nuu) { } DecisionRuleImpl(const _Tg &g, const PartitionY &yp, int nuu, - const Vector &ys, double sigma) + const ConstVector &ys, double sigma) : ctraits::Tpol(yp.ny(), yp.nys()+nuu), ysteady(ys), ypart(yp), nu(nuu) { fillTensors(g, sigma); @@ -130,7 +130,7 @@ public: { return ysteady; } - TwoDMatrix *simulate(emethod em, int np, const Vector &ystart, + TwoDMatrix *simulate(emethod em, int np, const ConstVector &ystart, ShockRealization &sr) const override; void evaluate(emethod em, Vector &out, const ConstVector &ys, const ConstVector &u) const override; @@ -280,7 +280,7 @@ DecisionRuleImpl::centralize(const DecisionRuleImpl &dr) template TwoDMatrix * -DecisionRuleImpl::simulate(emethod em, int np, const Vector &ystart, +DecisionRuleImpl::simulate(emethod em, int np, const ConstVector &ystart, ShockRealization &sr) const { KORD_RAISE_IF(ysteady.length() != ystart.length(), @@ -303,7 +303,7 @@ DecisionRuleImpl::simulate(emethod em, int np, const Vector &ystart, dy = ystart_pred; dy.add(-1.0, ysteady_pred); sr.get(0, u); - Vector out(*res, 0); + Vector out{res->getCol(0)}; eval(em, out, dyu); // perform all other steps of simulations @@ -312,11 +312,11 @@ DecisionRuleImpl::simulate(emethod em, int np, const Vector &ystart, int i = 1; while (i < np) { - ConstVector ym(*res, i-1); + ConstVector ym{res->getCol(i-1)}; ConstVector dym(ym, ypart.nstat, ypart.nys()); dy = dym; sr.get(i, u); - Vector out(*res, i); + Vector out{res->getCol(i)}; eval(em, out, dyu); if (!out.isFinite()) { @@ -335,7 +335,7 @@ DecisionRuleImpl::simulate(emethod em, int np, const Vector &ystart, and leave the padded columns to zero. */ for (int j = 0; j < i; j++) { - Vector col(*res, j); + Vector col{res->getCol(j)}; col.add(1.0, ysteady); } @@ -415,17 +415,17 @@ class FoldDecisionRule : public DecisionRuleImpl friend class UnfoldDecisionRule; public: FoldDecisionRule(const ctraits::Tpol &pol, const PartitionY &yp, int nuu, - const Vector &ys) + const ConstVector &ys) : DecisionRuleImpl(pol, yp, nuu, ys) { } FoldDecisionRule(ctraits::Tpol &pol, const PartitionY &yp, int nuu, - const Vector &ys) + const ConstVector &ys) : DecisionRuleImpl(pol, yp, nuu, ys) { } FoldDecisionRule(const ctraits::Tg &g, const PartitionY &yp, int nuu, - const Vector &ys, double sigma) + const ConstVector &ys, double sigma) : DecisionRuleImpl(g, yp, nuu, ys, sigma) { } @@ -445,17 +445,17 @@ class UnfoldDecisionRule : public DecisionRuleImpl friend class FoldDecisionRule; public: UnfoldDecisionRule(const ctraits::Tpol &pol, const PartitionY &yp, int nuu, - const Vector &ys) + const ConstVector &ys) : DecisionRuleImpl(pol, yp, nuu, ys) { } UnfoldDecisionRule(ctraits::Tpol &pol, const PartitionY &yp, int nuu, - const Vector &ys) + const ConstVector &ys) : DecisionRuleImpl(pol, yp, nuu, ys) { } UnfoldDecisionRule(const ctraits::Tg &g, const PartitionY &yp, int nuu, - const Vector &ys, double sigma) + const ConstVector &ys, double sigma) : DecisionRuleImpl(g, yp, nuu, ys, sigma) { } diff --git a/dynare++/kord/dynamic_model.hh b/dynare++/kord/dynamic_model.hh index b50977c50..58021f01d 100644 --- a/dynare++/kord/dynamic_model.hh +++ b/dynare++/kord/dynamic_model.hh @@ -107,9 +107,9 @@ public: virtual Vector&getSteady() = 0; virtual void solveDeterministicSteady() = 0; - virtual void evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) = 0; - virtual void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, - const Vector &yyp, const Vector &xx) = 0; + virtual void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) = 0; + virtual void evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy, + const ConstVector &yyp, const Vector &xx) = 0; virtual void calcDerivativesAtSteady() = 0; }; diff --git a/dynare++/kord/global_check.cc b/dynare++/kord/global_check.cc index 6cfb304ba..94b957e82 100644 --- a/dynare++/kord/global_check.cc +++ b/dynare++/kord/global_check.cc @@ -52,7 +52,7 @@ ResidFunction::~ResidFunction() |hss|. */ void -ResidFunction::setYU(const Vector &ys, const Vector &xx) +ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx) { // delete |y| and |u| dependent data /* NB: code shared with the destructor */ @@ -195,9 +195,9 @@ GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y, ConstTwoDMatrix ysmat(y, first_row, 0, model.npred()+model.nboth(), y.ncols()); for (int j = 0; j < y.ncols(); j++) { - ConstVector yj(ysmat, j); - ConstVector xj(x, j); - Vector outj(out, j); + ConstVector yj{ysmat.getCol(j)}; + ConstVector xj{x.getCol(j)}; + Vector outj{out.getCol(j)}; check(*quad, lev, yj, xj, outj); } @@ -221,7 +221,7 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix, TwoDMatrix y_mat(model.numeq(), 2*m *model.nexog()+1); for (int j = 0; j < 2*m*model.nexog()+1; j++) { - Vector yj(y_mat, j); + Vector yj{y_mat.getCol(j)}; yj = (const Vector &) model.getSteady(); } @@ -246,7 +246,7 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix, TwoDMatrix res(model.nexog(), 2*m+1); JournalRecord rec(journal); rec << "Shock value error" << endrec; - ConstVector err0(errors, 0); + ConstVector err0{errors.getCol(0)}; char shock[9]; char erbuf[17]; for (int ishock = 0; ishock < model.nexog(); ishock++) @@ -256,14 +256,14 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix, for (int j = 0; j < 2*m+1; j++) { int jj; - Vector error(err_out, j); + Vector error{err_out.getCol(j)}; if (j != m) { if (j < m) jj = 1 + 2*m*ishock+j; else jj = 1 + 2*m*ishock+j-1; - ConstVector coljj(errors, jj); + ConstVector coljj{errors.getCol(jj)}; error = coljj; } else @@ -347,7 +347,7 @@ GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const char *prefix, qmcpit end = qmc.end(m); for (qmcpit run = beg; run != end; ++run, icol++) { - Vector ycol(ymat, icol); + Vector ycol{ymat.getCol(icol)}; Vector x(run.point()); x.mult(2*M_PI); ycol[0] = 1; @@ -372,7 +372,7 @@ GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const char *prefix, model.npred()+model.nboth()); for (int icol = 0; icol < ymat.ncols(); icol++) { - Vector ycol(ymat, icol); + Vector ycol{ymat.getCol(icol)}; ycol.add(1.0, ys); } diff --git a/dynare++/kord/global_check.hh b/dynare++/kord/global_check.hh index e49d23923..e35be844e 100644 --- a/dynare++/kord/global_check.hh +++ b/dynare++/kord/global_check.hh @@ -84,7 +84,7 @@ public: return std::make_unique(*this); } void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; - void setYU(const Vector &ys, const Vector &xx); + void setYU(const ConstVector &ys, const ConstVector &xx); }; /* This is a |ResidFunction| wrapped with |GaussConverterFunction|. */ @@ -104,7 +104,7 @@ public: return std::make_unique(*this); } void - setYU(const Vector &ys, const Vector &xx) + setYU(const ConstVector &ys, const ConstVector &xx) { ((ResidFunction *) func)->setYU(ys, xx); } diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc index 54ab42795..c9bcef5be 100644 --- a/dynare++/kord/korder.cc +++ b/dynare++/kord/korder.cc @@ -307,8 +307,8 @@ KOrder::sylvesterSolve(ctraits::Ttensor &der) const TwoDMatrix gs_y(*(gs().get(Symmetry(1, 0, 0, 0)))); GeneralSylvester sylv(der.getSym()[0], ny, ypart.nys(), ypart.nstat+ypart.npred, - matA.getData().base(), matB.getData().base(), - gs_y.getData().base(), der.getData().base()); + matA.getData(), matB.getData(), + gs_y.getData(), der.getData()); sylv.solve(); } else if (ypart.nys() > 0 && ypart.nyss() == 0) diff --git a/dynare++/kord/normal_conjugate.cc b/dynare++/kord/normal_conjugate.cc index 09534af36..a8e231b44 100644 --- a/dynare++/kord/normal_conjugate.cc +++ b/dynare++/kord/normal_conjugate.cc @@ -18,12 +18,12 @@ NormalConj::NormalConj(const ConstTwoDMatrix &ydata) { mu.zeros(); for (int i = 0; i < ydata.numCols(); i++) - mu.add(1.0/ydata.numCols(), ConstVector(ydata, i)); + mu.add(1.0/ydata.numCols(), ydata.getCol(i)); lambda.zeros(); for (int i = 0; i < ydata.numCols(); i++) { - Vector diff(ConstVector(ydata, i)); + Vector diff{ydata.getCol(i)}; diff.add(-1, mu); lambda.addOuter(diff); } diff --git a/dynare++/kord/tests.cc b/dynare++/kord/tests.cc index 68239d655..ee2d4633e 100644 --- a/dynare++/kord/tests.cc +++ b/dynare++/kord/tests.cc @@ -13,6 +13,13 @@ struct Rand static bool discrete(double prob); // answers true with given probability }; +TwoDMatrix +make_matrix(int rows, int cols, const double *p) +{ + return TwoDMatrix{rows, cols, + ConstVector{std::shared_ptr{p, [](const double *arr) {}}, rows*cols}}; +} + void Rand::init(int n1, int n2, int n3, int n4, int n5) { @@ -345,9 +352,9 @@ public: bool run() const override { - TwoDMatrix gy(8, 4, gy_data); - TwoDMatrix gu(8, 3, gu_data); - TwoDMatrix v(3, 3, vdata); + TwoDMatrix gy{make_matrix(8, 4, gy_data)}; + TwoDMatrix gu{make_matrix(8, 3, gu_data)}; + TwoDMatrix v{make_matrix(3, 3, vdata)}; double err = korder_unfold_fold(4, 3, 2, 3, 1, 2, gy, gu, v); @@ -368,9 +375,9 @@ public: bool run() const override { - TwoDMatrix gy(30, 20, gy_data2); - TwoDMatrix gu(30, 10, gu_data2); - TwoDMatrix v(10, 10, vdata2); + TwoDMatrix gy{make_matrix(30, 20, gy_data2)}; + TwoDMatrix gu{make_matrix(30, 10, gu_data2)}; + TwoDMatrix v{make_matrix(10, 10, vdata2)}; v.mult(0.001); gu.mult(.01); double err = korder_unfold_fold(4, 4, 5, 12, 8, 5, @@ -392,9 +399,9 @@ public: bool run() const override { - TwoDMatrix gy(30, 20, gy_data2); - TwoDMatrix gu(30, 10, gu_data2); - TwoDMatrix v(10, 10, vdata2); + TwoDMatrix gy{make_matrix(30, 20, gy_data2)}; + TwoDMatrix gu{make_matrix(30, 10, gu_data2)}; + TwoDMatrix v{make_matrix(10, 10, vdata2)}; v.mult(0.001); gu.mult(.01); double err = korder_unfold_fold(4, 3, 5, 12, 8, 5, diff --git a/dynare++/src/dynare3.cc b/dynare++/src/dynare3.cc index bb55d2c96..ebbf90654 100644 --- a/dynare++/src/dynare3.cc +++ b/dynare++/src/dynare3.cc @@ -199,7 +199,7 @@ Dynare::solveDeterministicSteady(Vector &steady) // evaluate system at given y_t=y_{t+1}=y_{t-1}, and given shocks x_t void -Dynare::evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) +Dynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) { ConstVector yym(yy, nstat(), nys()); ConstVector yyp(yy, nstat()+npred(), nyss()); @@ -210,8 +210,8 @@ Dynare::evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) // exogenous x_t, all three vectors yym, yy, and yyp have the // respective lengths of y^*_{t-1}, y_t, y^{**}_{t+1} void -Dynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, - const Vector &yyp, const Vector &xx) +Dynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy, + const ConstVector &yyp, const Vector &xx) { ogdyn::DynareAtomValues dav(model->getAtoms(), model->getParams(), yym, yy, yyp, xx); DynareEvalLoader del(model->getAtoms(), out); diff --git a/dynare++/src/dynare3.hh b/dynare++/src/dynare3.hh index 9faed9f25..fc6f0a073 100644 --- a/dynare++/src/dynare3.hh +++ b/dynare++/src/dynare3.hh @@ -228,9 +228,9 @@ public: { solveDeterministicSteady(*ysteady); } - void evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) override; - void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, - const Vector &yyp, const Vector &xx) override; + void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) override; + void evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy, + const ConstVector &yyp, const Vector &xx) override; void calcDerivatives(const Vector &yy, const Vector &xx); void calcDerivativesAtSteady() override; diff --git a/dynare++/src/dynare_atoms.hh b/dynare++/src/dynare_atoms.hh index 1a781453a..e6309a2f9 100644 --- a/dynare++/src/dynare_atoms.hh +++ b/dynare++/src/dynare_atoms.hh @@ -112,7 +112,7 @@ namespace ogdyn { } DynareAtomValues(const ogp::FineAtoms &a, const Vector &pvals, const ConstVector &ym, - const Vector &y, const ConstVector &yp, const Vector &x) + const ConstVector &y, const ConstVector &yp, const Vector &x) : atoms(a), paramvals(pvals), yym(ym), yy(y), yyp(yp), xx(x) { } diff --git a/dynare++/src/main.cc b/dynare++/src/main.cc index 7f5667968..39dba693c 100644 --- a/dynare++/src/main.cc +++ b/dynare++/src/main.cc @@ -142,7 +142,7 @@ main(int argc, char **argv) if (params.num_condper > 0 && params.num_condsim > 0) { SimResultsDynamicStats rescond(dynare.numeq(), params.num_condper, 0); - ConstVector det_ss(app.getSS(), 0); + Vector det_ss{app.getSS().getCol(0)}; rescond.simulate(params.num_condsim, app.getFoldDecisionRule(), det_ss, dynare.getVcov(), journal); rescond.writeMat(matfd, params.prefix); } diff --git a/dynare++/sylv/cc/BlockDiagonal.cc b/dynare++/sylv/cc/BlockDiagonal.cc index 88db35091..ecaca4b5a 100644 --- a/dynare++/sylv/cc/BlockDiagonal.cc +++ b/dynare++/sylv/cc/BlockDiagonal.cc @@ -6,9 +6,10 @@ #include #include +#include -BlockDiagonal::BlockDiagonal(const double *d, int d_size) - : QuasiTriangular(d, d_size), +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]) { for (int i = 0; i < d_size; i++) diff --git a/dynare++/sylv/cc/BlockDiagonal.hh b/dynare++/sylv/cc/BlockDiagonal.hh index 3b5c18546..a2ab6233a 100644 --- a/dynare++/sylv/cc/BlockDiagonal.hh +++ b/dynare++/sylv/cc/BlockDiagonal.hh @@ -14,7 +14,7 @@ class BlockDiagonal : public QuasiTriangular int *const row_len; int *const col_len; public: - BlockDiagonal(const double *d, int d_size); + BlockDiagonal(ConstVector d, int d_size); BlockDiagonal(int p, const BlockDiagonal &b); BlockDiagonal(const BlockDiagonal &b); BlockDiagonal(const QuasiTriangular &t); diff --git a/dynare++/sylv/cc/GeneralMatrix.cc b/dynare++/sylv/cc/GeneralMatrix.cc index 912a4657f..7a659c973 100644 --- a/dynare++/sylv/cc/GeneralMatrix.cc +++ b/dynare++/sylv/cc/GeneralMatrix.cc @@ -16,8 +16,6 @@ #include #include -int GeneralMatrix::md_length = 32; - GeneralMatrix::GeneralMatrix(const GeneralMatrix &m) : data(m.rows*m.cols), rows(m.rows), cols(m.cols), ld(m.rows) { @@ -53,35 +51,59 @@ GeneralMatrix::GeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, in } GeneralMatrix::GeneralMatrix(GeneralMatrix &m, int i, int j, int nrows, int ncols) - : data(m.base()+m.ld*j+i, m.ld*(ncols-1)+nrows), rows(nrows), cols(ncols), ld(m.ld) + : data(m.data, m.ld*j+i, m.ld*(ncols-1)+nrows), rows(nrows), cols(ncols), ld(m.ld) { } -GeneralMatrix::GeneralMatrix(const GeneralMatrix &a, const GeneralMatrix &b) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b) : data(a.rows*b.cols), rows(a.rows), cols(b.cols), ld(a.rows) { gemm("N", a, "N", b, 1.0, 0.0); } -GeneralMatrix::GeneralMatrix(const GeneralMatrix &a, const GeneralMatrix &b, const char *dum) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, const char *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 GeneralMatrix &a, const char *dum, const GeneralMatrix &b) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const char *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 GeneralMatrix &a, const char *dum1, - const GeneralMatrix &b, const char *dum2) +GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix &a, const char *dum1, + const ConstGeneralMatrix &b, const char *dum2) : data(a.cols*b.rows), rows(a.cols), cols(b.rows), ld(a.cols) { gemm("T", a, "T", b, 1.0, 0.0); } +Vector +GeneralMatrix::getRow(int row) +{ + return Vector{data, row, ld, cols}; +} + +Vector +GeneralMatrix::getCol(int col) +{ + return Vector{data, ld*col, rows}; +} + +ConstVector +GeneralMatrix::getRow(int row) const +{ + return ConstVector{data, row, ld, cols}; +} + +ConstVector +GeneralMatrix::getCol(int col) const +{ + return ConstVector{data, ld*col, rows}; +} + void GeneralMatrix::place(const ConstGeneralMatrix &m, int i, int j) { @@ -357,14 +379,25 @@ ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix &m) { } +ConstVector +ConstGeneralMatrix::getRow(int row) const +{ + return ConstVector{data, row, ld, cols}; +} + +ConstVector +ConstGeneralMatrix::getCol(int col) const +{ + return ConstVector{data, ld*col, rows}; +} + double ConstGeneralMatrix::getNormInf() const { double norm = 0.0; for (int i = 0; i < numRows(); i++) { - ConstVector rowi(data.base()+i, ld, cols); - double normi = rowi.getNorm1(); + double normi = getRow(i).getNorm1(); if (norm < normi) norm = normi; } @@ -377,8 +410,7 @@ ConstGeneralMatrix::getNorm1() const double norm = 0.0; for (int j = 0; j < numCols(); j++) { - ConstVector colj(data.base()+ld *j, 1, rows); - double normj = colj.getNorm1(); + double normj = getCol(j).getNorm1(); if (norm < normj) norm = normj; } @@ -581,12 +613,12 @@ SVDDecomp::solve(const GeneralMatrix &B, GeneralMatrix &X) const GeneralMatrix Bprime(UTB, m-minmn, 0, minmn-nz, B.numCols()); // solve sigma for (int i = 0; i < minmn-nz; i++) - Vector(i, Bprime).mult(1.0/sigma[i]); + Bprime.getRow(i).mult(1.0/sigma[i]); // solve VT X.zeros(); //- copy Bprime to right place of X for (int i = 0; i < minmn-nz; i++) - Vector(n-minmn+i, X) = ConstVector(i, Bprime); + X.getRow(n-minmn+i) = Bprime.getRow(i); //- multiply with VT X.multLeftTrans(VT); } diff --git a/dynare++/sylv/cc/GeneralMatrix.hh b/dynare++/sylv/cc/GeneralMatrix.hh index 968061e4d..e0eb4c59f 100644 --- a/dynare++/sylv/cc/GeneralMatrix.hh +++ b/dynare++/sylv/cc/GeneralMatrix.hh @@ -8,6 +8,8 @@ #include "Vector.hh" #include +#include +#include class GeneralMatrix; @@ -20,10 +22,11 @@ protected: int cols; int ld; public: - ConstGeneralMatrix(const double *d, int m, int n) - : data(d, m*n), rows(m), cols(n), ld(m) + ConstGeneralMatrix(ConstVector d, int m, int n) + : data(std::move(d)), rows(m), cols(n), ld(m) { } + // Implicit conversion from ConstGeneralMatrix is ok, since it's cheap ConstGeneralMatrix(const GeneralMatrix &m); ConstGeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols); ConstGeneralMatrix(const ConstGeneralMatrix &m, int i, int j, int nrows, int ncols); @@ -59,6 +62,8 @@ public: { return data; } + ConstVector getRow(int row) const; + ConstVector getCol(int col) const; double getNormInf() const; double getNorm1() const; @@ -121,29 +126,30 @@ public: : data(m*n), rows(m), cols(n), ld(m) { } - GeneralMatrix(const double *d, int m, int n) - : data(d, m*n), rows(m), cols(n), ld(m) + GeneralMatrix(const ConstVector &d, int m, int n) + : data(d), rows(m), cols(n), ld(m) { } - GeneralMatrix(double *d, int m, int n) - : data(d, m*n), rows(m), cols(n), ld(m) + GeneralMatrix(Vector &d, int m, int n) + : data(d), rows(m), cols(n), ld(m) { } GeneralMatrix(const GeneralMatrix &m); - GeneralMatrix(const ConstGeneralMatrix &m); + // We don't want implict conversion from ConstGeneralMatrix, since it's expensive + explicit GeneralMatrix(const ConstGeneralMatrix &m); GeneralMatrix(const GeneralMatrix &m, const char *dummy); // transpose GeneralMatrix(const ConstGeneralMatrix &m, const char *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 GeneralMatrix &a, const GeneralMatrix &b); + GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b); /* this = a*b' */ - GeneralMatrix(const GeneralMatrix &a, const GeneralMatrix &b, const char *dum); + GeneralMatrix(const ConstGeneralMatrix &a, const ConstGeneralMatrix &b, const char *dum); /* this = a'*b */ - GeneralMatrix(const GeneralMatrix &a, const char *dum, const GeneralMatrix &b); + GeneralMatrix(const ConstGeneralMatrix &a, const char *dum, const ConstGeneralMatrix &b); /* this = a'*b */ - GeneralMatrix(const GeneralMatrix &a, const char *dum1, - const GeneralMatrix &b, const char *dum2); + GeneralMatrix(const ConstGeneralMatrix &a, const char *dum1, + const ConstGeneralMatrix &b, const char *dum2); virtual ~GeneralMatrix() = default; GeneralMatrix &operator=(const GeneralMatrix &m) = default; @@ -193,6 +199,10 @@ public: { return data; } + Vector getRow(int row); + Vector getCol(int col); + ConstVector getRow(int row) const; + ConstVector getCol(int col) const; double getNormInf() const @@ -432,7 +442,7 @@ private: } /* number of rows/columns for copy used in gemm_partial_* */ - static int md_length; + static constexpr int md_length = 23; }; class SVDDecomp @@ -473,8 +483,8 @@ public: void solve(const Vector &b, Vector &x) const { - GeneralMatrix xmat(x.base(), x.length(), 1); - solve(GeneralMatrix(b.base(), b.length(), 1), xmat); + GeneralMatrix xmat(x, x.length(), 1); + solve(GeneralMatrix(b, b.length(), 1), xmat); } private: void construct(const GeneralMatrix &A); diff --git a/dynare++/sylv/cc/GeneralSylvester.cc b/dynare++/sylv/cc/GeneralSylvester.cc index bc5c031e1..52ae68f1c 100644 --- a/dynare++/sylv/cc/GeneralSylvester.cc +++ b/dynare++/sylv/cc/GeneralSylvester.cc @@ -11,8 +11,8 @@ #include GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, const double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, const ConstVector &dd, const SylvParams &ps) : pars(ps), mem_driver(pars, 1, m, n, ord), order(ord), a(da, n), @@ -23,8 +23,8 @@ GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols, } GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, Vector &dd, const SylvParams &ps) : pars(ps), mem_driver(pars, 0, m, n, ord), order(ord), a(da, n), @@ -35,8 +35,8 @@ GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols, } GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, const double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, const ConstVector &dd, bool alloc_for_check) : pars(alloc_for_check), mem_driver(pars, 1, m, n, ord), order(ord), a(da, n), @@ -47,8 +47,8 @@ GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols, } GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, Vector &dd, bool alloc_for_check) : pars(alloc_for_check), mem_driver(pars, 0, m, n, ord), order(ord), a(da, n), @@ -68,7 +68,7 @@ GeneralSylvester::init() pars.rcondA1 = rcond1; pars.rcondAI = rcondinf; bdecomp = std::make_unique(ainvb); - cdecomp = std::make_unique(c.getData().base(), c.numRows(), *(pars.bs_norm)); + cdecomp = std::make_unique(c.getData(), c.numRows(), *(pars.bs_norm)); cdecomp->check(pars, c); cdecomp->infoToPars(pars); if (*(pars.method) == SylvParams::recurse) @@ -105,7 +105,7 @@ GeneralSylvester::solve() } void -GeneralSylvester::check(const double *ds) +GeneralSylvester::check(const ConstVector &ds) { if (!solved) throw SYLV_MES_EXCEPTION("Cannot run check on system, which is not solved yet."); @@ -117,8 +117,7 @@ GeneralSylvester::check(const double *ds) dcheck.multLeft(b.numRows()-b.numCols(), b, d); dcheck.multRightKron(c, order); dcheck.multAndAdd(a, d); - ConstVector dv(ds, d.numRows()*d.numCols()); - dcheck.getData().add(-1.0, dv); + dcheck.getData().add(-1.0, ds); // calculate relative norms pars.mat_err1 = dcheck.getNorm1()/d.getNorm1(); pars.mat_errI = dcheck.getNormInf()/d.getNormInf(); diff --git a/dynare++/sylv/cc/GeneralSylvester.hh b/dynare++/sylv/cc/GeneralSylvester.hh index 489169ae5..b1cd69a30 100644 --- a/dynare++/sylv/cc/GeneralSylvester.hh +++ b/dynare++/sylv/cc/GeneralSylvester.hh @@ -28,21 +28,21 @@ class GeneralSylvester public: /* construct with my copy of d*/ GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, const double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, const ConstVector &dd, const SylvParams &ps); GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, const double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, const ConstVector &dd, bool alloc_for_check = false); /* construct with provided storage for d */ GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, Vector &dd, bool alloc_for_check = false); GeneralSylvester(int ord, int n, int m, int zero_cols, - const double *da, const double *db, - const double *dc, double *dd, + const ConstVector &da, const ConstVector &db, + const ConstVector &dc, Vector &dd, const SylvParams &ps); virtual ~GeneralSylvester() = default; int @@ -71,7 +71,7 @@ public: return pars; } void solve(); - void check(const double *ds); + void check(const ConstVector &ds); private: void init(); }; diff --git a/dynare++/sylv/cc/KronUtils.cc b/dynare++/sylv/cc/KronUtils.cc index b3b9740bc..46ab2dd55 100644 --- a/dynare++/sylv/cc/KronUtils.cc +++ b/dynare++/sylv/cc/KronUtils.cc @@ -18,7 +18,7 @@ KronUtils::multAtLevel(int level, const QuasiTriangular &t, } else if (0 == level && 0 < x.getDepth()) { - GeneralMatrix tmp(x.base(), x.getN(), power(x.getM(), x.getDepth())); + GeneralMatrix tmp(x, x.getN(), power(x.getM(), x.getDepth())); t.multLeftOther(tmp); } else if (0 == level && 0 == x.getDepth()) @@ -46,7 +46,7 @@ KronUtils::multAtLevelTrans(int level, const QuasiTriangular &t, } else if (0 == level && 0 < x.getDepth()) { - GeneralMatrix tmp(x.base(), x.getN(), power(x.getM(), x.getDepth())); + GeneralMatrix tmp(x, x.getN(), power(x.getM(), x.getDepth())); t.multLeftOtherTrans(tmp); } else if (level == 0 && 0 == x.getDepth()) diff --git a/dynare++/sylv/cc/KronVector.cc b/dynare++/sylv/cc/KronVector.cc index db27fa13c..8fc940f2c 100644 --- a/dynare++/sylv/cc/KronVector.cc +++ b/dynare++/sylv/cc/KronVector.cc @@ -5,6 +5,8 @@ #include "KronVector.hh" #include "SylvException.hh" +#include + int power(int m, int depth) { @@ -77,8 +79,8 @@ ConstKronVector::ConstKronVector(const Vector &v, int mm, int nn, int dp) throw SYLV_MES_EXCEPTION("Bad conversion KronVector from Vector."); } -ConstKronVector::ConstKronVector(const ConstVector &v, int mm, int nn, int dp) - : ConstVector(v), m(mm), n(nn), depth(dp) +ConstKronVector::ConstKronVector(ConstVector v, int mm, int nn, int dp) + : ConstVector(std::move(v)), m(mm), n(nn), depth(dp) { if (length() != power(m, depth)*n) throw SYLV_MES_EXCEPTION("Bad conversion KronVector from Vector."); diff --git a/dynare++/sylv/cc/KronVector.hh b/dynare++/sylv/cc/KronVector.hh index f58d0f67d..88e005d44 100644 --- a/dynare++/sylv/cc/KronVector.hh +++ b/dynare++/sylv/cc/KronVector.hh @@ -16,13 +16,12 @@ protected: int n{0}; int depth{0}; public: - KronVector() : Vector((double *) nullptr, 0) - { - } + KronVector() = 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 - KronVector(const ConstKronVector &v); // new instance and copy + // 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=(const ConstKronVector &v); KronVector &operator=(const Vector &v); @@ -50,10 +49,11 @@ protected: int n; int depth; public: + // Implicit conversion from KronVector is ok, since it's cheap ConstKronVector(const KronVector &v); ConstKronVector(const ConstKronVector &v); ConstKronVector(const Vector &v, int mm, int nn, int dp); - ConstKronVector(const ConstVector &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); int diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc index 618c81db5..b0fda9f97 100644 --- a/dynare++/sylv/cc/QuasiTriangular.cc +++ b/dynare++/sylv/cc/QuasiTriangular.cc @@ -420,7 +420,7 @@ QuasiTriangular::QuasiTriangular(const QuasiTriangular &t) { } -QuasiTriangular::QuasiTriangular(const double *d, int d_size) +QuasiTriangular::QuasiTriangular(const ConstVector &d, int d_size) : SqSylvMatrix(d, d_size), diagonal(getData().base(), d_size) { } @@ -679,8 +679,8 @@ void QuasiTriangular::multaKron(KronVector &x, const ConstKronVector &b) const { int id = b.getN()*power(b.getM(), b.getDepth()-1); - ConstGeneralMatrix b_resh(b.base(), id, b.getM()); - GeneralMatrix x_resh(x.base(), id, b.getM()); + ConstGeneralMatrix b_resh(b, id, b.getM()); + GeneralMatrix x_resh(x, id, b.getM()); x_resh.multAndAdd(b_resh, ConstGeneralMatrix(*this), "trans"); } @@ -689,8 +689,8 @@ void QuasiTriangular::multaKronTrans(KronVector &x, const ConstKronVector &b) const { int id = b.getN()*power(b.getM(), b.getDepth()-1); - ConstGeneralMatrix b_resh(b.base(), id, b.getM()); - GeneralMatrix x_resh(x.base(), id, b.getM()); + ConstGeneralMatrix b_resh(b, id, b.getM()); + GeneralMatrix x_resh(x, id, b.getM()); x_resh.multAndAdd(b_resh, ConstGeneralMatrix(*this)); } diff --git a/dynare++/sylv/cc/QuasiTriangular.hh b/dynare++/sylv/cc/QuasiTriangular.hh index 773508020..f8d3c8e7e 100644 --- a/dynare++/sylv/cc/QuasiTriangular.hh +++ b/dynare++/sylv/cc/QuasiTriangular.hh @@ -399,7 +399,7 @@ public: protected: Diagonal diagonal; public: - QuasiTriangular(const double *d, int d_size); + QuasiTriangular(const ConstVector &d, int d_size); QuasiTriangular(double r, const QuasiTriangular &t); QuasiTriangular(double r, const QuasiTriangular &t, double rr, const QuasiTriangular &tt); diff --git a/dynare++/sylv/cc/QuasiTriangularZero.cc b/dynare++/sylv/cc/QuasiTriangularZero.cc index 665b6cd90..fca2614b0 100644 --- a/dynare++/sylv/cc/QuasiTriangularZero.cc +++ b/dynare++/sylv/cc/QuasiTriangularZero.cc @@ -8,14 +8,15 @@ #include "SylvException.hh" #include +#include -QuasiTriangularZero::QuasiTriangularZero(int num_zeros, const double *d, +QuasiTriangularZero::QuasiTriangularZero(int num_zeros, const ConstVector &d, int d_size) : QuasiTriangular(SqSylvMatrix(GeneralMatrix(d, num_zeros+d_size, d_size), - num_zeros, 0, d_size).getData().base(), + num_zeros, 0, d_size).getData(), d_size), nz(num_zeros), - ru(GeneralMatrix(d, num_zeros+d_size, d_size), 0, 0, num_zeros, d_size) + ru(GeneralMatrix(std::move(d), num_zeros+d_size, d_size), 0, 0, num_zeros, d_size) { } @@ -49,7 +50,7 @@ QuasiTriangularZero::QuasiTriangularZero(int p, const QuasiTriangularZero &t) } QuasiTriangularZero::QuasiTriangularZero(const SchurDecompZero &decomp) - : QuasiTriangular(decomp.getT().getData().base(), + : QuasiTriangular(decomp.getT().getData(), decomp.getT().numRows()), nz(decomp.getZeroCols()), ru(decomp.getRU()) diff --git a/dynare++/sylv/cc/QuasiTriangularZero.hh b/dynare++/sylv/cc/QuasiTriangularZero.hh index 154bc9ff3..853e04a1a 100644 --- a/dynare++/sylv/cc/QuasiTriangularZero.hh +++ b/dynare++/sylv/cc/QuasiTriangularZero.hh @@ -15,7 +15,7 @@ class QuasiTriangularZero : public QuasiTriangular int nz; // number of zero columns GeneralMatrix ru; // data in right upper part (nz,d_size) public: - QuasiTriangularZero(int num_zeros, const double *d, int d_size); + QuasiTriangularZero(int num_zeros, const ConstVector &d, int d_size); QuasiTriangularZero(double r, const QuasiTriangularZero &t); QuasiTriangularZero(double r, const QuasiTriangularZero &t, double rr, const QuasiTriangularZero &tt); diff --git a/dynare++/sylv/cc/SchurDecomp.cc b/dynare++/sylv/cc/SchurDecomp.cc index f4f85ba3d..7a101d47d 100644 --- a/dynare++/sylv/cc/SchurDecomp.cc +++ b/dynare++/sylv/cc/SchurDecomp.cc @@ -21,7 +21,7 @@ SchurDecomp::SchurDecomp(const SqSylvMatrix &m) dgees("V", "N", nullptr, &rows, auxt.base(), &rows, &sdim, wr.data(), wi.data(), q.base(), &rows, work.data(), &lwork, nullptr, &info); - t_storage = std::make_unique(auxt.base(), rows); + t_storage = std::make_unique(auxt.getData(), rows); t = t_storage.get(); } diff --git a/dynare++/sylv/cc/SimilarityDecomp.cc b/dynare++/sylv/cc/SimilarityDecomp.cc index 3de344500..0da2ebf7e 100644 --- a/dynare++/sylv/cc/SimilarityDecomp.cc +++ b/dynare++/sylv/cc/SimilarityDecomp.cc @@ -11,7 +11,7 @@ #include -SimilarityDecomp::SimilarityDecomp(const double *d, int d_size, double log10norm) +SimilarityDecomp::SimilarityDecomp(const ConstVector &d, int d_size, double log10norm) { SchurDecomp sd(SqSylvMatrix(d, d_size)); q = std::make_unique(sd.getQ()); diff --git a/dynare++/sylv/cc/SimilarityDecomp.hh b/dynare++/sylv/cc/SimilarityDecomp.hh index b831c821d..7dc7074c6 100644 --- a/dynare++/sylv/cc/SimilarityDecomp.hh +++ b/dynare++/sylv/cc/SimilarityDecomp.hh @@ -18,7 +18,7 @@ class SimilarityDecomp std::unique_ptr invq; using diag_iter = BlockDiagonal::diag_iter; public: - SimilarityDecomp(const double *d, int d_size, double log10norm = 3.0); + SimilarityDecomp(const ConstVector &d, int d_size, double log10norm = 3.0); virtual ~SimilarityDecomp() = default; const SqSylvMatrix & getQ() const diff --git a/dynare++/sylv/cc/SylvMatrix.cc b/dynare++/sylv/cc/SylvMatrix.cc index 7744e895c..6af0480c0 100644 --- a/dynare++/sylv/cc/SylvMatrix.cc +++ b/dynare++/sylv/cc/SylvMatrix.cc @@ -69,7 +69,7 @@ SylvMatrix::multRightKron(const SqSylvMatrix &m, int order) KronVector auxrow(m.numRows(), m.numRows(), order-1); for (int i = 0; i < rows; i++) { - Vector rowi(data.base()+i, rows, cols); + Vector rowi{getRow(i)}; KronVector rowikron(rowi, m.numRows(), m.numRows(), order-1); auxrow = rowi; // copy data m.multVecKronTrans(rowikron, auxrow); @@ -85,7 +85,7 @@ SylvMatrix::multRightKronTrans(const SqSylvMatrix &m, int order) KronVector auxrow(m.numRows(), m.numRows(), order-1); for (int i = 0; i < rows; i++) { - Vector rowi(data.base()+i, rows, cols); + Vector rowi{getRow(i)}; KronVector rowikron(rowi, m.numRows(), m.numRows(), order-1); auxrow = rowi; // copy data m.multVecKron(rowikron, auxrow); @@ -169,7 +169,7 @@ SqSylvMatrix::SqSylvMatrix(const GeneralMatrix &a, const GeneralMatrix &b) } void -SqSylvMatrix::multVecKron(KronVector &x, const KronVector &d) const +SqSylvMatrix::multVecKron(KronVector &x, const ConstKronVector &d) const { x.zeros(); if (d.getDepth() == 0) @@ -196,7 +196,7 @@ SqSylvMatrix::multVecKron(KronVector &x, const KronVector &d) const } void -SqSylvMatrix::multVecKronTrans(KronVector &x, const KronVector &d) const +SqSylvMatrix::multVecKronTrans(KronVector &x, const ConstKronVector &d) const { x.zeros(); if (d.getDepth() == 0) diff --git a/dynare++/sylv/cc/SylvMatrix.hh b/dynare++/sylv/cc/SylvMatrix.hh index 54a59509b..b23ddd9e9 100644 --- a/dynare++/sylv/cc/SylvMatrix.hh +++ b/dynare++/sylv/cc/SylvMatrix.hh @@ -17,11 +17,11 @@ public: : GeneralMatrix(m, n) { } - SylvMatrix(const double *d, int m, int n) + SylvMatrix(const ConstVector &d, int m, int n) : GeneralMatrix(d, m, n) { } - SylvMatrix(double *d, int m, int n) + SylvMatrix(Vector &d, int m, int n) : GeneralMatrix(d, m, n) { } @@ -68,10 +68,10 @@ public: SqSylvMatrix(int m) : SylvMatrix(m, m) { } - SqSylvMatrix(const double *d, int m) : SylvMatrix(d, m, m) + SqSylvMatrix(const ConstVector &d, int m) : SylvMatrix(d, m, m) { } - SqSylvMatrix(double *d, int m) : SylvMatrix(d, m, m) + SqSylvMatrix(Vector &d, int m) : SylvMatrix(d, m, m) { } SqSylvMatrix(const SqSylvMatrix &m) = default; @@ -91,9 +91,9 @@ public: return *this; } /* x = (this \otimes this..\otimes this)*d */ - void multVecKron(KronVector &x, const KronVector &d) const; + void multVecKron(KronVector &x, const ConstKronVector &d) const; /* x = (this' \otimes this'..\otimes this')*d */ - void multVecKronTrans(KronVector &x, const KronVector &d) const; + void multVecKronTrans(KronVector &x, const ConstKronVector &d) const; /* a = inv(this)*a, b=inv(this)*b */ void multInvLeft2(GeneralMatrix &a, GeneralMatrix &b, double &rcond1, double &rcondinf) const; diff --git a/dynare++/sylv/cc/SymSchurDecomp.cc b/dynare++/sylv/cc/SymSchurDecomp.cc index b11ee363c..0c7466e4c 100644 --- a/dynare++/sylv/cc/SymSchurDecomp.cc +++ b/dynare++/sylv/cc/SymSchurDecomp.cc @@ -76,7 +76,7 @@ SymSchurDecomp::getFactor(GeneralMatrix &f) const f = q; for (int i = 0; i < f.numCols(); i++) { - Vector fi(f, i); + Vector fi{f.getCol(i)}; fi.mult(std::sqrt(lambda[i])); } } diff --git a/dynare++/sylv/cc/TriangularSylvester.cc b/dynare++/sylv/cc/TriangularSylvester.cc index 8ba2e8bc2..3092b9d1f 100644 --- a/dynare++/sylv/cc/TriangularSylvester.cc +++ b/dynare++/sylv/cc/TriangularSylvester.cc @@ -390,7 +390,7 @@ TriangularSylvester::multEigVector(KronVector &eig, const Vector &feig, { KronVector eigi(eig, i); eigi.zeros(); - eigi.add(&feig[2*i], aux); + eigi.addComplex({ feig[2*i], feig[2*i+1] }, aux); } } } diff --git a/dynare++/sylv/cc/Vector.cc b/dynare++/sylv/cc/Vector.cc index c58bc089d..c902c52f6 100644 --- a/dynare++/sylv/cc/Vector.cc +++ b/dynare++/sylv/cc/Vector.cc @@ -19,13 +19,13 @@ using namespace std; ZeroPad zero_pad; Vector::Vector(const Vector &v) - : len(v.length()), data(new double[len]), destroy(true) + : len(v.len), data{new double[len], [](double *arr) { delete[] arr; }} { - copy(v.base(), v.skip()); + copy(v.base(), v.s); } Vector::Vector(const ConstVector &v) - : len(v.length()), data(new double[len]), destroy(true) + : len(v.length()), data{new double[len], [](double *arr) { delete[] arr; }} { copy(v.base(), v.skip()); } @@ -36,9 +36,9 @@ Vector::operator=(const Vector &v) if (this == &v) return *this; - if (v.length() != length()) + 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) @@ -48,23 +48,32 @@ Vector::operator=(const Vector &v) << ", data-v.data=" << (unsigned long) (data-v.data) << ", len=" << len << std::endl; throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors."); - } - copy(v.base(), v.skip()); + } */ + copy(v.base(), v.s); + return *this; +} + +Vector & +Vector::operator=(Vector &&v) +{ + if (v.len != len) + throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths."); + copy(v.base(), v.s); return *this; } Vector & Vector::operator=(const ConstVector &v) { - if (v.length() != length()) + if (v.length() != 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()))) throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors."); - + */ copy(v.base(), v.skip()); return *this; } @@ -72,37 +81,48 @@ Vector::operator=(const ConstVector &v) void Vector::copy(const double *d, int inc) { - blas_int n = length(); - blas_int incy = skip(); + blas_int n = len; + blas_int incy = s; blas_int inc2 = inc; dcopy(&n, d, &inc2, base(), &incy); } -Vector::Vector(Vector &v, int off, int l) - : len(l), s(v.skip()), data(v.base()+off*v.skip()) +Vector::Vector(Vector &v, int off_arg, int l) + : len(l), off{v.off+off_arg*v.s}, s(v.s), data{v.data} { - if (off < 0 || off + length() > v.length()) + if (off_arg < 0 || off_arg + len > v.len) throw SYLV_MES_EXCEPTION("Subvector not contained in supvector."); } -Vector::Vector(const Vector &v, int off, int l) - : len(l), data(new double[len]), destroy(true) +Vector::Vector(const Vector &v, int off_arg, int l) + : len(l), data{new double[len], [](double *arr) { delete[] arr; }} { - if (off < 0 || off + length() > v.length()) + if (off_arg < 0 || off_arg + len > v.len) throw SYLV_MES_EXCEPTION("Subvector not contained in supvector."); - copy(v.base()+off*v.skip(), v.skip()); + copy(v.base()+off_arg*v.s, v.s); } -Vector::Vector(GeneralMatrix &m, int col) - : len(m.numRows()), data(&(m.get(0, col))) +Vector::Vector(Vector &v, int off_arg, int skip, int l) + : len(l), off{v.off+off_arg*v.s}, s(v.s*skip), data{v.data} { } -Vector::Vector(int row, GeneralMatrix &m) - : len(m.numCols()), s(m.getLD()), data(&(m.get(row, 0))) +Vector::Vector(const Vector &v, int off_arg, int skip, int l) + : len(l), data{new double[len], [](double *arr) { delete[] arr; }} { + copy(v.base()+off_arg*v.s, v.s*skip); } +#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE) +Vector::Vector(mxArray *p) + : len{static_cast(mxGetNumberOfElements(p))}, + data{std::shared_ptr{mxGetPr(p), [](double *arr) { }}} +{ + if (!mxIsDouble(p)) + throw SYLV_MES_EXCEPTION("This is not a MATLAB array of doubles."); +} +#endif + bool Vector::operator==(const Vector &y) const { @@ -142,18 +162,18 @@ Vector::operator>=(const Vector &y) const void Vector::zeros() { - if (skip() == 1) + if (s == 1) { double *p = base(); - for (int i = 0; i < length()/ZeroPad::length; + for (int i = 0; i < len/ZeroPad::length; i++, p += ZeroPad::length) memcpy(p, zero_pad.getBase(), sizeof(double)*ZeroPad::length); - for (; p < base()+length(); p++) + for (; p < base()+len; p++) *p = 0.0; } else { - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) operator[](i) = 0.0; } } @@ -161,23 +181,17 @@ Vector::zeros() void Vector::nans() { - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) operator[](i) = std::numeric_limits::quiet_NaN(); } void Vector::infs() { - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) operator[](i) = std::numeric_limits::infinity(); } -Vector::~Vector() -{ - if (destroy) - delete[] data; -} - void Vector::rotatePair(double alpha, double beta1, double beta2, int i) { @@ -195,32 +209,32 @@ Vector::add(double r, const Vector &v) void Vector::add(double r, const ConstVector &v) { - blas_int n = length(); + blas_int n = len; blas_int incx = v.skip(); - blas_int incy = skip(); + blas_int incy = s; daxpy(&n, &r, v.base(), &incx, base(), &incy); } void -Vector::add(const double *z, const Vector &v) +Vector::addComplex(const std::complex &z, const Vector &v) { - add(z, ConstVector(v)); + addComplex(z, ConstVector(v)); } void -Vector::add(const double *z, const ConstVector &v) +Vector::addComplex(const std::complex &z, const ConstVector &v) { - blas_int n = length()/2; + blas_int n = len/2; blas_int incx = v.skip(); - blas_int incy = skip(); - zaxpy(&n, z, v.base(), &incx, base(), &incy); + blas_int incy = s; + zaxpy(&n, reinterpret_cast(z), v.base(), &incx, base(), &incy); } void Vector::mult(double r) { - blas_int n = length(); - blas_int incx = skip(); + blas_int n = len; + blas_int incx = s; dscal(&n, &r, base(), &incx); } @@ -283,71 +297,74 @@ Vector::print() const { auto ff = std::cout.flags(); std::cout << std::setprecision(4); - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) std::cout << i << '\t' << std::setw(8) << operator[](i) << std::endl; std::cout.flags(ff); } -ConstVector::ConstVector(const Vector &v, int off, int l) - : BaseConstVector(l, v.skip(), v.base() + v.skip()*off) +ConstVector::ConstVector(const Vector &v) + : len{v.len}, off{v.off}, s{v.s}, data{v.data} { - if (off < 0 || off + length() > v.length()) +} + +ConstVector::ConstVector(const ConstVector &v, int off_arg, int l) + : len{l}, off{v.off+off_arg*v.s}, s{v.s}, data{v.data} +{ + if (off_arg < 0 || off_arg + len > v.len) throw SYLV_MES_EXCEPTION("Subvector not contained in supvector."); } -ConstVector::ConstVector(const ConstVector &v, int off, int l) - : BaseConstVector(l, v.skip(), v.base() + v.skip()*off) -{ - if (off < 0 || off + length() > v.length()) - throw SYLV_MES_EXCEPTION("Subvector not contained in supvector."); -} - -ConstVector::ConstVector(const double *d, int skip, int l) - : BaseConstVector(l, skip, d) +ConstVector::ConstVector(const ConstVector &v, int off_arg, int skip, int l) + : len(l), off{v.off+off_arg*v.s}, s{v.s*skip}, data{v.data} { } -ConstVector::ConstVector(const ConstGeneralMatrix &m, int col) - : BaseConstVector(m.numRows(), 1, &(m.get(0, col))) +ConstVector::ConstVector(std::shared_ptr d, int skip, int l) + : len{l}, s{skip}, data{std::move(d)} { } -ConstVector::ConstVector(int row, const ConstGeneralMatrix &m) - : BaseConstVector(m.numCols(), m.getLD(), &(m.get(row, 0))) +#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE) +ConstVector::ConstVector(const mxArray *p) + : len{static_cast(mxGetNumberOfElements(p))}, + data{std::shared_ptr{mxGetPr(p), [](const double *arr) { }}} { + if (!mxIsDouble(p)) + throw SYLV_MES_EXCEPTION("This is not a MATLAB array of doubles."); } +#endif bool ConstVector::operator==(const ConstVector &y) const { - if (length() != y.length()) + if (len != y.len) return false; - if (length() == 0) + if (len == 0) return true; int i = 0; - while (i < length() && operator[](i) == y[i]) + while (i < len && operator[](i) == y[i]) i++; - return i == length(); + return i == len; } bool ConstVector::operator<(const ConstVector &y) const { - int i = std::min(length(), y.length()); + int i = std::min(len, y.len); int ii = 0; while (ii < i && operator[](ii) == y[ii]) ii++; if (ii < i) return operator[](ii) < y[ii]; else - return length() < y.length(); + return len < y.len; } double ConstVector::getNorm() const { double s = 0; - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) s += operator[](i)*operator[](i); return sqrt(s); } @@ -356,7 +373,7 @@ double ConstVector::getMax() const { double r = 0; - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) if (abs(operator[](i)) > r) r = abs(operator[](i)); return r; @@ -366,7 +383,7 @@ double ConstVector::getNorm1() const { double norm = 0.0; - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) norm += abs(operator[](i)); return norm; } @@ -374,11 +391,11 @@ ConstVector::getNorm1() const double ConstVector::dot(const ConstVector &y) const { - if (length() != y.length()) + if (len != y.len) throw SYLV_MES_EXCEPTION("Vector has different length in ConstVector::dot."); - blas_int n = length(); - blas_int incx = skip(); - blas_int incy = y.skip(); + blas_int n = len; + blas_int incx = s; + blas_int incy = y.s; return ddot(&n, base(), &incx, y.base(), &incy); } @@ -386,9 +403,9 @@ bool ConstVector::isFinite() const { int i = 0; - while (i < length() && isfinite(operator[](i))) + while (i < len && isfinite(operator[](i))) i++; - return i == length(); + return i == len; } void @@ -396,7 +413,7 @@ ConstVector::print() const { auto ff = std::cout.flags(); std::cout << std::setprecision(4); - for (int i = 0; i < length(); i++) + for (int i = 0; i < len; i++) std::cout << i << '\t' << std::setw(8) << operator[](i) << std::endl; std::cout.flags(ff); } diff --git a/dynare++/sylv/cc/Vector.hh b/dynare++/sylv/cc/Vector.hh index 34e2ef422..025dc4c56 100644 --- a/dynare++/sylv/cc/Vector.hh +++ b/dynare++/sylv/cc/Vector.hh @@ -10,65 +10,67 @@ * members, and methods are thus duplicated */ #include +#include +#include + +#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE) +# include +#endif class GeneralMatrix; class ConstVector; class Vector { + friend class ConstVector; protected: int len{0}; - int s{1}; - double *data{nullptr}; - bool destroy{false}; + int off{0}; // offset to double* pointer + int s{1}; // stride (also called "skip" in some places) + std::shared_ptr data; public: Vector() = default; - Vector(int l) : len(l), data(new double[l]), destroy(true) - { - } - Vector(Vector &v) : len(v.length()), s(v.skip()), data(v.base()) + Vector(int l) : len{l}, data{new double[l], [](double *arr) { delete[] arr; }} { } + Vector(Vector &v) = default; Vector(const Vector &v); - Vector(const ConstVector &v); - Vector(const double *d, int l) - : len(l), data(new double[len]), destroy(true) - { - copy(d, 1); - } - Vector(double *d, int l) - : len(l), data(d) + Vector(Vector &&v) = default; + // We don't want implict conversion from ConstVector, since it's expensive + explicit Vector(const ConstVector &v); + Vector(std::shared_ptr d, int l) + : len(l), data{std::move(d)} { } - Vector(double *d, int skip, int l) - : len(l), s(skip), data(d) - { - } - Vector(Vector &v, int off, int l); - Vector(const Vector &v, int off, int l); - Vector(GeneralMatrix &m, int col); - Vector(int row, GeneralMatrix &m); + Vector(Vector &v, int off_arg, int l); + Vector(const Vector &v, int off_arg, int l); + Vector(Vector &v, int off_arg, int skip, int l); + Vector(const Vector &v, int off_arg, int skip, int l); +#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE) + explicit Vector(mxArray *p); +#endif Vector &operator=(const Vector &v); + Vector &operator=(Vector &&v); Vector &operator=(const ConstVector &v); double & operator[](int i) { - return data[s*i]; + return data.get()[off+s*i]; } const double & operator[](int i) const { - return data[s*i]; + return data.get()[off+s*i]; } const double * base() const { - return data; + return data.get() + off; } double * base() { - return data; + return data.get() + off; } int length() const @@ -90,20 +92,19 @@ public: bool operator>(const Vector &y) const; bool operator>=(const Vector &y) const; - virtual ~Vector(); + virtual ~Vector() = default; void zeros(); void nans(); void infs(); - bool - toBeDestroyed() const - { - return destroy; - } void rotatePair(double alpha, double beta1, double beta2, int i); + // Computes this += r*v void add(double r, const Vector &v); + // Computes this += r*v void add(double r, const ConstVector &v); - void add(const double *z, const Vector &v); - void add(const double *z, const ConstVector &v); + // Computes this += z*v (where this and v are intepreted as complex vectors) + void addComplex(const std::complex &z, const Vector &v); + // Computes this += z*v (where this and v are intepreted as complex vectors) + void addComplex(const std::complex &z, const ConstVector &v); void mult(double r); double getNorm() const; double getMax() const; @@ -133,31 +134,43 @@ public: } private: void copy(const double *d, int inc); - Vector &operator=(int); // must not be used (not implemented) - Vector &operator=(double); // must not be used (not implemented) }; -class BaseConstVector +class ConstGeneralMatrix; + +class ConstVector { protected: int len; - int s; - const double *data; + int off{0}; // offset to double* pointer + int s{1}; // stride (also called "skip" in some places) + std::shared_ptr data; public: - BaseConstVector(int l, int si, const double *d) : len(l), s(si), data(d) + // Implicit conversion from Vector is ok, since it's cheap + ConstVector(const Vector &v); + ConstVector(const ConstVector &v) = default; + ConstVector(ConstVector &&v) = default; + ConstVector(std::shared_ptr d, int l) : len{l}, data{std::move(d)} { } - BaseConstVector(const BaseConstVector &v) = default; - BaseConstVector &operator=(const BaseConstVector &v) = default; + ConstVector(const ConstVector &v, int off_arg, int l); + ConstVector(const ConstVector &v, int off_arg, int skip, int l); + ConstVector(std::shared_ptr d, int skip, int l); +#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE) + explicit ConstVector(const mxArray *p); +#endif + virtual ~ConstVector() = default; + ConstVector &operator=(const ConstVector &v) = delete; + ConstVector &operator=(ConstVector &&v) = delete; const double & operator[](int i) const { - return data[s*i]; + return data.get()[off+s*i]; } const double * base() const { - return data; + return data.get() + off; } int length() const @@ -169,28 +182,6 @@ public: { return s; } -}; - -class ConstGeneralMatrix; - -class ConstVector : public BaseConstVector -{ -public: - ConstVector(const Vector &v) : BaseConstVector(v.length(), v.skip(), v.base()) - { - } - ConstVector(const ConstVector &v) = default; - ConstVector(const double *d, int l) : BaseConstVector(l, 1, d) - { - } - ConstVector(const Vector &v, int off, int l); - ConstVector(const ConstVector &v, int off, int l); - ConstVector(const double *d, int skip, int l); - ConstVector(const ConstGeneralMatrix &m, int col); - ConstVector(int row, const ConstGeneralMatrix &m); - - virtual ~ConstVector() = default; - ConstVector &operator=(const ConstVector &v) = default; /** Exact equality. */ bool operator==(const ConstVector &y) const; bool diff --git a/dynare++/sylv/matlab/gensylv.cc b/dynare++/sylv/matlab/gensylv.cc index 2cd7691b2..92bd0481e 100644 --- a/dynare++/sylv/matlab/gensylv.cc +++ b/dynare++/sylv/matlab/gensylv.cc @@ -8,23 +8,23 @@ void gen_sylv_solve(int order, int n, int m, int zero_cols, - const double *A, const double *B, - const double *C, double *X) + const ConstVector &A, const ConstVector &B, + const ConstVector &C, Vector &X) { GeneralSylvester s(order, n, m, zero_cols, A, B, C, X, false); s.solve(); } -void +mxArray * gen_sylv_solve_and_check(int order, int n, int m, int zero_cols, - const double *A, const double *B, - const double *C, const double *D, - double *X, mxArray * &p) + const ConstVector &A, const ConstVector &B, + const ConstVector &C, const ConstVector &D, + Vector &X) { GeneralSylvester s(order, n, m, zero_cols, A, B, C, X, true); s.solve(); s.check(D); - p = s.getParams().createStructArray(); + return s.getParams().createStructArray(); } extern "C" { @@ -35,7 +35,7 @@ extern "C" { if (nhrs != 5 || nlhs > 3 || nlhs < 2) DYN_MEX_FUNC_ERR_MSG_TXT("Gensylv: Must have exactly 5 input args and either 2 or 3 output args."); - auto order = (int) mxGetScalar(prhs[0]); + auto order = static_cast(mxGetScalar(prhs[0])); const mxArray *const A = prhs[1]; const mxArray *const B = prhs[2]; const mxArray *const C = prhs[3]; @@ -55,28 +55,24 @@ extern "C" { DYN_MEX_FUNC_ERR_MSG_TXT("Matrix C must be square."); if (Bdims[0] < Bdims[1]) DYN_MEX_FUNC_ERR_MSG_TXT("Matrix B must not have more columns than rows."); - if (Ddims[1] != (mwSize) power(Cdims[0], order)) + if (Ddims[1] != static_cast(power(Cdims[0], order))) DYN_MEX_FUNC_ERR_MSG_TXT("Matrix D has wrong number of columns."); - int n = Adims[0]; - int m = Cdims[0]; - int zero_cols = Bdims[0] - Bdims[1]; + auto n = static_cast(Adims[0]); + auto m = static_cast(Cdims[0]); + auto zero_cols = static_cast(Bdims[0] - Bdims[1]); mxArray *X = mxCreateDoubleMatrix(Ddims[0], Ddims[1], mxREAL); // copy D to X - Vector Xvec((double *)mxGetPr(X), power(m, order)*n); - ConstVector Dvec((double *)mxGetPr(D), power(m, order)*n); + ConstVector Avec{A}, Bvec{B}, Cvec{C}, Dvec{D}; + Vector Xvec{X}; Xvec = Dvec; // solve (or solve and check) 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, Avec, Bvec, Cvec, Xvec); 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]); + plhs[2] = gen_sylv_solve_and_check(order, n, m, zero_cols, Avec, Bvec, Cvec, Dvec, Xvec); } catch (const SylvException &e) { diff --git a/dynare++/sylv/testing/MMMatrix.cc b/dynare++/sylv/testing/MMMatrix.cc index 8a6f52823..d9c425990 100644 --- a/dynare++/sylv/testing/MMMatrix.cc +++ b/dynare++/sylv/testing/MMMatrix.cc @@ -24,12 +24,12 @@ MMMatrixIn::MMMatrixIn(const char *fname) if (2 != sscanf(buffer, "%d %d", &rows, &cols)) throw MMException("Couldn't parse rows and cols\n"); // read in data - data = (double *) operator new[](rows*cols*sizeof(double)); + 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) { - if (1 != sscanf(buffer, "%lf", &data[i])) + if (1 != sscanf(buffer, "%lf", const_cast(data.get())+i)) throw MMException(string("Couldn't parse float number ")+buffer+"\n"); i++; } @@ -42,11 +42,6 @@ MMMatrixIn::MMMatrixIn(const char *fname) fclose(fd); } -MMMatrixIn::~MMMatrixIn() -{ - operator delete [](data); -} - void MMMatrixOut::write(const char *fname, int rows, int cols, const double *data) { diff --git a/dynare++/sylv/testing/MMMatrix.hh b/dynare++/sylv/testing/MMMatrix.hh index 542b7931c..fc44ce082 100644 --- a/dynare++/sylv/testing/MMMatrix.hh +++ b/dynare++/sylv/testing/MMMatrix.hh @@ -10,6 +10,7 @@ #include #include +#include using namespace std; @@ -32,16 +33,16 @@ public: class MMMatrixIn : public MallocAllocator { - double *data; + std::shared_ptr data; int rows; int cols; public: MMMatrixIn(const char *fname); - ~MMMatrixIn(); - const double * + ~MMMatrixIn() = default; + ConstVector getData() const { - return data; + return ConstVector{data, size()}; } int size() const diff --git a/dynare++/sylv/testing/tests.cc b/dynare++/sylv/testing/tests.cc index 3b8e4d9bb..2e768a864 100644 --- a/dynare++/sylv/testing/tests.cc +++ b/dynare++/sylv/testing/tests.cc @@ -33,6 +33,7 @@ public: { strncpy(name, n, 100); } + virtual ~TestRunnable() = default; bool test() const; virtual bool run() const = 0; const char * @@ -117,7 +118,7 @@ TestRunnable::quasi_solve(bool trans, const char *mname, const char *vname) printf(" Wrong quasi triangular dimensions, rows must be >= cols.\n"); return false; } - ConstVector v(mmv.getData(), mmv.row()); + ConstVector v{mmv.getData()}; Vector x(v.length()); double eig_min = 1.0e20; if (trans) @@ -158,9 +159,9 @@ TestRunnable::mult_kron(bool trans, const char *mname, const char *vname, SylvMemoryDriver memdriver(1, m, n, depth); QuasiTriangular t(mmt.getData(), mmt.row()); - Vector vraw(mmv.getData(), mmv.row()); + Vector vraw{mmv.getData()}; KronVector v(vraw, m, n, depth); - Vector craw(mmc.getData(), mmc.row()); + Vector craw{mmc.getData()}; KronVector c(craw, m, n, depth); if (trans) t.multKronTrans(v); @@ -192,9 +193,9 @@ TestRunnable::level_kron(bool trans, const char *mname, const char *vname, SylvMemoryDriver memdriver(1, m, n, depth); QuasiTriangular t(mmt.getData(), mmt.row()); - Vector vraw(mmv.getData(), mmv.row()); + Vector vraw{mmv.getData()}; ConstKronVector v(vraw, m, n, depth); - Vector craw(mmc.getData(), mmc.row()); + Vector craw{mmc.getData()}; KronVector c(craw, m, n, depth); KronVector x(v); if (trans) @@ -229,9 +230,9 @@ TestRunnable::kron_power(const char *m1name, const char *m2name, const char *vna SylvMemoryDriver memdriver(2, m, n, depth); QuasiTriangular t1(mmt1.getData(), mmt1.row()); QuasiTriangular t2(mmt2.getData(), mmt2.row()); - Vector vraw(mmv.getData(), mmv.row()); + Vector vraw{mmv.getData()}; ConstKronVector v(vraw, m, n, depth); - Vector craw(mmc.getData(), mmc.row()); + Vector craw{mmc.getData()}; KronVector c(craw, m, n, depth); KronVector x(v); memdriver.setStackMode(true); @@ -267,14 +268,14 @@ TestRunnable::lin_eval(const char *m1name, const char *m2name, const char *vname QuasiTriangular t1(mmt1.getData(), mmt1.row()); QuasiTriangular t2(mmt2.getData(), mmt2.row()); TriangularSylvester ts(t2, t1); - Vector vraw1(mmv.getData(), length); + ConstVector vraw1{mmv.getData(), 0, length}; ConstKronVector v1(vraw1, m, n, depth); - Vector vraw2(mmv.getData()+length, length); + ConstVector vraw2{mmv.getData(), length, length}; ConstKronVector v2(vraw2, m, n, depth); - Vector craw1(mmc.getData(), length); - KronVector c1(craw1, m, n, depth); - Vector craw2(mmc.getData()+length, length); - KronVector c2(craw2, m, n, depth); + ConstVector craw1{mmc.getData(), 0, length}; + ConstKronVector c1(craw1, m, n, depth); + ConstVector craw2{mmc.getData(), length, length}; + ConstKronVector c2(craw2, m, n, depth); KronVector x1(m, n, depth); KronVector x2(m, n, depth); memdriver.setStackMode(true); @@ -313,14 +314,14 @@ TestRunnable::qua_eval(const char *m1name, const char *m2name, const char *vname QuasiTriangular t1(mmt1.getData(), mmt1.row()); QuasiTriangular t2(mmt2.getData(), mmt2.row()); TriangularSylvester ts(t2, t1); - Vector vraw1(mmv.getData(), length); + ConstVector vraw1{mmv.getData(), 0, length}; ConstKronVector v1(vraw1, m, n, depth); - Vector vraw2(mmv.getData()+length, length); + ConstVector vraw2{mmv.getData(), length, length}; ConstKronVector v2(vraw2, m, n, depth); - Vector craw1(mmc.getData(), length); - KronVector c1(craw1, m, n, depth); - Vector craw2(mmc.getData()+length, length); - KronVector c2(craw2, m, n, depth); + ConstVector craw1{mmc.getData(), 0, length}; + ConstKronVector c1(craw1, m, n, depth); + ConstVector craw2{mmc.getData(), length, length}; + ConstKronVector c2(craw2, m, n, depth); KronVector x1(m, n, depth); KronVector x2(m, n, depth); memdriver.setStackMode(true); @@ -356,7 +357,7 @@ TestRunnable::tri_sylv(const char *m1name, const char *m2name, const char *vname QuasiTriangular t1(mmt1.getData(), mmt1.row()); QuasiTriangular t2(mmt2.getData(), mmt2.row()); TriangularSylvester ts(t2, t1); - Vector vraw(mmv.getData(), length); + Vector vraw{mmv.getData()}; ConstKronVector v(vraw, m, n, depth); KronVector d(v); // copy of v SylvParams pars; @@ -458,7 +459,7 @@ TestRunnable::block_diag(const char *aname, double log10norm) int n = mma.row(); SylvMemoryDriver memdriver(3, n, n, 2); SqSylvMatrix orig(mma.getData(), n); - SimilarityDecomp dec(orig.base(), orig.numRows(), log10norm); + SimilarityDecomp dec(orig.getData(), orig.numRows(), log10norm); dec.getB().printInfo(); SqSylvMatrix check(dec.getQ(), dec.getB()); check.multRight(dec.getInvQ()); @@ -506,7 +507,7 @@ TestRunnable::iter_sylv(const char *m1name, const char *m2name, const char *vnam QuasiTriangular t1(mmt1.getData(), mmt1.row()); QuasiTriangular t2(mmt2.getData(), mmt2.row()); IterativeSylvester is(t2, t1); - Vector vraw(mmv.getData(), length); + Vector vraw{mmv.getData()}; ConstKronVector v(vraw, m, n, depth); KronVector d(v); // copy of v SylvParams pars; diff --git a/dynare++/tl/cc/fs_tensor.cc b/dynare++/tl/cc/fs_tensor.cc index 08e389f17..13e29450b 100644 --- a/dynare++/tl/cc/fs_tensor.cc +++ b/dynare++/tl/cc/fs_tensor.cc @@ -181,7 +181,7 @@ UFSTensor::UFSTensor(const UFSTensor &t, const ConstVector &x) for (int i = 0; i < ncols(); i++) { ConstTwoDMatrix tpart(t, i *nvar(), nvar()); - Vector outcol(*this, i); + Vector outcol{getCol(i)}; tpart.multaVec(outcol, x); } } diff --git a/dynare++/tl/cc/kron_prod.cc b/dynare++/tl/cc/kron_prod.cc index 0eb9a65aa..d15f746a5 100644 --- a/dynare++/tl/cc/kron_prod.cc +++ b/dynare++/tl/cc/kron_prod.cc @@ -182,8 +182,8 @@ KronProdAI::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const if (in.getLD() == in.nrows()) { - ConstTwoDMatrix in_resh(in.nrows()*id_cols, a.nrows(), in.getData().base()); - TwoDMatrix out_resh(in.nrows()*id_cols, a.ncols(), out.getData().base()); + ConstTwoDMatrix in_resh(in.nrows()*id_cols, a.nrows(), in.getData()); + TwoDMatrix out_resh(in.nrows()*id_cols, a.ncols(), out.getData()); out_resh.mult(in_resh, a); } else @@ -293,7 +293,7 @@ KronProdAll::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const } else { - last = new TwoDMatrix(in.nrows(), in.ncols(), in.getData().base()); + last = new TwoDMatrix(in.nrows(), in.ncols(), in.getData()); } // perform intermediate multiplications IAI @@ -352,7 +352,7 @@ KronProdAll::multRows(const IntSequence &irows) const the |row| as ConstVector of this vector, which sheduled for deallocation. */ if (matlist[j]) - row = new ConstVector(irows[j], *(matlist[j])); + row = new ConstVector(matlist[j]->getRow(irows[j])); else { auto *aux = new Vector(ncols(j)); diff --git a/dynare++/tl/cc/ps_tensor.cc b/dynare++/tl/cc/ps_tensor.cc index a97dd9681..b72118ceb 100644 --- a/dynare++/tl/cc/ps_tensor.cc +++ b/dynare++/tl/cc/ps_tensor.cc @@ -390,7 +390,7 @@ FPSTensor::FPSTensor(const TensorDimens &td, const Equivalence &e, const Permuta auto su = a.getMap().upper_bound(c); for (auto srun = sl; srun != su; ++srun) { - Vector out_row((*srun).second.first, *this); + Vector out_row{getRow((*srun).second.first)}; out_row.add((*srun).second.second, *row_prod); } delete row_prod; diff --git a/dynare++/tl/cc/pyramid_prod.cc b/dynare++/tl/cc/pyramid_prod.cc index 30eb37faf..4cc5c8af7 100644 --- a/dynare++/tl/cc/pyramid_prod.cc +++ b/dynare++/tl/cc/pyramid_prod.cc @@ -68,10 +68,10 @@ USubTensor::addKronColumn(int i, const vector &ts, IntSequence ind(pindex, lastdim, lastdim+t->dimen()); lastdim += t->dimen(); index in(t, ind); - tmpcols.emplace_back(*t, *in); + tmpcols.push_back(t->getCol(*in)); } URSingleTensor kronmult(tmpcols); - Vector coli(*this, i); + Vector coli{getCol(i)}; coli.add(1.0, kronmult.getData()); } diff --git a/dynare++/tl/cc/stack_container.cc b/dynare++/tl/cc/stack_container.cc index 768f49fb8..6af9f5633 100644 --- a/dynare++/tl/cc/stack_container.cc +++ b/dynare++/tl/cc/stack_container.cc @@ -216,7 +216,7 @@ FoldedStackContainer::multAndAddSparse3(const FSSparseTensor &t, const EquivalenceSet &eset = ebundle.get(out.dimen()); for (Tensor::index run = out.begin(); run != out.end(); ++run) { - Vector outcol(out, *run); + Vector outcol{out.getCol(*run)}; FRSingleTensor sumcol(t.nvar(), t.dimen()); sumcol.zeros(); for (const auto & it : eset) diff --git a/dynare++/tl/cc/stack_container.hh b/dynare++/tl/cc/stack_container.hh index 742d40cf5..074eb6900 100644 --- a/dynare++/tl/cc/stack_container.hh +++ b/dynare++/tl/cc/stack_container.hh @@ -249,7 +249,7 @@ public: const _Ttype *t = getMatrix(i, s); Tensor::index ind(t, coor); Vector subres(*res, stack_offsets[i], stack_sizes[i]); - subres = ConstVector(ConstGeneralMatrix(*t), *ind); + subres = ConstGeneralMatrix(*t).getCol(*ind); i++; } if (iu != -1) diff --git a/dynare++/tl/cc/twod_matrix.cc b/dynare++/tl/cc/twod_matrix.cc index 38d1e885d..a88896ff5 100644 --- a/dynare++/tl/cc/twod_matrix.cc +++ b/dynare++/tl/cc/twod_matrix.cc @@ -64,17 +64,13 @@ TwoDMatrix::copyRow(int from, int to) void TwoDMatrix::copyRow(const ConstTwoDMatrix &m, int from, int to) { - ConstVector fr_row(from, m); - Vector to_row(to, *this); - to_row = fr_row; + getRow(to) = m.getRow(from); } void TwoDMatrix::addRow(double d, const ConstTwoDMatrix &m, int from, int to) { - ConstVector fr_row(from, m); - Vector to_row(to, *this); - to_row.add(d, fr_row); + getRow(to).add(d, m.getRow(from)); } void @@ -87,17 +83,13 @@ TwoDMatrix::copyColumn(int from, int to) void TwoDMatrix::copyColumn(const ConstTwoDMatrix &m, int from, int to) { - ConstVector fr_col(m, from); - Vector to_col(*this, to); - to_col = fr_col; + getCol(to) = m.getCol(from); } void TwoDMatrix::addColumn(double d, const ConstTwoDMatrix &m, int from, int to) { - ConstVector fr_col(m, from); - Vector to_col(*this, to); - to_col.add(d, fr_col); + getCol(to).add(d, m.getCol(from)); } void diff --git a/dynare++/tl/cc/twod_matrix.hh b/dynare++/tl/cc/twod_matrix.hh index 2f2435c9d..43c0f3bc8 100644 --- a/dynare++/tl/cc/twod_matrix.hh +++ b/dynare++/tl/cc/twod_matrix.hh @@ -19,6 +19,7 @@ #include #include +#include class TwoDMatrix; @@ -29,10 +30,11 @@ class TwoDMatrix; class ConstTwoDMatrix : public ConstGeneralMatrix { public: - ConstTwoDMatrix(int m, int n, const double *d) - : ConstGeneralMatrix(d, m, n) + ConstTwoDMatrix(int m, int n, ConstVector d) + : ConstGeneralMatrix(std::move(d), m, n) { } + // Implicit conversion from TwoDMatrix is ok, since it's cheap ConstTwoDMatrix(const TwoDMatrix &m); ConstTwoDMatrix(const TwoDMatrix &m, int first_col, int num); ConstTwoDMatrix(const ConstTwoDMatrix &m, int first_col, int num); @@ -71,11 +73,11 @@ public: : GeneralMatrix(r, c) { } - TwoDMatrix(int r, int c, double *d) + TwoDMatrix(int r, int c, Vector &d) : GeneralMatrix(d, r, c) { } - TwoDMatrix(int r, int c, const double *d) + TwoDMatrix(int r, int c, const ConstVector &d) : GeneralMatrix(d, r, c) { } @@ -83,6 +85,11 @@ public: : GeneralMatrix(m) { } + // We don't want implict conversion from ConstGeneralMatrix, since it's expensive + explicit TwoDMatrix(const ConstGeneralMatrix &m) + : GeneralMatrix(m) + { + } TwoDMatrix(const GeneralMatrix &m, const char *dummy) : GeneralMatrix(m, dummy) { diff --git a/dynare++/tl/testing/monoms.cc b/dynare++/tl/testing/monoms.cc index bf3977dde..3fbe378c0 100644 --- a/dynare++/tl/testing/monoms.cc +++ b/dynare++/tl/testing/monoms.cc @@ -124,7 +124,7 @@ Monom1Vector::deriv(int dim) const = new FGSTensor(len, TensorDimens(Symmetry(dim), IntSequence(1, nx))); for (Tensor::index it = res->begin(); it != res->end(); ++it) { - Vector outcol(*res, *it); + Vector outcol{res->getCol(*it)}; deriv(it.getCoor(), outcol); } return res; @@ -214,7 +214,7 @@ Monom2Vector::deriv(const Symmetry &s) const FGSTensor *t = new FGSTensor(len, TensorDimens(s, nvs)); for (Tensor::index it = t->begin(); it != t->end(); ++it) { - Vector col(*t, *it); + Vector col{t->getCol(*it)}; deriv(s, it.getCoor(), col); } return t; @@ -392,7 +392,7 @@ Monom4Vector::deriv(const Symmetry &s) const FGSTensor *res = new FGSTensor(len, TensorDimens(s, nvs)); for (Tensor::index run = res->begin(); run != res->end(); ++run) { - Vector col(*res, *run); + Vector col{res->getCol(*run)}; deriv(s, run.getCoor(), col); } return res; diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc index 162e98374..45c720d55 100644 --- a/mex/sources/k_order_perturbation/dynamic_m.cc +++ b/mex/sources/k_order_perturbation/dynamic_m.cc @@ -45,7 +45,7 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam if (retVal != 0) throw DynareException(__FILE__, __LINE__, "Trouble calling " + DynamicMFilename); - residual = Vector(mxGetPr(plhs[0]), residual.skip(), (int) mxGetM(plhs[0])); + residual = Vector{plhs[0]}; copyDoubleIntoTwoDMatData(mxGetPr(plhs[1]), g1, (int) mxGetM(plhs[1]), (int) mxGetN(plhs[1])); if (g2 != nullptr) unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2); diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc index 4b95563fa..6da182bb8 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cc +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc @@ -89,15 +89,15 @@ KordpDynare::solveDeterministicSteady() } void -KordpDynare::evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) noexcept(false) +KordpDynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) noexcept(false) { // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals out.zeros(); } void -KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, - const Vector &yyp, const Vector &xx) noexcept(false) +KordpDynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy, + const ConstVector &yyp, const Vector &xx) noexcept(false) { // This method is only called when checking the residuals at steady state (Approximation::check), so return zero residuals out.zeros(); diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh index 93b1d90b7..35d00e20f 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.hh +++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh @@ -230,9 +230,9 @@ public: } void solveDeterministicSteady(); - void evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) noexcept(false); - void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, - const Vector &yyp, const Vector &xx) noexcept(false); + void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) noexcept(false); + void evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy, + const ConstVector &yyp, const Vector &xx) noexcept(false); void calcDerivativesAtSteady(); unique_ptr dynamicModelFile; DynamicModel * diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc index 79a921a67..c59fe8325 100644 --- a/mex/sources/k_order_perturbation/k_order_perturbation.cc +++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc @@ -111,23 +111,18 @@ extern "C" { qz_criterium = (double) mxGetScalar(mxFldp); mxFldp = mxGetField(M_, 0, "params"); - double *dparams = mxGetPr(mxFldp); - int npar = (int) mxGetM(mxFldp); - Vector modParams(dparams, npar); + Vector modParams{mxFldp}; if (!modParams.isFinite()) DYN_MEX_FUNC_ERR_MSG_TXT("The parameters vector contains NaN or Inf"); mxFldp = mxGetField(M_, 0, "Sigma_e"); - dparams = mxGetPr(mxFldp); - npar = (int) mxGetN(mxFldp); - TwoDMatrix vCov(npar, npar, dparams); + int npar = static_cast(mxGetN(mxFldp)); + TwoDMatrix vCov(npar, npar, Vector{mxFldp}); if (!vCov.isFinite()) DYN_MEX_FUNC_ERR_MSG_TXT("The covariance matrix of shocks contains NaN or Inf"); mxFldp = mxGetField(dr, 0, "ys"); // and not in order of dr.order_var - dparams = mxGetPr(mxFldp); - const int nSteady = (int) mxGetM(mxFldp); - Vector ySteady(dparams, nSteady); + Vector ySteady{mxFldp}; if (!ySteady.isFinite()) DYN_MEX_FUNC_ERR_MSG_TXT("The steady state vector contains NaN or Inf"); @@ -152,7 +147,7 @@ extern "C" { const int nPar = (int) mxGetScalar(mxFldp); mxFldp = mxGetField(dr, 0, "order_var"); - dparams = mxGetPr(mxFldp); + auto dparams = mxGetPr(mxFldp); npar = (int) mxGetM(mxFldp); if (npar != nEndo) DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect number of input var_order vars."); @@ -163,11 +158,10 @@ extern "C" { // the lag, current and lead blocks of the jacobian respectively mxFldp = mxGetField(M_, 0, "lead_lag_incidence"); - dparams = mxGetPr(mxFldp); npar = (int) mxGetN(mxFldp); int nrows = (int) mxGetM(mxFldp); - TwoDMatrix llincidence(nrows, npar, dparams); + TwoDMatrix llincidence(nrows, npar, Vector{mxFldp}); if (npar != nEndo) { ostringstream strstrm; @@ -176,8 +170,7 @@ extern "C" { } //get NNZH =NNZD(2) = the total number of non-zero Hessian elements mxFldp = mxGetField(M_, 0, "NNZDerivatives"); - dparams = mxGetPr(mxFldp); - Vector NNZD(dparams, (int)mxGetM(mxFldp)); + Vector NNZD{mxFldp}; if (NNZD[kOrder-1] == -1) DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the stoch_simul command"); @@ -207,19 +200,19 @@ extern "C" { const mxArray *g1 = prhs[3]; int m = (int) mxGetM(g1); int n = (int) mxGetN(g1); - g1m = new TwoDMatrix(m, n, mxGetPr(g1)); + g1m = new TwoDMatrix(m, n, ConstVector{g1}); if (nrhs > 4) { const mxArray *g2 = prhs[4]; int m = (int) mxGetM(g2); int n = (int) mxGetN(g2); - g2m = new TwoDMatrix(m, n, mxGetPr(g2)); + g2m = new TwoDMatrix(m, n, ConstVector{g2}); if (nrhs > 5) { const mxArray *g3 = prhs[5]; int m = (int) mxGetM(g3); int n = (int) mxGetN(g3); - g3m = new TwoDMatrix(m, n, mxGetPr(g3)); + g3m = new TwoDMatrix(m, n, ConstVector{g3}); } } }