stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Estimation DLL: first version of matrix library + still unfinished decision rules computation

time-shift
Sébastien Villemot 2010-02-09 19:15:59 +01:00
parent cf8f793c55
commit fb97797405
23 changed files with 1670 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@ -122,6 +122,11 @@ ylwrap
# Windows
/windows/dynare-version.nsi
# Estimation DLL tests
/mex/sources/estimation/tests/test-dr
/mex/sources/estimation/libmat/tests/test-qr
/mex/sources/estimation/libmat/tests/test-gsd
# Misc
!/matlab/swz/c-code/Makefile
!/mex/sources/kalman/cc/Makefile

2
Makefile.am
View File

@ -2,7 +2,7 @@ SUBDIRS = preprocessor doc tests
if HAVE_BLAS
if HAVE_LAPACK
SUBDIRS += dynare++
SUBDIRS += dynare++ mex/sources/estimation/tests mex/sources/estimation/libmat/tests
endif
endif

2
configure.ac
View File

@ -171,6 +171,8 @@ AC_CONFIG_FILES([Makefile
dynare++/integ/testing/Makefile
dynare++/kord/Makefile
dynare++/src/Makefile
mex/sources/estimation/tests/Makefile
mex/sources/estimation/libmat/tests/Makefile
])
AC_ARG_ENABLE([matlab], AS_HELP_STRING([--disable-matlab], [disable compilation of MEX files for MATLAB]), [], [enable_matlab=yes])

12
mex/sources/dynblas.h
View File

@ -8,7 +8,7 @@
* and MATLAB_VERSION (for version 7.4, define it to 0x0704).
*
*
* Copyright (C) 2009 Dynare Team
* Copyright (C) 2009-2010 Dynare Team
*
* This file is part of Dynare.
*
@ -62,6 +62,12 @@ extern "C" {
CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
BLDOU c, CONST_BLINT ldc);
#define dsymm FORTRAN_WRAPPER(dsymm)
void dsymm(BLCHAR side, BLCHAR uplo, CONST_BLINT m, CONST_BLINT n,
CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
BLDOU c, CONST_BLINT ldc);
#define dgemv FORTRAN_WRAPPER(dgemv)
void dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha,
CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx,
@ -99,6 +105,10 @@ extern "C" {
double ddot(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU y,
CONST_BLINT incy);
#define dsyr FORTRAN_WRAPPER(dsyr)
void dsyr(BLCHAR uplo, CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU x,
CONST_BLINT incx, BLDOU a, CONST_BLINT lda);
#ifdef __cplusplus
} /* extern "C" */
#endif

11
mex/sources/dynlapack.h
View File

@ -8,7 +8,7 @@
* and MATLAB_VERSION (for version 7.4, define it to 0x0704).
*
*
* Copyright (C) 2009 Dynare Team
* Copyright (C) 2009-2010 Dynare Team
*
* This file is part of Dynare.
*
@ -110,6 +110,15 @@ extern "C" {
LAINT isuppz, LADOU work, CONST_LAINT lwork, LAINT iwork, CONST_LAINT liwork,
LAINT info);
#define dgeqrf FORTRAN_WRAPPER(dgeqrf)
void dgeqrf(CONST_LAINT m, CONST_LAINT n, LADOU a, CONST_LAINT lda,
LADOU tau, LADOU work, CONST_LAINT lwork, LAINT info);
#define dormqr FORTRAN_WRAPPER(dormqr)
void dormqr(LACHAR side, LACHAR trans, CONST_LAINT m, CONST_LAINT n, CONST_LAINT k,
CONST_LADOU a, CONST_LAINT lda, CONST_LADOU tau, LADOU c, CONST_LAINT ldc,
LADOU work, CONST_LAINT lwork, LAINT info);
#ifdef __cplusplus
} /* extern "C" */
#endif

195
mex/sources/estimation/DecisionRules.cc Normal file
View File

