diff --git a/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc
index 0b9e5d01b1cbd83f59cf08f453139aacddef6d6e..fbacf73bb4b1a57cba07f80144c4f64c10a6a4ab 100644
--- a/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc
+++ b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
*
* This file is part of Dynare.
*
@@ -24,7 +24,7 @@
double GeneralizedSchurDecomposition::criterium_static;
-GeneralizedSchurDecomposition::GeneralizedSchurDecomposition(size_t n_arg, double criterium_arg) :
+GeneralizedSchurDecomposition::GeneralizedSchurDecomposition(ptrdiff_t n_arg, double criterium_arg) :
n(n_arg), criterium(criterium_arg)
{
alphar = new double[n];
diff --git a/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh
index f776ba3330a7f86d4eda22326c23cde261412515..cd4dc3d5782542c0fe6d72b7c05b333f597142bd 100644
--- a/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh
+++ b/mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
*
* This file is part of Dynare.
*
@@ -19,13 +19,14 @@
#include <dynlapack.h>
-#include "Vector.hh"
-#include "Matrix.hh"
+#include <Eigen/Core>
+
+using namespace Eigen;
class GeneralizedSchurDecomposition
{
private:
- const size_t n;
+ const ptrdiff_t n;
const double criterium;
lapack_int lwork;
double *alphar, *alphai, *beta, *vsl, *work;
@@ -40,38 +41,44 @@ public:
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(ptrdiff_t n_arg, double criterium_arg);
virtual ~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>
+ template<typename Mat1, typename Mat2, typename Mat3>
+ void compute(PlainObjectBase<Mat1> &S, PlainObjectBase<Mat2> &T,
+ PlainObjectBase<Mat3> &Z, ptrdiff_t &sdim) throw (GSDException);
+
/*!
\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);
+ template<typename Mat1, typename Mat2, typename Mat3, typename Mat4, typename Mat5>
+ void compute(const MatrixBase<Mat1> &D, const MatrixBase<Mat2> &E,
+ PlainObjectBase<Mat3> &S, PlainObjectBase<Mat4> &T,
+ PlainObjectBase<Mat5> &Z, ptrdiff_t &sdim) throw (GSDException);
+
+ template<typename Vec1, typename Vec2>
+ void getGeneralizedEigenvalues(PlainObjectBase<Vec1> &eig_real,
+ PlainObjectBase<Vec2> &eig_cmplx);
};
std::ostream &operator<<(std::ostream &out, const GeneralizedSchurDecomposition::GSDException &e);
-template<class Mat1, class Mat2, class Mat3>
+template<typename Mat1, typename Mat2, typename Mat3>
void
-GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException)
+GeneralizedSchurDecomposition::compute(PlainObjectBase<Mat1> &S, PlainObjectBase<Mat2> &T, PlainObjectBase<Mat3> &Z, ptrdiff_t &sdim) throw (GSDException)
{
- assert(S.getRows() == n && S.getCols() == n
- && T.getRows() == n && T.getCols() == n
- && Z.getRows() == n && Z.getCols() == n);
+ assert(S.rows() == n && S.cols() == n
+ && T.rows() == n && T.cols() == n
+ && Z.rows() == n && Z.cols() == n);
lapack_int n2 = n;
lapack_int info, sdim2;
- lapack_int lds = S.getLd(), ldt = T.getLd(), ldz = Z.getLd();
+ lapack_int lds = S.outerStride(), ldt = T.outerStride(), ldz = Z.outerStride();
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,
+ dgges("N", "V", "S", &selctg, &n2, S.data(), &lds, T.data(), &ldt,
+ &sdim2, alphar, alphai, beta, vsl, &n2, Z.data(), &ldz,
work, &lwork, bwork, &info);
if (info != 0)
@@ -80,26 +87,27 @@ GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim)
sdim = sdim2;
}
-template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5>
+template<typename Mat1, typename Mat2, typename Mat3, typename Mat4, typename Mat5>
void
-GeneralizedSchurDecomposition::compute(const Mat1 &D, const Mat2 &E,
- Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException)
+GeneralizedSchurDecomposition::compute(const MatrixBase<Mat1> &D, const MatrixBase<Mat2> &E,
+ PlainObjectBase<Mat3> &S, PlainObjectBase<Mat4> &T, PlainObjectBase<Mat5> &Z, ptrdiff_t &sdim) throw (GSDException)
{
- assert(D.getRows() == n && D.getCols() == n
- && E.getRows() == n && E.getCols() == n);
+ assert(D.rows() == n && D.cols() == n
+ && E.rows() == n && E.cols() == n);
S = D;
T = E;
compute(S, T, Z, sdim);
}
-template<class Vec1, class Vec2>
+template<typename Vec1, typename Vec2>
void
-GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx)
+GeneralizedSchurDecomposition::getGeneralizedEigenvalues(PlainObjectBase<Vec1> &eig_real,
+ PlainObjectBase<Vec2> &eig_cmplx)
{
- assert(eig_real.getSize() == n && eig_cmplx.getSize() == n);
+ assert(eig_real.size() == n && eig_cmplx.size() == n);
double *par = alphar, *pai = alphai, *pb = beta,
- *per = eig_real.getData(), *pei = eig_cmplx.getData();
+ *per = eig_real.data(), *pei = eig_cmplx.data();
while (par < alphar + n)
{
*per = *par / *pb;
@@ -110,7 +118,7 @@ GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &e
par++;
pai++;
pb++;
- per += eig_real.getStride();
- pei += eig_cmplx.getStride();
+ per += eig_real.innerStride();
+ pei += eig_cmplx.innerStride();
}
}
diff --git a/mex/sources/estimation/libmat/tests/Makefile.am b/mex/sources/estimation/libmat/tests/Makefile.am
index 33d4a988389e990495db478fc84a649002d8494a..d2524b684293688ab30297153ad83dd3e7b48d26 100644
--- a/mex/sources/estimation/libmat/tests/Makefile.am
+++ b/mex/sources/estimation/libmat/tests/Makefile.am
@@ -4,7 +4,7 @@ test_qr_SOURCES = ../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_SOURCES = ../GeneralizedSchurDecomposition.cc test-gsd.cc
test_gsd_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_gsd_CPPFLAGS = -I.. -I../../../
diff --git a/mex/sources/estimation/libmat/tests/test-gsd.cc b/mex/sources/estimation/libmat/tests/test-gsd.cc
index fd8dea7f651e6d6e6d96344d9054337ead5a7055..ac0c9758ba5dc9c11a8cd7ab57d8502bad668c7c 100644
--- a/mex/sources/estimation/libmat/tests/test-gsd.cc
+++ b/mex/sources/estimation/libmat/tests/test-gsd.cc
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
*
* This file is part of Dynare.
*
@@ -24,26 +24,22 @@
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);
+ Matrix3d D, E;
+ D << 1, 2, 3,
+ 4, 5, 6,
+ 7, 8, 9;
+
+ E << 1, -3, 4,
+ -7, 9, 1,
+ -3, 4, 0;
std::cout << "D =" << std::endl << D << std::endl;
std::cout << "E =" << std::endl << E << std::endl;
- GeneralizedSchurDecomposition GSD(n, 1.00001);
+ GeneralizedSchurDecomposition GSD(3, 1.00001);
- Matrix S(n), T(n), Z(n);
- size_t sdim;
+ Matrix3d S, T, Z;
+ ptrdiff_t sdim;
GSD.compute(D, E, S, T, Z, sdim);
@@ -51,7 +47,7 @@ main(int argc, char **argv)
std::cout << "T =" << std::endl << T << std::endl;
std::cout << "Z =" << std::endl << Z << std::endl;
- Vector eig_real(n), eig_cmplx(n);
+ Vector3d eig_real, eig_cmplx;
GSD.getGeneralizedEigenvalues(eig_real, eig_cmplx);
std::cout << "Real part of generalized eigenvalues: " << std::endl << eig_real << std::endl;