Commit 449553e7 authored by Sébastien Villemot's avatar Sébastien Villemot

Dynare++: various improvements

parent 4c11e9e9
......@@ -112,11 +112,11 @@ SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
// todo: check |l>1|, |l>=d|
// todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()|
int cum = 0;
for (auto &si : SymmetrySet(l-1, d+1))
for (const auto &si : SymmetrySet(l-1, d+1))
{
if (si[d] <= d-1)
{
IntSequence lev((const IntSequence &) si, 0, d);
IntSequence lev(si, 0, d);
lev.add(1);
levels.push_back(lev);
IntSequence levpts(d);
......@@ -170,11 +170,11 @@ int
SmolyakQuadrature::calcNumEvaluations(int lev) const
{
int cum = 0;
for (auto &si : SymmetrySet(lev-1, dim+1))
for (const auto &si : SymmetrySet(lev-1, dim+1))
{
if (si[dim] <= dim-1)
{
IntSequence lev((const IntSequence &) si, 0, dim);
IntSequence lev(si, 0, dim);
lev.add(1);
IntSequence levpts(dim);
for (int i = 0; i < dim; i++)
......
......@@ -127,7 +127,7 @@ Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
for (int i = 0; i < dim; i++)
r *= point[i];
r = pow(r, 1.0/dim);
r *= pow(1.0 + 1.0/dim, (double) dim);
r *= pow(1.0 + 1.0/dim, static_cast<double>(dim));
out[0] = r;
}
......@@ -377,7 +377,9 @@ public:
run() const override
{
GeneralMatrix m(2, 2);
m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
m.zeros();
m.get(0, 0) = 1;
m.get(1, 1) = 1;
return smolyak_normal_moments(m, 4, 4);
}
};
......@@ -395,8 +397,12 @@ public:
{
GeneralMatrix m(3, 3);
m.zeros();
m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
m.get(0, 0) = 1;
m.get(0, 2) = 0.5;
m.get(1, 1) = 1;
m.get(1, 0) = 0.5;
m.get(2, 2) = 2;
m.get(2, 1) = 4;
return smolyak_normal_moments(m, 8, 8);
}
};
......@@ -413,7 +419,9 @@ public:
run() const override
{
GeneralMatrix m(2, 2);
m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
m.zeros();
m.get(0, 0) = 1;
m.get(1, 1) = 1;
return product_normal_moments(m, 4, 4);
}
};
......@@ -431,8 +439,12 @@ public:
{
GeneralMatrix m(3, 3);
m.zeros();
m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
m.get(0, 0) = 1;
m.get(0, 2) = 0.5;
m.get(1, 1) = 1;
m.get(1, 0) = 0.5;
m.get(2, 2) = 2;
m.get(2, 1) = 4;
return product_normal_moments(m, 8, 8);
}
};
......@@ -450,7 +462,8 @@ public:
run() const override
{
Function1Trans f1(6);
Vector res(1); res[0] = 1.0;
Vector res(1);
res[0] = 1.0;
return smolyak_product_cube(f1, res, 1e-2, 13);
}
};
......@@ -488,10 +501,8 @@ main()
int nvmax = 0;
for (const auto &test : all_tests)
{
if (dmax < test->dim)
dmax = test->dim;
if (nvmax < test->nvar)
nvmax = test->nvar;
dmax = std::max(dmax, test->dim);
nvmax = std::max(nvmax, test->nvar);
}
TLStatic::init(dmax, nvmax); // initialize library
......
......@@ -41,7 +41,7 @@ BlockDiagonal::BlockDiagonal(int p, const BlockDiagonal &b)
void
BlockDiagonal::setZerosToRU(diag_iter edge)
{
int iedge = (*edge).getIndex();
int iedge = edge->getIndex();
for (int i = 0; i < iedge; i++)
for (int j = iedge; j < numCols(); j++)
get(i, j) = 0.0;
......@@ -61,24 +61,24 @@ BlockDiagonal::setZeroBlockEdge(diag_iter edge)
{
setZerosToRU(edge);
int iedge = (*edge).getIndex();
int iedge = edge->getIndex();
for (diag_iter run = diag_begin(); run != edge; ++run)
{
int ind = (*run).getIndex();
int ind = run->getIndex();
if (row_len[ind] > iedge)
{
row_len[ind] = iedge;
if (!(*run).isReal())
if (!run->isReal())
row_len[ind+1] = iedge;
}
}
for (diag_iter run = edge; run != diag_end(); ++run)
{
int ind = (*run).getIndex();
int ind = run->getIndex();
if (col_len[ind] < iedge)
{
col_len[ind] = iedge;
if (!(*run).isReal())
if (!run->isReal())
col_len[ind+1] = iedge;
}
}
......@@ -136,7 +136,7 @@ BlockDiagonal::findBlockStart(const_diag_iter from) const
{
++from;
while (from != diag_end()
&& col_len[(*from).getIndex()] != (*from).getIndex())
&& col_len[from->getIndex()] != from->getIndex())
++from;
}
return from;
......@@ -150,12 +150,11 @@ BlockDiagonal::getLargestBlock() const
const_diag_iter end = findBlockStart(start);
while (start != diag_end())
{
int si = (*start).getIndex();
int si = start->getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
if (largest < ei-si)
largest = ei-si;
ei = end->getIndex();
largest = std::max(largest, ei-si);
start = end;
end = findBlockStart(start);
}
......@@ -177,21 +176,21 @@ void
BlockDiagonal::multKronBlock(const_diag_iter start, const_diag_iter end,
KronVector &x, Vector &work) const
{
int si = (*start).getIndex();
int si = start->getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
ei = end->getIndex();
savePartOfX(si, ei, x, work);
for (const_diag_iter di = start; di != end; ++di)
{
int jbar = (*di).getIndex();
if ((*di).isReal())
int jbar = di->getIndex();
if (di->isReal())
{
KronVector xi(x, jbar);
xi.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
xi.add(*(di->getAlpha()), wi);
for (const_row_iter ri = row_begin(*di); ri != row_end(*di); ++ri)
{
int col = ri.getCol();
......@@ -207,10 +206,10 @@ BlockDiagonal::multKronBlock(const_diag_iter start, const_diag_iter end,
xii.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
Vector wii(work, (jbar+1-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
xi.add((*di).getBeta1(), wii);
xii.add((*di).getBeta2(), wi);
xii.add(*((*di).getAlpha()), wii);
xi.add(*(di->getAlpha()), wi);
xi.add(di->getBeta1(), wii);
xii.add(di->getBeta2(), wi);
xii.add(*(di->getAlpha()), wii);
for (const_row_iter ri = row_begin(*di); ri != row_end(*di); ++ri)
{
int col = ri.getCol();
......@@ -226,21 +225,21 @@ void
BlockDiagonal::multKronBlockTrans(const_diag_iter start, const_diag_iter end,
KronVector &x, Vector &work) const
{
int si = (*start).getIndex();
int si = start->getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
ei = end->getIndex();
savePartOfX(si, ei, x, work);
for (const_diag_iter di = start; di != end; ++di)
{
int jbar = (*di).getIndex();
if ((*di).isReal())
int jbar = di->getIndex();
if (di->isReal())
{
KronVector xi(x, jbar);
xi.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
xi.add(*(di->getAlpha()), wi);
for (const_col_iter ci = col_begin(*di); ci != col_end(*di); ++ci)
{
int row = ci.getRow();
......@@ -256,10 +255,10 @@ BlockDiagonal::multKronBlockTrans(const_diag_iter start, const_diag_iter end,
xii.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
Vector wii(work, (jbar+1-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
xi.add((*di).getBeta2(), wii);
xii.add((*di).getBeta1(), wi);
xii.add(*((*di).getAlpha()), wii);
xi.add(*(di->getAlpha()), wi);
xi.add(di->getBeta2(), wii);
xii.add(di->getBeta1(), wi);
xii.add(*(di->getAlpha()), wii);
for (const_col_iter ci = col_begin(*di); ci != col_end(*di); ++ci)
{
int row = ci.getRow();
......@@ -310,10 +309,10 @@ BlockDiagonal::printInfo() const
const_diag_iter end = findBlockStart(start);
while (start != diag_end())
{
int si = (*start).getIndex();
int si = start->getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
ei = end->getIndex();
std::cout << ' ' << ei-si;
num_blocks++;
start = end;
......
......@@ -268,9 +268,7 @@ GeneralMatrix::gemm(const std::string &transa, const ConstGeneralMatrix &a,
if (opa_rows != numRows()
|| opb_cols != numCols()
|| opa_cols != opb_rows)
{
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix multiplication.");
}
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix multiplication.");
blas_int m = opa_rows;
blas_int n = opb_cols;
......@@ -279,10 +277,8 @@ GeneralMatrix::gemm(const std::string &transa, const ConstGeneralMatrix &a,
blas_int ldb = b.ld;
blas_int ldc = ld;
if (lda > 0 && ldb > 0 && ldc > 0)
{
dgemm(transa.c_str(), transb.c_str(), &m, &n, &k, &alpha, a.data.base(), &lda,
b.data.base(), &ldb, &beta, data.base(), &ldc);
}
dgemm(transa.c_str(), transb.c_str(), &m, &n, &k, &alpha, a.data.base(), &lda,
b.data.base(), &ldb, &beta, data.base(), &ldc);
else if (numRows()*numCols() > 0)
{
if (beta == 0.0)
......@@ -299,14 +295,14 @@ GeneralMatrix::gemm_partial_left(const std::string &trans, const ConstGeneralMat
int icol;
for (icol = 0; icol + md_length < cols; icol += md_length)
{
GeneralMatrix incopy((const GeneralMatrix &)*this, 0, icol, rows, md_length);
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, md_length);
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), 0, icol, rows, md_length);
GeneralMatrix inplace(*this, 0, icol, rows, md_length);
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
}
if (cols > icol)
{
GeneralMatrix incopy((const GeneralMatrix &)*this, 0, icol, rows, cols - icol);
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, cols - icol);
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), 0, icol, rows, cols - icol);
GeneralMatrix inplace(*this, 0, icol, rows, cols - icol);
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
}
}
......@@ -318,14 +314,14 @@ GeneralMatrix::gemm_partial_right(const std::string &trans, const ConstGeneralMa
int irow;
for (irow = 0; irow + md_length < rows; irow += md_length)
{
GeneralMatrix incopy((const GeneralMatrix &)*this, irow, 0, md_length, cols);
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, md_length, cols);
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), irow, 0, md_length, cols);
GeneralMatrix inplace(*this, irow, 0, md_length, cols);
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
}
if (rows > irow)
{
GeneralMatrix incopy((const GeneralMatrix &)*this, irow, 0, rows - irow, cols);
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, rows - irow, cols);
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), irow, 0, rows - irow, cols);
GeneralMatrix inplace(*this, irow, 0, rows - irow, cols);
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
}
}
......@@ -366,8 +362,7 @@ ConstGeneralMatrix::getNormInf() const
for (int i = 0; i < numRows(); i++)
{
double normi = getRow(i).getNorm1();
if (norm < normi)
norm = normi;
norm = std::max(normi, norm);
}
return norm;
}
......@@ -379,8 +374,7 @@ ConstGeneralMatrix::getNorm1() const
for (int j = 0; j < numCols(); j++)
{
double normj = getCol(j).getNorm1();
if (norm < normj)
norm = normj;
norm = std::max(normj, norm);
}
return norm;
}
......@@ -548,7 +542,7 @@ SVDDecomp::construct(const GeneralMatrix &A)
// query for optimal lwork
dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &tmpwork,
&lwork, iwork.data(), &info);
lwork = (lapack_int) tmpwork;
lwork = static_cast<lapack_int>(tmpwork);
Vector work(lwork);
// do the decomposition
dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work.base(),
......@@ -567,7 +561,7 @@ SVDDecomp::solve(const ConstGeneralMatrix &B, GeneralMatrix &X) const
// reciprocal condition number for determination of zeros in the
// end of sigma
double rcond = 1e-13;
constexpr double rcond = 1e-13;
// determine nz=number of zeros in the end of sigma
int nz = 0;
......
......@@ -96,7 +96,7 @@ GeneralSylvester::solve()
d.multLeftI(bdecomp->getQ());
d.multRightKron(cdecomp->getInvQ(), order);
clock_t end = clock();
pars.cpu_time = ((double) (end-start))/CLOCKS_PER_SEC;
pars.cpu_time = static_cast<double>(end-start)/CLOCKS_PER_SEC;
solved = true;
}
......
......@@ -24,7 +24,7 @@ KronUtils::multAtLevel(int level, const QuasiTriangular &t,
}
else if (0 == level && 0 == x.getDepth())
{
Vector b((const Vector &)x);
Vector b(const_cast<const KronVector &>(x));
t.multVec(x, b);
}
else // 0 < level == depth
......@@ -36,13 +36,11 @@ KronUtils::multAtLevelTrans(int level, const QuasiTriangular &t,
KronVector &x)
{
if (0 < level && level < x.getDepth())
{
for (int i = 0; i < x.getM(); i++)
{
KronVector xi(x, i);
multAtLevelTrans(level, t, xi);
}
}
for (int i = 0; i < x.getM(); i++)
{
KronVector xi(x, i);
multAtLevelTrans(level, t, xi);
}
else if (0 == level && 0 < x.getDepth())
{
GeneralMatrix tmp(x, x.getN(), power(x.getM(), x.getDepth()));
......@@ -50,7 +48,7 @@ KronUtils::multAtLevelTrans(int level, const QuasiTriangular &t,
}
else if (level == 0 && 0 == x.getDepth())
{
Vector b((const Vector &)x);
Vector b(const_cast<const KronVector &>(x));
t.multVecTrans(x, b);
}
else // 0 < level == depth
......
......@@ -557,8 +557,7 @@ QuasiTriangular::solvePre(Vector &x, double &eig_min)
}
else
eig_size = *di->getAlpha()*(*di->getAlpha());
if (eig_size < eig_min)
eig_min = eig_size;
eig_min = std::min(eig_min, eig_size);
}
blas_int nn = diagonal.getSize();
......@@ -627,7 +626,7 @@ QuasiTriangular::multVecTrans(Vector &x, const ConstVector &b) const
void
QuasiTriangular::multaVec(Vector &x, const ConstVector &b) const
{
Vector tmp((const Vector &)x); // new copy
Vector tmp(const_cast<const Vector &>(x)); // new copy
multVec(x, b);
x.add(1.0, tmp);
}
......@@ -635,7 +634,7 @@ QuasiTriangular::multaVec(Vector &x, const ConstVector &b) const
void
QuasiTriangular::multaVecTrans(Vector &x, const ConstVector &b) const
{
Vector tmp((const Vector &)x); // new copy
Vector tmp(const_cast<const Vector &>(x)); // new copy
multVecTrans(x, b);
x.add(1.0, tmp);
}
......@@ -663,7 +662,7 @@ QuasiTriangular::multaKronTrans(KronVector &x, const ConstKronVector &b) const
void
QuasiTriangular::multKron(KronVector &x) const
{
KronVector b((const KronVector &)x); // make copy
KronVector b(const_cast<const KronVector &>(x)); // make copy
x.zeros();
multaKron(x, b);
}
......@@ -671,7 +670,7 @@ QuasiTriangular::multKron(KronVector &x) const
void
QuasiTriangular::multKronTrans(KronVector &x) const
{
KronVector b((const KronVector &)x); // make copy
KronVector b(const_cast<const KronVector &>(x)); // make copy
x.zeros();
multaKronTrans(x, b);
}
......
......@@ -41,7 +41,7 @@ public:
std::unique_ptr<QuasiTriangular>
clone(int p, const QuasiTriangular &t) const override
{
return std::make_unique<QuasiTriangularZero>(p, (const QuasiTriangularZero &) t);
return std::make_unique<QuasiTriangularZero>(p, dynamic_cast<const QuasiTriangularZero &>(t));
}
std::unique_ptr<QuasiTriangular>
clone(double r) const override
......@@ -51,7 +51,7 @@ public:
std::unique_ptr<QuasiTriangular>
clone(double r, double rr, const QuasiTriangular &tt) const override
{
return std::make_unique<QuasiTriangularZero>(r, *this, rr, (const QuasiTriangularZero &) tt);
return std::make_unique<QuasiTriangularZero>(r, *this, rr, dynamic_cast<const QuasiTriangularZero &>(tt));
}
void print() const override;
};
......
......@@ -53,8 +53,8 @@ SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
--itadd;
lapack_int n = getDim(), ldt = getT().getLD(), ldq = getQ().getLD();
lapack_int ifst = (*it).getIndex() + 1;
lapack_int ilst = (*itadd).getIndex() + 1;
lapack_int ifst = it->getIndex() + 1;
lapack_int ilst = itadd->getIndex() + 1;
std::vector<double> work(n);
lapack_int info;
dtrexc("V", &n, getT().base(), &ldt, getQ().base(), &ldq, &ifst, &ilst, work.data(),
......@@ -87,7 +87,7 @@ SchurDecompEig::orderEigen()
{
diag_iter least = getT().findNextLargerBlock(run, getT().diag_end(),
last_size);
last_size = (*least).getSize();
last_size = least->getSize();
if (run == least)
++run;
else
......
......@@ -27,8 +27,8 @@ void
SimilarityDecomp::getXDim(diag_iter start, diag_iter end,
int &rows, int &cols) const
{
int si = (*start).getIndex();
int ei = (*end).getIndex();
int si = start->getIndex();
int ei = end->getIndex();
cols = b->numRows() - ei;
rows = ei - si;
}
......@@ -41,12 +41,12 @@ bool
SimilarityDecomp::solveX(diag_iter start, diag_iter end,
GeneralMatrix &X, double norm) const
{
int si = (*start).getIndex();
int ei = (*end).getIndex();
int si = start->getIndex();
int ei = end->getIndex();
SqSylvMatrix A((const GeneralMatrix &)*b, si, si, X.numRows());
SqSylvMatrix B((const GeneralMatrix &)*b, ei, ei, X.numCols());
GeneralMatrix C((const GeneralMatrix &)*b, si, ei, X.numRows(), X.numCols());
SqSylvMatrix A(const_cast<const BlockDiagonal &>(*b), si, si, X.numRows());
SqSylvMatrix B(const_cast<const BlockDiagonal &>(*b), ei, ei, X.numCols());
GeneralMatrix C(const_cast<const BlockDiagonal &>(*b), si, ei, X.numRows(), X.numCols());
lapack_int isgn = -1;
lapack_int m = A.numRows();
......@@ -73,8 +73,8 @@ void
SimilarityDecomp::updateTransform(diag_iter start, diag_iter end,
GeneralMatrix &X)
{
int si = (*start).getIndex();
int ei = (*end).getIndex();
int si = start->getIndex();
int ei = end->getIndex();
SqSylvMatrix iX(q->numRows());
iX.setUnit();
......@@ -92,7 +92,7 @@ SimilarityDecomp::bringGuiltyBlock(diag_iter start, diag_iter &end)
{
double av = b->getAverageDiagSize(start, end);
diag_iter guilty = b->findClosestDiagBlock(end, b->diag_end(), av);
SchurDecompEig sd((QuasiTriangular&)*b); // works on b including diagonal structure
SchurDecompEig sd(*b); // works on b including diagonal structure
end = sd.bubbleEigen(guilty, end); // iterators are valid
++end;
q->multRight(sd.getQ());
......
......@@ -117,7 +117,7 @@ mxArray *
SylvParams::IntParamItem::createMatlabArray() const
{
mxArray *res = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
*((int *) mxGetData(res)) = value;
*reinterpret_cast<int *>(mxGetData(res)) = value;
return res;
}
......
......@@ -34,7 +34,7 @@ SymSchurDecomp::SymSchurDecomp(const ConstGeneralMatrix &mata)
double *w = lambda.base();
double *z = q.base();
lapack_int ldz = q.getLD();
std::vector<lapack_int> isuppz(2*std::max(1, (int) m));
std::vector<lapack_int> isuppz(2*std::max(1, static_cast<int>(m)));
double tmpwork;
lapack_int lwork = -1;
lapack_int tmpiwork;
......@@ -44,7 +44,7 @@ SymSchurDecomp::SymSchurDecomp(const ConstGeneralMatrix &mata)
// query for lwork and liwork
dsyevr("V", "A", "U", &n, a, &lda, vl, vu, il, iu, &abstol,
&m, w, z, &ldz, isuppz.data(), &tmpwork, &lwork, &tmpiwork, &liwork, &info);
lwork = (int) tmpwork;
lwork = static_cast<lapack_int>(tmpwork);
liwork = tmpiwork;
// allocate work arrays
std::vector<double> work(lwork);
......
......@@ -64,7 +64,7 @@ TriangularSylvester::solvi(double r, KronVector &d, double &eig_min) const
for (const_diag_iter di = matrixF->diag_begin();
di != matrixF->diag_end();
++di)
if ((*di).isReal())
if (di->isReal())
solviRealAndEliminate(r, di, d, eig_min);
else
solviComplexAndEliminate(r, di, d, eig_min);
......@@ -106,7 +106,7 @@ TriangularSylvester::solviip(double alpha, double betas,
const_diag_iter di = matrixF->diag_begin();
const_diag_iter dsi = matrixFF->diag_begin();
for (; di != matrixF->diag_end(); ++di, ++dsi)
if ((*di).isReal())
if (di->isReal())
solviipRealAndEliminate(alpha, betas, di, dsi, d, eig_min);
else
solviipComplexAndEliminate(alpha, betas, di, dsi, d, eig_min);
......@@ -118,14 +118,14 @@ TriangularSylvester::solviRealAndEliminate(double r, const_diag_iter di,
KronVector &d, double &eig_min) const
{
// di is real
int jbar = (*di).getIndex();
double f = *((*di).getAlpha());
int jbar = di->getIndex();
double f = *(di->getAlpha());
KronVector dj(d, jbar);
// solve system
if (std::abs(r*f) > diag_zero)
solvi(r*f, dj, eig_min);
// calculate y
KronVector y((const KronVector &)dj);
KronVector y(const_cast<const KronVector &>(dj));
KronUtils::multKron(*matrixF, *matrixK, y);
y.mult(r);
double divisor = 1.0;
......@@ -137,8 +137,7 @@ TriangularSylvester::solviEliminateReal(const_diag_iter di, KronVector &d,
const KronVector &y, double divisor) const
{
for (const_row_iter ri = matrixF->row_begin(*di);
ri != matrixF->row_end(*di);
++ri)
ri != matrixF->row_end(*di); ++ri)
{
KronVector dk(d, ri.getCol());
dk.add(-(*ri)/divisor, y);
......@@ -150,12 +149,12 @@ TriangularSylvester::solviComplexAndEliminate(double r, const_diag_iter di,
KronVector &d, double &eig_min) const
{
// di is complex
int jbar = (*di).getIndex();
int jbar = di->getIndex();
// pick data
double alpha = *(*di).getAlpha();
double beta1 = (*di).getBeta2();
double beta2 = -(*di).getBeta1();
double aspbs = (*di).getDeterminant();
double alpha = *di->getAlpha();
double beta1 = di->getBeta2();
double beta2 = -di->getBeta1();
double aspbs = di->getDeterminant();
KronVector dj(d, jbar);
KronVector djj(d, jbar+1);
// solve
......@@ -177,8 +176,7 @@ TriangularSylvester::solviEliminateComplex(const_diag_iter di, KronVector &d,
double divisor) const
{
for (const_row_iter ri = matrixF->row_begin(*di);
ri != matrixF->row_end(*di);
++ri)
ri != matrixF->row_end(*di); ++ri)
{
KronVector dk(d, ri.getCol());
dk.add(-ri.a()/divisor, y1);
......@@ -192,17 +190,17 @@ TriangularSylvester::solviipRealAndEliminate(double alpha, double betas,
KronVector &d, double &eig_min) const
{
// di, and dsi are real
int jbar = (*di).getIndex();
int jbar = di->getIndex();
double aspbs = alpha*alpha+betas;
// pick data
double f = *((*di).getAlpha());
double f = *(di->getAlpha());
double fs = f*f;
KronVector dj(d, jbar);
// solve
if (fs*aspbs > diag_zero_sq)
solviip(f*alpha, fs*betas, dj, eig_min);
KronVector y1((const KronVector &)dj);
KronVector y2((const KronVector &)dj);
KronVector y1(const_cast<const KronVector &>(dj));
KronVector y2(const_cast<const KronVector &>(dj));
KronUtils::multKron(*matrixF, *matrixK, y1);
y1.mult(2*alpha);
KronUtils::multKron(*matrixFF, *matrixKK, y2);
......@@ -234,32 +232,32 @@ TriangularSylvester::solviipComplexAndEliminate(double alpha, double betas,
KronVector &d, double &eig_min) const
{
// di, and dsi are complex
int jbar = (*di).getIndex();
int jbar = di->getIndex();
double aspbs = alpha*alpha+betas;
// pick data
double gamma = *((*di).getAlpha());
double delta1 = (*di).getBeta2(); // swap because of transpose
double delta2 = -(*di).getBeta1();
double gspds = (*di).getDeterminant();
double gamma = *(di->getAlpha());
double delta1 = di->getBeta2(); // swap because of transpose
double delta2 = -di->getBeta1();
double gspds = di->getDeterminant();
KronVector dj(d, jbar);