@ -0,0 +1,195 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cassert>
#include <algorithm>
#include "DecisionRules.hh"
DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
const std::vector<size_t> &zeta_fwrd_arg,
const std::vector<size_t> &zeta_back_arg,
const std::vector<size_t> &zeta_mixed_arg,
const std::vector<size_t> &zeta_static_arg,
double qz_criterium) :
n(n_arg), p(p_arg), zeta_fwrd(zeta_fwrd_arg), zeta_back(zeta_back_arg),
zeta_mixed(zeta_mixed_arg), zeta_static(zeta_static_arg),
n_fwrd(zeta_fwrd.size()), n_back(zeta_back.size()),
n_mixed(zeta_mixed.size()), n_static(zeta_static.size()),
n_back_mixed(n_back+n_mixed), n_fwrd_mixed(n_fwrd+n_mixed),
n_dynamic(n-n_static),
S(n, n_static),
A(n_back_mixed + n + n_fwrd_mixed),
D(n_fwrd + n_back + 2*n_mixed),
E(n_fwrd + n_back + 2*n_mixed),
Z(n_fwrd + n_back + 2*n_mixed),
QR(n, n_static),
GSD(n_fwrd + n_back + 2*n_mixed, qz_criterium),
LU1(n_fwrd_mixed),
LU2(n_back_mixed),
LU3(n_static),
g_y_back(n_back_mixed),
g_y_back_tmp(n_back_mixed),
g_y_static(n_static, n_back_mixed),
A0s(n_static),
A0d(n_static, n_dynamic),
g_y_dynamic(n_dynamic, n_back_mixed),
g_y_static_tmp(n_static, n_back_mixed)
{
assert(n == n_back + n_fwrd + n_mixed + n_static);
set_union(zeta_fwrd.begin(), zeta_fwrd.end(),
zeta_mixed.begin(), zeta_mixed.end(),
back_inserter(zeta_fwrd_mixed));
set_union(zeta_back.begin(), zeta_back.end(),
zeta_mixed.begin(), zeta_mixed.end(),
back_inserter(zeta_back_mixed));
set_union(zeta_back_mixed.begin(), zeta_back_mixed.end(),
zeta_fwrd.begin(), zeta_fwrd.end(),
back_inserter(zeta_dynamic));
// Compute beta_back and pi_back
for(size_t i = 0; i < n_back_mixed; i++)
if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_back_mixed[i])
== zeta_mixed.end())
beta_back.push_back(i);
else
pi_back.push_back(i);
// Compute beta_fwrd and pi_fwrd
for(size_t i = 0; i < n_fwrd_mixed; i++)
if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_fwrd_mixed[i])
== zeta_mixed.end())
beta_fwrd.push_back(i);
else
pi_fwrd.push_back(i);
}
void
DecisionRules::compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
{
assert(jacobian.getRows() == n
&& jacobian.getCols() == (n_back_mixed + n + n_fwrd_mixed + p));
assert(g_y.getRows() == n && g_y.getCols() == n);
assert(g_u.getRows() == n && g_u.getCols() == p);
// Construct S, perform QR decomposition and get A = Q*jacobian
for(size_t i = 0; i < n_static; i++)
mat::col_copy(jacobian, n_back_mixed+zeta_static[i], S, i);
A = MatrixConstView(jacobian, 0, 0, n, n_back_mixed + n + n_fwrd_mixed);
QR.computeAndMultByQ(S, "L", "N", A);
// Construct matrix D
for (size_t j = 0; j < n_back_mixed; j++)
mat::col_copy(A, n_back_mixed + zeta_back_mixed[j], n_static, n - n_static,
D, j, 0);
MatrixView(D, 0, n_back_mixed, n - n_static, n_fwrd_mixed) = MatrixView(A, n_static, n_back_mixed + n, n - n_static, n_fwrd_mixed);
MatrixView(D, n - n_static, 0, n_mixed, n_fwrd + n_back + 2*n_mixed).setAll(0.0);
for (size_t i = 0; i < n_mixed; i++)
D(n - n_static + i, beta_back[i]) = 1.0;
// Construct matrix E
MatrixView(E, 0, 0, n - n_static, n_back_mixed) = MatrixView(A, n_static, 0, n - n_static, n_back_mixed);
for (size_t j = 0; j < n_fwrd_mixed; j++)
mat::col_copy(A, n_back_mixed + zeta_fwrd_mixed[j], n_static, n - n_static,
E, n_back_mixed + j, 0);
mat::negate(MatrixView(E, 0, 0, n - n_static, n_fwrd + n_back + 2*n_mixed));
MatrixView(E, n - n_static, 0, n_mixed, n_fwrd + n_back + 2*n_mixed).setAll(0.0);
for (size_t i = 0; i < n_mixed; i++)
E(n - n_static + i, n_back_mixed + beta_fwrd[i]) = 1.0;
// Perform the generalized Schur
size_t sdim;
GSD.compute(D, E, Z, sdim);
if (n_back_mixed != sdim)
throw BlanchardKahnException(true, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim);
// Compute DR for forward variables w.r. to endogenous
MatrixView Z21(Z, n_back_mixed, 0, n_fwrd_mixed, n_back_mixed),
Z22(Z, n_back_mixed, n_back_mixed, n_fwrd_mixed, n_fwrd_mixed);
try
{
LU1.invMult("N", Z22, Z21);
}
catch (LUSolver::LUException &e)
{
throw BlanchardKahnException(false, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim);
}
mat::negate(Z21);
const MatrixView &g_y_fwrd = Z21;
for (size_t i = 0; i < n_fwrd_mixed; i++)
mat::row_copy(g_y_fwrd, i, g_y, zeta_fwrd_mixed[i]);
// Compute DR for backward variables w.r. to endogenous
MatrixView Z11(Z, 0, 0, n_back_mixed, n_back_mixed),
T11(D, 0, 0, n_back_mixed, n_back_mixed),
S11(E, 0, 0, n_back_mixed, n_back_mixed);
mat::set_identity(g_y_back);
mat::transpose(Z11);
g_y_back_tmp = Z11;
LU2.invMult("N", g_y_back_tmp, g_y_back);
g_y_back_tmp = g_y_back;
blas::gemm("N", "N", 1.0, S11, g_y_back_tmp, 0.0, g_y_back);
LU2.invMult("N", T11, g_y_back);
g_y_back_tmp = g_y_back;
blas::gemm("N", "N", 1.0, Z11, g_y_back_tmp, 0.0, g_y_back);
// TODO: avoid to copy mixed variables again, rather test it...
for (size_t i = 0; i < n_back_mixed; i++)
mat::row_copy(g_y_back, i, g_y, zeta_back_mixed[i]);
// Compute DR for static variables w.r. to endogenous
g_y_static = MatrixView(A, 0, 0, n_static, n_back_mixed);
for (size_t i = 0; i < n_dynamic; i++)
{
mat::row_copy(g_y, zeta_dynamic[i], g_y_dynamic, i);
mat::col_copy(A, n_back_mixed + zeta_dynamic[i], A0d, i);
}
blas::gemm("N", "N", 1.0, A0d, g_y_dynamic, 1.0, g_y_static);
blas::gemm("N", "N", 1.0, g_y_fwrd, g_y_back, 0.0, g_y_static_tmp);
blas::gemm("N", "N", 1.0, MatrixView(A, 0, n_back_mixed + n, n_static, n_fwrd_mixed),
g_y_static_tmp, 1.0, g_y_static);
for (size_t i = 0; i < n_static; i++)
mat::col_copy(A, n_back_mixed + zeta_static[i], A0s, i);
LU3.invMult("N", A0s, g_y_static);
mat::negate(g_y_static);
for (size_t i = 0; i < n_static; i++)
mat::row_copy(g_y_static, i, g_y, zeta_static[i]);
// Compute DR for all endogenous w.r. to shocks
}
std::ostream &
operator<<(std::ostream &out, const DecisionRules::BlanchardKahnException &e)
{
if (e.order)
out << "The Blanchard-Kahn order condition is not satisfied: you have " << e.n_fwrd_vars << " forward variables for " << e.n_explosive_eigenvals << " explosive eigenvalues";
else
out << "The Blanchard Kahn rank condition is not satisfied";
out << std::endl;
return out;
}

64
mex/sources/estimation/DecisionRules.hh Normal file
View File

