From a2f2814cb2f0cb48f717c2c31303966be68a9eb6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 24 Jan 2019 13:08:05 +0100
Subject: [PATCH] Dynare++: ensure that all calls to BLAS/LAPACK use the
correct LD values
Many BLAS/LAPACK calls were making the assumption that LD==rows when passing
matrices. In some cases this was correct (because of implementation details,
in particular because how the copy constructor of GeneralMatrix is implemented).
But in other cases it was a bug.
With this commit, the actual value for LD is systematically used (this fixes
existing bugs and prevent possible future ones if the implementation details
were changed).
---
dynare++/integ/cc/vector_function.cc | 4 ++--
dynare++/kord/decision_rule.cc | 4 ++--
dynare++/kord/first_order.cc | 9 ++++++---
dynare++/kord/korder.cc | 10 ++++++----
dynare++/sylv/cc/GeneralMatrix.cc | 8 ++++----
dynare++/sylv/cc/QuasiTriangular.cc | 11 ++++++-----
dynare++/sylv/cc/SchurDecomp.cc | 5 +++--
dynare++/sylv/cc/SchurDecompEig.cc | 4 ++--
dynare++/sylv/cc/SimilarityDecomp.cc | 3 ++-
dynare++/sylv/cc/SylvMatrix.cc | 12 ++++++------
10 files changed, 39 insertions(+), 31 deletions(-)
diff --git a/dynare++/integ/cc/vector_function.cc b/dynare++/integ/cc/vector_function.cc
index 92afd6ee76..d93881a152 100644
--- a/dynare++/integ/cc/vector_function.cc
+++ b/dynare++/integ/cc/vector_function.cc
@@ -120,12 +120,12 @@ GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix &vcov)
{
A = vcov;
- lapack_int rows = A.numRows();
+ lapack_int rows = A.numRows(), lda = A.getLD();
for (int i = 0; i < rows; i++)
for (int j = i+1; j < rows; j++)
A.get(i, j) = 0.0;
lapack_int info;
- dpotrf("L", &rows, A.base(), &rows, &info);
+ dpotrf("L", &rows, A.base(), &lda, &info);
// todo: raise if |info!=1|
}
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index e1998d0d77..f61d937b5a 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -563,13 +563,13 @@ void
RandomShockRealization::choleskyFactor(const TwoDMatrix &v)
{
factor = v;
- lapack_int rows = factor.nrows();
+ lapack_int rows = factor.nrows(), lda = factor.getLD();
for (int i = 0; i < rows; i++)
for (int j = i+1; j < rows; j++)
factor.get(i, j) = 0.0;
lapack_int info;
- dpotrf("L", &rows, factor.base(), &rows, &info);
+ dpotrf("L", &rows, factor.base(), &lda, &info);
KORD_RAISE_IF(info != 0,
"Info!=0 in RandomShockRealization::choleskyFactor");
}
diff --git a/dynare++/kord/first_order.cc b/dynare++/kord/first_order.cc
index 7986481337..733d51670d 100644
--- a/dynare++/kord/first_order.cc
+++ b/dynare++/kord/first_order.cc
@@ -145,6 +145,7 @@ FirstOrder::solve(const TwoDMatrix &fd)
matD.place(fyplus, 0, ypart.nys()+ypart.nstat);
for (int i = 0; i < ypart.nboth; i++)
matD.get(ypart.ny()+i, ypart.npred+i) = 1.0;
+ lapack_int ldb = matD.getLD();
// form matrix $E$
TwoDMatrix matE(n, n);
@@ -155,18 +156,20 @@ FirstOrder::solve(const TwoDMatrix &fd)
for (int i = 0; i < ypart.nboth; i++)
matE.get(ypart.ny()+i, ypart.nys()+ypart.nstat+i) = -1.0;
matE.mult(-1.0);
+ lapack_int lda = matE.getLD();
// solve generalized Schur
TwoDMatrix vsl(n, n);
TwoDMatrix vsr(n, n);
+ lapack_int ldvsl = vsl.getLD(), ldvsr = vsr.getLD();
lapack_int lwork = 100*n+16;
Vector work(lwork);
auto *bwork = new lapack_int[n];
lapack_int info;
lapack_int sdim2 = sdim;
- dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &n,
- matD.getData().base(), &n, &sdim2, alphar.base(), alphai.base(),
- beta.base(), vsl.getData().base(), &n, vsr.getData().base(), &n,
+ dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &lda,
+ matD.getData().base(), &ldb, &sdim2, alphar.base(), alphai.base(),
+ beta.base(), vsl.getData().base(), &ldvsl, vsr.getData().base(), &ldvsr,
work.base(), &lwork, bwork, &info);
if (info)
{
diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc
index c9bcef5be6..89893f7bcf 100644
--- a/dynare++/kord/korder.cc
+++ b/dynare++/kord/korder.cc
@@ -17,9 +17,9 @@ void
PLUMatrix::calcPLU()
{
lapack_int info;
- lapack_int rows = nrows();
+ lapack_int rows = nrows(), lda = ld;
inv = (const Vector &) getData();
- dgetrf(&rows, &rows, inv.base(), &rows, ipiv, &info);
+ dgetrf(&rows, &rows, inv.base(), &lda, ipiv, &info);
}
/* Here we just call the LAPACK machinery to multiply by the inverse. */
@@ -30,11 +30,13 @@ PLUMatrix::multInv(TwoDMatrix &m) const
KORD_RAISE_IF(m.nrows() != ncols(),
"The matrix is not square in PLUMatrix::multInv");
lapack_int info;
+ lapack_int lda = ld;
lapack_int mcols = m.ncols();
lapack_int mrows = m.nrows();
+ lapack_int ldb = m.getLD();
double *mbase = m.getData().base();
- dgetrs("N", &mrows, &mcols, inv.base(), &mrows, ipiv,
- mbase, &mrows, &info);
+ dgetrs("N", &mrows, &mcols, inv.base(), &lda, ipiv,
+ mbase, &ldb, &info);
KORD_RAISE_IF(info != 0,
"Info!=0 in PLUMatrix::multInv");
}
diff --git a/dynare++/sylv/cc/GeneralMatrix.cc b/dynare++/sylv/cc/GeneralMatrix.cc
index 7a659c9734..7ba83d783a 100644
--- a/dynare++/sylv/cc/GeneralMatrix.cc
+++ b/dynare++/sylv/cc/GeneralMatrix.cc
@@ -448,7 +448,7 @@ ConstGeneralMatrix::multVecTrans(double a, Vector &x, double b,
blas_int mm = rows;
blas_int nn = cols;
double alpha = b;
- blas_int lda = rows;
+ blas_int lda = ld;
blas_int incx = d.skip();
double beta = a;
blas_int incy = x.skip();
@@ -471,9 +471,9 @@ ConstGeneralMatrix::multInvLeft(const char *trans, int mrows, int mcols, int mld
GeneralMatrix inv(*this);
std::vector<lapack_int> ipiv(rows);
lapack_int info;
- lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld;
- dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv.data(), &info);
- dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv.data(), d,
+ lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld, lda = inv.ld;
+ dgetrf(&rows2, &rows2, inv.getData().base(), &lda, ipiv.data(), &info);
+ dgetrs(trans, &rows2, &mcols2, inv.base(), &lda, ipiv.data(), d,
&mld2, &info);
}
}
diff --git a/dynare++/sylv/cc/QuasiTriangular.cc b/dynare++/sylv/cc/QuasiTriangular.cc
index b0fda9f970..f3a20d2fbc 100644
--- a/dynare++/sylv/cc/QuasiTriangular.cc
+++ b/dynare++/sylv/cc/QuasiTriangular.cc
@@ -432,8 +432,9 @@ QuasiTriangular::QuasiTriangular(int p, const QuasiTriangular &t)
blas_int d_size = diagonal.getSize();
double alpha = 1.0;
double beta = 0.0;
+ blas_int lda = t.getLD(), ldb = t.getLD(), ldc = ld;
dgemm("N", "N", &d_size, &d_size, &d_size, &alpha, aux.base(),
- &d_size, t.getData().base(), &d_size, &beta, getData().base(), &d_size);
+ &lda, t.getData().base(), &ldb, &beta, getData().base(), &ldc);
}
QuasiTriangular::QuasiTriangular(const SchurDecomp &decomp)
@@ -590,7 +591,7 @@ QuasiTriangular::solvePre(Vector &x, double &eig_min)
}
blas_int nn = diagonal.getSize();
- blas_int lda = diagonal.getSize();
+ blas_int lda = ld;
blas_int incx = x.skip();
dtrsv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
}
@@ -616,7 +617,7 @@ QuasiTriangular::solvePreTrans(Vector &x, double &eig_min)
}
blas_int nn = diagonal.getSize();
- blas_int lda = diagonal.getSize();
+ blas_int lda = ld;
blas_int incx = x.skip();
dtrsv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
}
@@ -627,7 +628,7 @@ QuasiTriangular::multVec(Vector &x, const ConstVector &b) const
{
x = b;
blas_int nn = diagonal.getSize();
- blas_int lda = diagonal.getSize();
+ blas_int lda = ld;
blas_int incx = x.skip();
dtrmv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
for (const_diag_iter di = diag_begin(); di != diag_end(); ++di)
@@ -645,7 +646,7 @@ QuasiTriangular::multVecTrans(Vector &x, const ConstVector &b) const
{
x = b;
blas_int nn = diagonal.getSize();
- blas_int lda = diagonal.getSize();
+ blas_int lda = ld;
blas_int incx = x.skip();
dtrmv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
for (const_diag_iter di = diag_begin(); di != diag_end(); ++di)
diff --git a/dynare++/sylv/cc/SchurDecomp.cc b/dynare++/sylv/cc/SchurDecomp.cc
index 7a101d47d4..8cc1abc751 100644
--- a/dynare++/sylv/cc/SchurDecomp.cc
+++ b/dynare++/sylv/cc/SchurDecomp.cc
@@ -13,13 +13,14 @@ SchurDecomp::SchurDecomp(const SqSylvMatrix &m)
{
lapack_int rows = m.numRows();
SqSylvMatrix auxt(m);
+ lapack_int lda = auxt.getLD(), ldvs = q.getLD();
lapack_int sdim;
std::vector<double> wr(rows), wi(rows);
lapack_int lwork = 6*rows;
std::vector<double> work(lwork);
lapack_int info;
- dgees("V", "N", nullptr, &rows, auxt.base(), &rows, &sdim,
- wr.data(), wi.data(), q.base(), &rows,
+ dgees("V", "N", nullptr, &rows, auxt.base(), &lda, &sdim,
+ wr.data(), wi.data(), q.base(), &ldvs,
work.data(), &lwork, nullptr, &info);
t_storage = std::make_unique<QuasiTriangular>(auxt.getData(), rows);
t = t_storage.get();
diff --git a/dynare++/sylv/cc/SchurDecompEig.cc b/dynare++/sylv/cc/SchurDecompEig.cc
index 3c5955d326..5193ed7bf1 100644
--- a/dynare++/sylv/cc/SchurDecompEig.cc
+++ b/dynare++/sylv/cc/SchurDecompEig.cc
@@ -54,12 +54,12 @@ SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
itadd = it;
--itadd;
- lapack_int n = getDim();
+ lapack_int n = getDim(), ldt = getT().getLD(), ldq = getQ().getLD();
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(), &n, getQ().base(), &n, &ifst, &ilst, work.data(),
+ dtrexc("V", &n, getT().base(), &ldt, getQ().base(), &ldq, &ifst, &ilst, work.data(),
&info);
if (info < 0)
throw SYLV_MES_EXCEPTION("Wrong argument to dtrexc.");
diff --git a/dynare++/sylv/cc/SimilarityDecomp.cc b/dynare++/sylv/cc/SimilarityDecomp.cc
index 0da2ebf7e8..93e5807fec 100644
--- a/dynare++/sylv/cc/SimilarityDecomp.cc
+++ b/dynare++/sylv/cc/SimilarityDecomp.cc
@@ -51,9 +51,10 @@ SimilarityDecomp::solveX(diag_iter start, diag_iter end,
lapack_int isgn = -1;
lapack_int m = A.numRows();
lapack_int n = B.numRows();
+ lapack_int lda = A.getLD(), ldb = B.getLD();
double scale;
lapack_int info;
- dtrsyl("N", "N", &isgn, &m, &n, A.base(), &m, B.base(), &n,
+ dtrsyl("N", "N", &isgn, &m, &n, A.base(), &lda, B.base(), &ldb,
C.base(), &m, &scale, &info);
if (info < -1)
throw SYLV_MES_EXCEPTION("Wrong parameter to LAPACK dtrsyl.");
diff --git a/dynare++/sylv/cc/SylvMatrix.cc b/dynare++/sylv/cc/SylvMatrix.cc
index 6af0480c03..1e0b3241bb 100644
--- a/dynare++/sylv/cc/SylvMatrix.cc
+++ b/dynare++/sylv/cc/SylvMatrix.cc
@@ -233,27 +233,27 @@ SqSylvMatrix::multInvLeft2(GeneralMatrix &a, GeneralMatrix &b,
Vector inv(data);
std::vector<lapack_int> ipiv(rows);
lapack_int info;
- lapack_int rows2 = rows;
- dgetrf(&rows2, &rows2, inv.base(), &rows2, ipiv.data(), &info);
+ lapack_int rows2 = rows, lda = ld;
+ dgetrf(&rows2, &rows2, inv.base(), &lda, ipiv.data(), &info);
// solve a
lapack_int acols = a.numCols();
double *abase = a.base();
- dgetrs("N", &rows2, &acols, inv.base(), &rows2, ipiv.data(),
+ dgetrs("N", &rows2, &acols, inv.base(), &lda, ipiv.data(),
abase, &rows2, &info);
// solve b
lapack_int bcols = b.numCols();
double *bbase = b.base();
- dgetrs("N", &rows2, &bcols, inv.base(), &rows2, ipiv.data(),
+ dgetrs("N", &rows2, &bcols, inv.base(), &lda, ipiv.data(),
bbase, &rows2, &info);
// condition numbers
std::vector<double> work(4*rows);
std::vector<lapack_int> iwork(rows);
double norm1 = getNorm1();
- dgecon("1", &rows2, inv.base(), &rows2, &norm1, &rcond1,
+ dgecon("1", &rows2, inv.base(), &lda, &norm1, &rcond1,
work.data(), iwork.data(), &info);
double norminf = getNormInf();
- dgecon("I", &rows2, inv.base(), &rows2, &norminf, &rcondinf,
+ dgecon("I", &rows2, inv.base(), &lda, &norminf, &rcondinf,
work.data(), iwork.data(), &info);
}
--
GitLab