Commit c711d34d authored by Sébastien Villemot's avatar Sébastien Villemot

Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes

- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
  they do not perform memory management, and several {Const,}Vector instances
  can transparently share the same underlying data

- make converting constructor from ConstVector to Vector explicit, since that
  entails memory allocation (but the reverse conversion is almost costless, so
  keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
  TwoDMatrix/ConstTwoDMatrix

- remove the constructors that were extracting a row/column from a matrix, and
  replace them by getRow() and getCol() methods on {Const,}GeneralMatrix

- rename and change the API of the complex version Vector::add(), so that it is
  explicit that it deals with complex numbers

- add constructors that take a MATLAB mxArray
parent 49c06f5c
......@@ -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;
......
......@@ -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;
......
......@@ -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;
}
......
......@@ -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<t>::Tpol(pol), ysteady(ys), ypart(yp), nu(nuu)
{
}
DecisionRuleImpl(_Tparent &pol, const PartitionY &yp, int nuu,
const Vector &ys)
const ConstVector &ys)
: ctraits<t>::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<t>::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<t>::centralize(const DecisionRuleImpl &dr)
template <int t>
TwoDMatrix *
DecisionRuleImpl<t>::simulate(emethod em, int np, const Vector &ystart,
DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
ShockRealization &sr) const
{
KORD_RAISE_IF(ysteady.length() != ystart.length(),
......@@ -303,7 +303,7 @@ DecisionRuleImpl<t>::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<t>::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<t>::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<KOrder::fold>
friend class UnfoldDecisionRule;
public:
FoldDecisionRule(const ctraits<KOrder::fold>::Tpol &pol, const PartitionY &yp, int nuu,
const Vector &ys)
const ConstVector &ys)
: DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys)
{
}
FoldDecisionRule(ctraits<KOrder::fold>::Tpol &pol, const PartitionY &yp, int nuu,
const Vector &ys)
const ConstVector &ys)
: DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys)
{
}
FoldDecisionRule(const ctraits<KOrder::fold>::Tg &g, const PartitionY &yp, int nuu,
const Vector &ys, double sigma)
const ConstVector &ys, double sigma)
: DecisionRuleImpl<KOrder::fold>(g, yp, nuu, ys, sigma)
{
}
......@@ -445,17 +445,17 @@ class UnfoldDecisionRule : public DecisionRuleImpl<KOrder::unfold>
friend class FoldDecisionRule;
public:
UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
const Vector &ys)
const ConstVector &ys)
: DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys)
{
}
UnfoldDecisionRule(ctraits<KOrder::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
const Vector &ys)
const ConstVector &ys)
: DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys)
{
}
UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tg &g, const PartitionY &yp, int nuu,
const Vector &ys, double sigma)
const ConstVector &ys, double sigma)
: DecisionRuleImpl<KOrder::unfold>(g, yp, nuu, ys, sigma)
{
}
......
......@@ -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;
};
......
......@@ -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);
}
......
......@@ -84,7 +84,7 @@ public:
return std::make_unique<ResidFunction>(*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<GResidFunction>(*this);
}
void
setYU(const Vector &ys, const Vector &xx)
setYU(const ConstVector &ys, const ConstVector &xx)
{
((ResidFunction *) func)->setYU(ys, xx);
}
......
......@@ -307,8 +307,8 @@ KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
TwoDMatrix gs_y(*(gs<unfold>().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)
......
......@@ -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);
}
......
......@@ -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<const double>{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,
......
......@@ -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);
......
......@@ -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;
......
......@@ -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)
{
}
......
......@@ -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);
}
......
......@@ -6,9 +6,10 @@
#include <iostream>
#include <cstring>
#include <utility>
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++)
......
......@@ -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);
......
......@@ -16,8 +16,6 @@
#include <limits>
#include <vector>
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);
}
......@@ -8,6 +8,8 @@
#include "Vector.hh"
#include <algorithm>
#include <memory>
#include <utility>
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)
{
}