@ -0,0 +1,64 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstdlib>
#include <vector>
#include "Vector.hh"
#include "Matrix.hh"
#include "QRDecomposition.hh"
#include "GeneralizedSchurDecomposition.hh"
#include "LUSolver.hh"
class DecisionRules
{
private:
const size_t n, p;
const std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
const size_t n_fwrd, n_back, n_mixed, n_static, n_back_mixed, n_fwrd_mixed, n_dynamic;
std::vector<size_t> zeta_fwrd_mixed, zeta_back_mixed, zeta_dynamic,
beta_back, beta_fwrd, pi_back, pi_fwrd;
Matrix S, A, D, E, Z;
QRDecomposition QR;
GeneralizedSchurDecomposition GSD;
LUSolver LU1, LU2, LU3;
Matrix g_y_back, g_y_back_tmp;
Matrix g_y_static, A0s, A0d, g_y_dynamic, g_y_static_tmp;
public:
class BlanchardKahnException
{
public:
//! True if the model fails the order condition. False if it fails the rank condition.
const bool order;
const int n_fwrd_vars, n_explosive_eigenvals;
BlanchardKahnException(bool order_arg, int n_fwrd_vars_arg, int n_explosive_eigenvals_arg) : order(order_arg), n_fwrd_vars(n_fwrd_vars_arg), n_explosive_eigenvals(n_explosive_eigenvals_arg) {};
};
/*!
The zetas are supposed to follow C convention (first vector index is zero).
*/
DecisionRules(size_t n_arg, size_t p_arg, const std::vector<size_t> &zeta_fwrd_arg,
const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
const std::vector<size_t> &zeta_static_arg, double qz_criterium);
/*!
\param jacobian First columns are backetermined vars at t-1 (in the order of zeta_back_mixed), then all vars at t (in the orig order), then forward vars at t+1 (in the order of zeta_fwrd_mixed), then exogenous vars.
*/
void compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
};
std::ostream &operator<<(std::ostream &out, const DecisionRules::BlanchardKahnException &e);

80
mex/sources/estimation/libmat/BlasBindings.hh Normal file
View File

@ -0,0 +1,80 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _BLAS_BINDINGS_HH
#define _BLAS_BINDINGS_HH
#include <dynblas.h>
#include "Vector.hh"
#include "Matrix.hh"
namespace blas
{
/* Level 2 */
//! Symmetric rank 1 operation: A = alpha*X*X' + A
template<class Mat, class Vec>
inline void
syr(const char *uplo, double alpha, Vec X, Mat A)
{
assert(X.getSize() == A.getCols() && A.getCols() == A.getRows());
blas_int n = X.getSize();
blas_int incx = X.getStride();
blas_int lda = A.getLd();
dsyr(uplo, &n, &alpha, X.getData(), &incx, A.getData(), &lda);
}
/* Level 3 */
//! General matrix multiplication
template<class Mat1, class Mat2, class Mat3>
inline void
gemm(const char *transa, const char *transb,
double alpha, const Mat1 &A, const Mat2 &B,
double beta, Mat3 &C)
{
assert(A.getRows() == C.getRows());
assert(A.getCols() == B.getRows());
assert(B.getCols() == C.getCols());
blas_int m = A.getRows(), n = B.getCols(), k = A.getCols();
blas_int lda = A.getLd(), ldb = B.getLd(), ldc = C.getLd();
dgemm(transa, transb, &m, &n, &k, &alpha, A.getData(), &lda,
B.getData(), &ldb, &beta, C.getData(), &ldc);
}
//! Symmetric matrix multiplication
template<class Mat1, class Mat2, class Mat3>
inline void
symm(const char *side, const char *uplo,
double alpha, const Mat1 &A, const Mat2 &B,
double beta, Mat3 &C)
{
assert(A.getRows() == A.getCols());
assert(A.getRows() == C.getRows());
assert(A.getCols() == B.getRows());
assert(B.getCols() == C.getCols());
blas_int m = A.getRows(), n = B.getCols();
blas_int lda = A.getLd(), ldb = B.getLd(), ldc = C.getLd();
dsymm(side, uplo, &m, &n, &alpha, A.getData(), &lda,
B.getData(), &ldb, &beta, C.getData(), &ldc);
}
} // End of namespace
#endif

73
mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc Normal file
View File

@ -0,0 +1,73 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "GeneralizedSchurDecomposition.hh"
#include <cassert>
#include <cstdlib>
double GeneralizedSchurDecomposition::criterium_static;
GeneralizedSchurDecomposition::GeneralizedSchurDecomposition(size_t n_arg, double criterium_arg)
: n(n_arg), criterium(criterium_arg)
{
alphar = new double[n];
alphai = new double[n];
beta = new double[n];
lwork = 16*n+16; // Same heuristic choice than mjdgges
work = new double[(int) lwork];
vsl = new double[n*n];
bwork = new lapack_int[n];
}
GeneralizedSchurDecomposition::~GeneralizedSchurDecomposition()
{
delete[] alphar;
delete[] alphai;
delete[] beta;
delete[] work;
delete[] vsl;
delete[] bwork;
}
lapack_int
GeneralizedSchurDecomposition::selctg(const double *alphar, const double *alphai, const double *beta)
{
return((*alphar * *alphar + *alphai * *alphai) < criterium_static * *beta * *beta);
}
std::ostream &
operator<<(std::ostream &out, const GeneralizedSchurDecomposition::GSDException &e)
{
out << "DGGES return code " << e.info << ": ";
if (e.info < 0)
out << "argument " << -e.info << " has an illegal value";
else if (e.info <= e.n)
out << "the QZ iteration failed";
else if (e.info == e.n + 1)
out << "other than QZ iteration failed in DHGEQZ";
else if (e.info == e.n + 2)
out << "after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Generalized Schur form no longer satisfy SELCTG=TRUE. This could also be caused due to scaling";
else if (e.info == e.n + 3)
out << "reordering failed in DTGSEN";
out << std::endl;
return out;
}

110
mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh Normal file
View File

@ -0,0 +1,110 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <dynlapack.h>
#include "Vector.hh"
#include "Matrix.hh"
class GeneralizedSchurDecomposition
{
private:
const size_t n;
const double criterium;
lapack_int lwork;
double *alphar, *alphai, *beta, *vsl, *work;
lapack_int *bwork;
static double criterium_static;
static lapack_int selctg(const double *alphar, const double *alphai, const double *beta);
public:
class GSDException
{
public:
const lapack_int info, n;
GSDException(lapack_int info_arg, lapack_int n_arg) : info(info_arg), n(n_arg) {};
};
//! \todo Replace heuristic choice for workspace size by a query to determine the optimal size
GeneralizedSchurDecomposition(size_t n_arg, double criterium_arg);
~GeneralizedSchurDecomposition();
//! \todo Add a lock around the modification of criterium_static for making it thread-safe
template<class Mat1, class Mat2, class Mat3>
void compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException);
template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5>
/*!
\param[out] sdim Number of non-explosive generalized eigenvalues
*/
void compute(const Mat1 &D, const Mat2 &E, Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException);
template<class Vec1, class Vec2>
void getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx);
};
std::ostream &operator<<(std::ostream &out, const GeneralizedSchurDecomposition::GSDException &e);
template<class Mat1, class Mat2, class Mat3>
void
GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException)
{
assert(S.getRows() == n && S.getCols() == n
&& T.getRows() == n && T.getCols() == n
&& Z.getRows() == n && Z.getCols() == n);
lapack_int n2 = n;
lapack_int info, sdim2;
lapack_int lds = S.getLd(), ldt = T.getLd(), ldz = Z.getLd();
criterium_static = criterium;
// Here we are forced to give space for left Schur vectors, even if we don't use them, because of a bug in dgges()
dgges("N", "V", "S", &selctg, &n2, S.getData(), &lds, T.getData(), &ldt,
&sdim2, alphar, alphai, beta, vsl, &n2, Z.getData(), &ldz,
work, &lwork, bwork, &info);
if (info != 0)
throw GSDException(info, n2);
sdim = sdim2;
}
template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5>
void
GeneralizedSchurDecomposition::compute(const Mat1 &D, const Mat2 &E,
Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException)
{
assert(D.getRows() == n && D.getCols() == n
&& E.getRows() == n && E.getCols() == n);
S = D;
T = E;
compute(S, T, Z, sdim);
}
template<class Vec1, class Vec2>
void
GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx)
{
assert(eig_real.getSize() == n && eig_cmplx.getSize() == n);
double *par = alphar, *pai = alphai, *pb = beta,
*per = eig_real.getData(), *pei = eig_cmplx.getData();
while (par < alphar + n)
{
*per = *par++ / *pb;
*pei = *pai++ / *pb++;
per += eig_real.getStride();
pei += eig_cmplx.getStride();
}
}

30
mex/sources/estimation/libmat/LUSolver.cc Normal file
View File

@ -0,0 +1,30 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "LUSolver.hh"
LUSolver::LUSolver(size_t dim_arg) : dim(dim_arg)
{
ipiv = new lapack_int[dim];
}
LUSolver::~LUSolver()
{
delete[] ipiv;
}

63
mex/sources/estimation/libmat/LUSolver.hh Normal file
View File

@ -0,0 +1,63 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstdlib>
#include <cassert>
#include <dynlapack.h>
class LUSolver
{
private:
const size_t dim;
lapack_int *ipiv;
public:
class LUException
{
public:
const lapack_int info;
LUException(lapack_int info_arg) : info(info_arg) {};
};
LUSolver(size_t dim_arg);
~LUSolver();
/*!
Computes A^(-1)*B (possibly transposing A).
The output is stored in B.
The input matrix A is modified.
*/
template<class Mat1, class Mat2>
void invMult(const char *trans, Mat1 &A, Mat2 &B) throw (LUException);
};
template<class Mat1, class Mat2>
void
LUSolver::invMult(const char *trans, Mat1 &A, Mat2 &B) throw (LUException)
{
assert(A.getRows() == dim && A.getCols() == dim);
assert(B.getRows() == dim);
lapack_int n = dim, lda = A.getLd(), info;
dgetrf(&n, &n, A.getData(), &lda, ipiv, &info);
if (info != 0)
throw LUException(info);
lapack_int nrhs = B.getCols(), ldb = B.getLd();
dgetrs(trans, &n, &nrhs, A.getData(), &lda, ipiv, B.getData(), &ldb, &info);
assert(info == 0);
}

112
mex/sources/estimation/libmat/Matrix.cc Normal file
View File

@ -0,0 +1,112 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "Matrix.hh"
#include <cstring>
#include <cassert>
Matrix::Matrix(size_t rows_arg, size_t cols_arg) : rows(rows_arg), cols(cols_arg)
{
data = new double[rows*cols];
}
Matrix::Matrix(size_t size_arg) : rows(size_arg), cols(size_arg)
{
data = new double[rows*cols];
}
Matrix::Matrix(const Matrix &arg) : rows(arg.rows), cols(arg.cols)
{
data = new double[rows*cols];
memcpy(data, arg.data, rows*cols*sizeof(double));
}
Matrix::~Matrix()
{
delete[] data;
}
Matrix &
Matrix::operator= (const Matrix &arg)
{
assert(rows == arg.rows && cols == arg.cols);
memcpy(data, arg.data, rows*cols*sizeof(double));
return *this;
}
std::ostream &
operator<<(std::ostream &out, const Matrix &M)
{
mat::print(out, M);
return out;
}
MatrixView::MatrixView(double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg)
: data(data_arg), rows(rows_arg), cols(cols_arg), ld(ld_arg)
{
}
MatrixView::MatrixView(Matrix &arg, size_t row_offset, size_t col_offset,
size_t rows_arg, size_t cols_arg) :
data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
{
assert(row_offset < arg.getRows()
&& row_offset + rows_arg <= arg.getRows()
&& col_offset < arg.getCols()
&& col_offset + cols_arg <= arg.getCols());
}
std::ostream &
operator<<(std::ostream &out, const MatrixView &M)
{
mat::print(out, M);
return out;
}
MatrixView &
MatrixView::operator= (const MatrixView &arg)
{
assert(rows == arg.getRows() && cols == arg.getCols());
for (size_t j = 0; j < cols; j++)
memcpy(data + j*ld, arg.getData() + j*arg.getLd(), rows*sizeof(double));
return *this;
}
MatrixConstView::MatrixConstView(const double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg)
: data(data_arg), rows(rows_arg), cols(cols_arg), ld(ld_arg)
{
}
MatrixConstView::MatrixConstView(const Matrix &arg, size_t row_offset, size_t col_offset,
size_t rows_arg, size_t cols_arg) :
data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
{
assert(row_offset < arg.getRows()
&& row_offset + rows_arg <= arg.getRows()
&& col_offset < arg.getCols()
&& col_offset + cols_arg <= arg.getCols());
}
std::ostream &
operator<<(std::ostream &out, const MatrixConstView &M)
{
mat::print(out, M);
return out;
}

348
mex/sources/estimation/libmat/Matrix.hh Normal file
View File

@ -0,0 +1,348 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _MATRIX_HH
#define _MATRIX_HH
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <cassert>
#include <cstring>
#include "Vector.hh"
/*
This header defines three matrix classes, which implement a "matrix concept"
(much like the concepts of the Standard Template Library or of Boost
Library). The first class is a matrix owning the data space for its
elements, and the other two are matrix "views" of another matrix, i.e. a
contiguous submatrix. This design philosophy is close to the design of the
GNU Scientific Library, but here using the syntactic power of C++ class and
templates, while achieving very high efficiency.
These classes can be used with various templated functions, including
wrappers around the BLAS primitives.
The expressions required to be valid for a class M implementing the "matrix concept" are:
- M.getRows(): return number of rows
- M.getCols(): return number of columns
- M.getLd(): return the leading dimension (here the offset between two columns in the data space, since we are in column-major order)
- M.getData(): return the pointer to the data space
- M(i,j): get an element of the matrix
The expressions required to be valid for a class M implementing the "mutable matrix concept" are (in addition to those of "matrix concept"):
- M = X: assignment operator
- M(i,j) = d: assign an element of the matrix
- M.setAll(d): set all the elements of the matrix
*/
//! A full matrix, implements the "mutable matrix concept"
/*! Owns the data space for the elements */
class Matrix
{
private:
//! Stored in column-major order, as in Fortran and MATLAB
double *data;
const size_t rows, cols;
public:
Matrix(size_t rows_arg, size_t cols_arg);
Matrix(size_t size_arg);
Matrix(const Matrix &arg);
~Matrix();
inline size_t getRows() const { return rows; }
inline size_t getCols() const { return cols; }
inline size_t getLd() const { return rows; }
inline double *getData() { return data; }
inline const double *getData() const { return data; }
inline void setAll(double val) { std::fill_n(data, rows*cols, val); }
inline double &operator() (size_t i, size_t j) { return data[i+j*rows]; }
inline const double &operator() (size_t i, size_t j) const { return data[i+j*rows]; }
//! Assignment operator, only works for matrices of same dimension
template<class Mat>
Matrix &
operator= (const Mat &arg)
{
assert(rows == arg.getRows() && cols == arg.getCols());
for (size_t j = 0; j < cols; j++)
memcpy(data + j*rows, arg.getData() + j*arg.getLd(), rows*sizeof(double));
return *this;
}
//! The copy assignment operator, which is not generated by the template assignment operator
/*! See C++ standard, §12.8.9, in the footnote */
Matrix &operator= (const Matrix &arg);
};
//! A contiguous submatrix of another matrix, implements the "mutable matrix concept"
/*! Does not own the data space for the elements, so depends on another matrix */
class MatrixView
{
private:
double *const data;
const size_t rows, cols, ld;
public:
MatrixView(double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg);
MatrixView(Matrix &arg, size_t row_offset, size_t col_offset,
size_t rows_arg, size_t cols_arg);
inline size_t getRows() const { return rows; }
inline size_t getCols() const { return cols; }
inline size_t getLd() const { return ld; }
inline double *getData() { return data; }
inline const double *getData() const { return data; }
inline void setAll(double val)
{
for (double *p = data; p < data + rows*ld; p += ld)
std::fill_n(p, cols, val);
}
inline double &operator() (size_t i, size_t j) { return data[i+j*ld]; }
inline const double &operator() (size_t i, size_t j) const { return data[i+j*ld]; }
//! Assignment operator, only works for matrices of same dimension
template<class Mat>
MatrixView &
operator= (const Mat &arg)
{
assert(rows == arg.getRows() && cols == arg.getCols());
for (size_t j = 0; j < cols; j++)
memcpy(data + j*ld, arg.getData() + j*arg.getLd(), rows*sizeof(double));
return *this;
}
//! The copy assignment operator, which is not generated by the template assignment operator
/*! See C++ standard, §12.8.9, in the footnote */
MatrixView &operator= (const MatrixView &arg);
};
//! Like MatrixView, but cannot be modified (implements the "matrix concept")
class MatrixConstView
{
private:
const double *const data;
const size_t rows, cols, ld;
public:
MatrixConstView(const double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg);
MatrixConstView(const Matrix &arg, size_t row_offset, size_t col_offset,
size_t rows_arg, size_t cols_arg);
inline size_t getRows() const { return rows; }
inline size_t getCols() const { return cols; }
inline size_t getLd() const { return ld; }
inline const double *getData() const { return data; }
inline const double &operator() (size_t i, size_t j) const { return data[i+j*ld]; }
};
std::ostream &operator<<(std::ostream &out, const Matrix &M);
std::ostream &operator<<(std::ostream &out, const MatrixView &M);
std::ostream &operator<<(std::ostream &out, const MatrixConstView &M);
namespace mat
{
template<class Mat>
void
print(std::ostream &out, const Mat &M)
{
for(size_t i = 0; i < M.getRows(); i++)
{
for(size_t j = 0; j < M.getCols(); j++)
out << std::setw(13) << std::right << M(i, j) << " ";
out << std::endl;
}
}
template<class Mat>
inline VectorView
get_col(Mat &M, size_t j)
{
return VectorView(M.getData()+j*M.getLd(), M.getRows(), 1);
}
template<class Mat>
inline VectorView
get_row(Mat &M, size_t i)
{
return VectorView(M.getData()+i, M.getCols(), M.getLd());
}
template<class Mat>
inline VectorConstView
get_col(const Mat &M, size_t j)
{
return VectorView(M.getData()+j*M.getLd(), M.getRows(), 1);
}
template<class Mat>
inline VectorConstView
get_row(const Mat &M, size_t i)
{
return VectorView(M.getData()+i, M.getCols(), M.getLd());
}
template<class Mat1, class Mat2>
inline void
col_copy(const Mat1 &src, size_t col_src, Mat2 &dest, size_t col_dest)
{
assert(src.getRows() == dest.getRows()
&& col_src < src.getCols() && col_dest < dest.getCols());
memcpy(dest.getData() + col_dest*dest.getLd(),
const_cast<double *>(src.getData()) + col_src*src.getLd(),
src.getRows()*sizeof(double));
}
template<class Mat1, class Mat2>
inline void
col_copy(const Mat1 &src, size_t col_src, size_t row_offset_src, size_t row_nb,
Mat2 &dest, size_t col_dest, size_t row_offset_dest)
{
assert(src.getRows() == dest.getRows()
&& col_src < src.getCols() && col_dest < dest.getCols()
&& row_offset_src < src.getRows() && row_offset_src+row_nb <= src.getRows()
&& row_offset_dest < dest.getRows() && row_offset_dest+row_nb <= dest.getRows());
memcpy(dest.getData() + row_offset_dest + col_dest*dest.getLd(),
src.getData() + row_offset_src + col_src*src.getLd(),
row_nb*sizeof(double));
}
template<class Mat1, class Mat2>
inline void
row_copy(const Mat1 &src, size_t row_src, Mat2 &dest, size_t row_dest)
{
assert(src.getCols() == dest.getCols()
&& row_src < src.getRows() && row_dest < dest.getRows());
const double *p1 = src.getData() + row_src;
double *p2 = dest.getData() + row_dest;
while (p1 < src.getData() + src.getRows() * src.getLd())
{
*p2 = *p1;
p1 += src.getLd();
p2 += dest.getLd();
}
}
template<class Mat>
inline void
col_set(Mat &M, size_t col, size_t row_offset, size_t row_nb, double val)
{
assert(col < M.getCols());
assert(row_offset < M.getRows() && row_offset + row_nb <= M.getRows());
std::fill_n(M.getData() + M.getLd()*col + row_offset, row_nb, val);
}
//! Copy under the diagonal the elements above the diagonal
template<class Mat>
inline void
copy_upper_to_lower(Mat &M)
{
size_t d = std::min(M.getCols(), M.getRows());
for (size_t i = 0; i < d; i++)
for (size_t j = 0; j < i; j++)
M(i, j) = M(j, i);
}
//! Copy above the diagonal the elements under the diagonal
template<class Mat>
inline void
copy_lower_to_upper(Mat &M)
{
size_t d = std::min(M.getCols(), M.getRows());
for (size_t i = 0; i < d; i++)
for (size_t j = 0; j < i; j++)
M(j, i) = M(i, j);
}
//! Fill the matrix with the identity matrix
template<class Mat>
inline void
set_identity(Mat &M)
{
M.setAll(0.0);
size_t d = std::min(M.getCols(), M.getRows());
for (size_t i = 0; i < d; i++)
M(i, i) = 1.0;
}
//! In-place transpose of a square matrix
template<class Mat>
inline void
transpose(Mat &M)
{
assert(M.getRows() == M.getCols());
for (size_t i = 0; i < M.getRows(); i++)
for (size_t j = 0; j < i; j++)
std::swap(M(i,j), M(j,i));
}
//! Computes m1 = m1 + m2
template<class Mat1, class Mat2>
void
add(Mat1 m1, const Mat2 m2)
{
assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
double *p1 = m1.getData();
const double *p2 = m2.getData();
while (p1 < m1.getData() + m1.getRows() * m1.getLd())
{
double *pp1, *pp2;
while (pp1 < p1 + m1.getCols())
*pp1++ += *pp2++;
p1 += m1.getLd();
p2 += m2.getLd();
}
}
//! Computes m1 = m1 m2
template<class Mat1, class Mat2>
void
sub(Mat1 m1, const Mat2 m2)
{
assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
double *p1 = m1.getData();
const double *p2 = m2.getData();
while (p1 < m1.getData() + m1.getRows() * m1.getLd())
{
double *pp1, *pp2;
while (pp1 < p1 + m1.getCols())
*pp1++ -= *pp2++;
p1 += m1.getLd();
p2 += m2.getLd();
}
}
//! Does m = -m
template<class Mat>
void
negate(Mat m)
{
double *p = m.getData();
while (p < m.getData() + m.getRows() * m.getLd())
{
double *pp;
while (pp < p + m.getCols())
{
*pp = -*pp;
pp++;
}
p += m.getLd();
}
}
} // End of namespace
#endif

36
mex/sources/estimation/libmat/QRDecomposition.cc Normal file
View File

@ -0,0 +1,36 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <algorithm> // For std::min
#include "QRDecomposition.hh"
QRDecomposition::QRDecomposition(size_t rows_arg, size_t cols_arg) :
rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), lwork(rows*cols),
H(rows), Q2(rows), v(rows)
{
work = new double[lwork];
tau = new double[mind];
}
QRDecomposition::~QRDecomposition()
{
delete[] work;
delete[] tau;
}

70
mex/sources/estimation/libmat/QRDecomposition.hh Normal file
View File

@ -0,0 +1,70 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <algorithm> // For std::min()
#include <dynlapack.h>
#include "Vector.hh"
#include "Matrix.hh"
#include "BlasBindings.hh"
class QRDecomposition
{
private:
const size_t rows, cols, mind;
lapack_int lwork;
double *work, *tau;
Matrix H, Q2;
Vector v;
public:
/*!
\todo Replace heuristic choice for workspace size by a query to determine the optimal size
*/
QRDecomposition(size_t rows_arg, size_t cols_arg);
~QRDecomposition();
//! Performs the QR decomposition of a matrix, and multiplies another matrix by Q
/*!
\param[in,out] A On input, the matrix to be decomposed. On output, equals to the output of dgeqrf
\param[in] side Specifies on which side is Q in the multiplication, either "L" or "R"
\param[in] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N"
\param[in,out] C The matrix to be multiplied by Q, modified in place
*/
template<class Mat1, class Mat2>
void computeAndMultByQ(Mat1 &A, const char *side, const char *trans, Mat2 &C);
};
template<class Mat1, class Mat2>
void
QRDecomposition::computeAndMultByQ(Mat1 &A, const char *side, const char *trans, Mat2 &C)
{
assert(A.getRows() == rows && A.getCols() == cols);
assert(C.getRows() == rows);
lapack_int m = rows, n = cols, lda = A.getLd();
lapack_int info;
dgeqrf(&m, &n, A.getData(), &lda, tau, work, &lwork, &info);
assert(info == 0);
n = C.getCols();
lapack_int k = mind, ldc = C.getLd();
dormqr(side, trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
work, &lwork, &info);
assert(info == 0);
}

98
mex/sources/estimation/libmat/Vector.cc Normal file
View File

@ -0,0 +1,98 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "Vector.hh"
#include <cstring>
#include <cassert>
Vector::Vector(size_t size_arg) : size(size_arg)
{
data = new double[size];
}
Vector::Vector(const Vector &arg) : size(arg.size)
{
memcpy(data, arg.data, size*sizeof(double));
}
Vector::~Vector()
{
delete[] data;
}
Vector &
Vector::operator= (const Vector &arg)
{
assert(size == arg.size);
memcpy(data, arg.data, size*sizeof(double));
return *this;
}
std::ostream &
operator<<(std::ostream &out, const Vector &V)
{
vec::print(out, V);
return out;
}
VectorView::VectorView(Vector &arg, size_t offset, size_t size_arg)
: data(arg.getData() + offset), size(size_arg), stride(1)
{
}
VectorView::VectorView(double *data_arg, size_t size_arg, size_t stride_arg)
: data(data_arg), size(size_arg), stride(stride_arg)
{
}
VectorView &
VectorView::operator= (const VectorView &arg)
{
assert(size == arg.getSize());
double *p1;
const double *p2;
for (p1 = data, p2 = arg.getData(); p1 < data + size*stride; p1 += stride, p2 += arg.getStride())
*p1 = *p2;
return *this;
}
std::ostream &
operator<<(std::ostream &out, const VectorView &V)
{
vec::print(out, V);
return out;
}
VectorConstView::VectorConstView(const Vector &arg, size_t offset, size_t size_arg)
: data(arg.getData() + offset), size(size_arg), stride(1)
{
}
VectorConstView::VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg)
: data(data_arg), size(size_arg), stride(stride_arg)
{
}
std::ostream &
operator<<(std::ostream &out, const VectorConstView &V)
{
vec::print(out, V);
return out;
}

206
mex/sources/estimation/libmat/Vector.hh Normal file
View File

@ -0,0 +1,206 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _VECTOR_HH
#define _VECTOR_HH
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cassert>
/*
This header defines three vector classes, which implement a "vector concept"
(much like the concepts of the Standard Template Library or of Boost
Library). The first class is a vector owning the data space for its elements,
and the other two are vector "views" of something else (either a subvector of
another vector, or a row or column of a matrix). This design philosophy is
close to the design of the GNU Scientific Library, but here using the
syntactic power of C++ class and templates, while achieving very high
efficiency.
These classes can be used with various templated functions, including
wrappers around the BLAS primitives.
The expressions required to be valid for a class V implementing the "vector concept" are:
- V.getSize(): return number of elements
- V.getStride(): return the stride, i.e. the offset between two elements in the data space
- V.getData(): return the pointer to the data space
- V(i): get an element of the vector
The expressions required to be valid for a class V implementing the "mutable vector concept" are (in addition to those of "vector concept"):
- V = X: assignment operator
- V(i) = d: assign an element of the vector
- V.setAll(d): set all the elements of the vector
*/
//! A full vector, implements the "mutable vector concept"
/*! Owns the data space for the elements */
class Vector
{
private:
double *data;
const size_t size;
public:
Vector(size_t size_arg);
Vector(const Vector &arg);
~Vector();
inline size_t getSize() const { return size; }
inline size_t getStride() const { return 1; }
inline double *getData() { return data; }
inline const double *getData() const { return data; }
inline void setAll(double val) { std::fill_n(data, size, val); }
inline double &operator() (size_t i) { return data[i]; }
inline const double &operator() (size_t i) const { return data[i]; }
//! Assignment operator, only works for vectors of same size
template<class Vec>
Vector &operator= (const Vec &arg)
{
assert(size == arg.getSize());
const double *p2 = arg.getData();
for(double *p1 = data; p1 < data + size; p1++)
{
*p1 = *p2;
p2 += arg.getStride();
}
return *this;
}
//! The copy assignment operator, which is not generated by the template assignment operator
/*! See C++ standard, §12.8.9, in the footnote */
Vector &operator= (const Vector &arg);
};
//! A vector view (i.e. a subvector of another vector, or a row or column of a matrix), implements the "mutable vector concept"
/*! Does not own the data space for the elements, so depends on another vector or matrix */
class VectorView
{
private:
double *const data;
const size_t size, stride;
public:
VectorView(Vector &arg, size_t offset, size_t size_arg);
VectorView(double *data_arg, size_t size_arg, size_t stride_arg);
inline size_t getSize() const { return size; }
inline size_t getStride() const { return stride; }
inline double *getData() { return data; }
inline const double *getData() const { return data; }
inline void setAll(double val)
{
for (double *p = data; p < data + size*stride; p += stride)
*p = val;
}
inline double &operator() (size_t i) { return data[i*stride]; }
inline const double &operator() (size_t i) const { return data[i*stride]; }
//! Assignment operator, only works for vectors of same size
template<class Vec>
VectorView &
operator= (const Vec &arg)
{
assert(size == arg.getSize());
double *p1;
const double *p2;
for (p1 = data, p2 = arg.getData(); p1 < data + size*stride; p1 += stride, p2 += arg.getStride())
*p1 = *p2;
return *this;
}
//! The copy assignment operator, which is not generated by the template assignment operator
/*! See C++ standard, §12.8.9, in the footnote */
VectorView &operator= (const VectorView &arg);
};
//! Like a VectorView, but the data cannot be modified (implements the "vector concept")
class VectorConstView
{
private:
const double *const data;
const size_t size, stride;
public:
VectorConstView(const Vector &arg, size_t offset, size_t size_arg);
VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg);
inline size_t getSize() const { return size; }
inline size_t getStride() const { return stride; }
inline const double *getData() const { return data; }
inline const double &operator() (size_t i) const { return data[i*stride]; }
};
std::ostream &operator<<(std::ostream &out, const Vector &V);
std::ostream &operator<<(std::ostream &out, const VectorView &V);
std::ostream &operator<<(std::ostream &out, const VectorConstView &V);
namespace vec
{
template<class Vec>
void
print(std::ostream &out, const Vec &V)
{
for (size_t i = 0; i < V.getSize(); i++)
out << V(i) << " ";
out << std::endl;
}
//! Computes v1 = v1 + v2
template<class Vec1, class Vec2>
void
add(Vec1 v1, const Vec2 v2)
{
assert(v1.getSize() == v2.getSize());
double *p1 = v1.getData();
const double *p2 = v2.getData();
while (p1 < v1.getData() + v1.getSize() * v1.getStride())
{
*p1 += *p2;
p1 += v1.getStride();
p2 += v2.getStride();
}
}
//! Computes v1 = v1 - v2
template<class Vec1, class Vec2>
void
sub(Vec1 v1, const Vec2 v2)
{
assert(v1.getSize() == v2.getSize());
double *p1 = v1.getData();
const double *p2 = v2.getData();
while (p1 < v1.getData() + v1.getSize() * v1.getStride())
{
*p1 -= *p2;
p1 += v1.getStride();
p2 += v2.getStride();
}
}
//! Does v = -v
template<class Vec>
void
negate(Vec v)
{
double *p = v.getData();
while (p < v.getData() + v.getSize() * v.getStride())
{
*p = -*p;
p += v.getStride();
}
}
} // End of namespace
#endif

9
mex/sources/estimation/libmat/tests/Makefile.am Normal file
View File

@ -0,0 +1,9 @@
check_PROGRAMS = test-qr test-gsd
test_qr_SOURCES = ../Matrix.cc ../Vector.cc ../QRDecomposition.cc test-qr.cc
test_qr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_qr_CPPFLAGS = -I.. -I../../../
test_gsd_SOURCES = ../Matrix.cc ../Vector.cc ../GeneralizedSchurDecomposition.cc test-gsd.cc
test_gsd_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_gsd_CPPFLAGS = -I.. -I../../../

61
mex/sources/estimation/libmat/tests/test-gsd.cc Normal file
View File

@ -0,0 +1,61 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include "GeneralizedSchurDecomposition.hh"
int
main(int argc, char **argv)
{
size_t n = 3;
double D_data[] = { 1, 2, 3,
4, 5, 6,
7, 8, 9 };
double E_data[] = { 1, -3, 4,
-7, 9, 1,
-3, 4, 0 };
MatrixView D(D_data, n, n, n), E(E_data, n, n, n);
// Need to transpose because internally matrices are in column-major order
mat::transpose(D);
mat::transpose(E);
std::cout << "Matrix D: " << std::endl << D << std::endl;
std::cout << "Matrix E: " << std::endl << E << std::endl;
GeneralizedSchurDecomposition GSD(n, 1.00001);
Matrix S(n), T(n), Z(n);
size_t sdim;
GSD.compute(D, E, S, T, Z, sdim);
std::cout << "Matrix S: " << std::endl << S << std::endl;
std::cout << "Matrix T: " << std::endl << T << std::endl;
std::cout << "Matrix Z: " << std::endl << Z << std::endl;
Vector eig_real(n), eig_cmplx(n);
GSD.getGeneralizedEigenvalues(eig_real, eig_cmplx);
std::cout << "Real part of generalized eigenvalues: " << std::endl << eig_real << std::endl;
std::cout << "Complex part of generalized eigenvalues: " << std::endl << eig_cmplx << std::endl;
return 0;
}

56
mex/sources/estimation/libmat/tests/test-qr.cc Normal file
View File

@ -0,0 +1,56 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include "BlasBindings.hh"
#include "QRDecomposition.hh"
int
main(int argc, char **argv)
{
size_t m = 4, n = 3;
Matrix S(m, n), Q(m), A(m, n), B(m);
QRDecomposition QRD(m, n);
for (size_t i = 0; i < m; i++)
for (size_t j = 0; j < n; j++)
S(i, j) = i*n + j + 1;
std::cout << "Matrix to be decomposed:" << std::endl << S << std::endl;
mat::set_identity(Q);
QRD.computeAndMultByQ(S, "L", "N", Q);
std::cout << "Matrix Q:" << std::endl << Q << std::endl;
blas::gemm("T", "N", 1.0, Q, Q, 0.0, B);
std::cout << "Matrix Q'*Q:" << std::endl << B << std::endl;
for (size_t j = 0; j < n; j++)
mat::col_set(S, j, j+1, m-j-1, 0);
std::cout << "Matrix R:" << std::endl << S << std::endl;
blas::gemm("N", "N", 1.0, Q, S, 0.0, A);
std::cout << "Product Q*R:" << std::endl << A << std::endl;
}

5
mex/sources/estimation/tests/Makefile.am Normal file
View File

@ -0,0 +1,5 @@
check_PROGRAMS = test-dr
test_dr_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../DecisionRules.cc test-dr.cc
test_dr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_dr_CPPFLAGS = -I.. -I../libmat -I../../

25
mex/sources/estimation/tests/test-dr.cc Normal file
View File

@ -0,0 +1,25 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "DecisionRules.hh"
int
main(int argc, char **argv)
{
}