From 0ed63ac01afdf2e7196f56cc37f2ed309ddbaf4d Mon Sep 17 00:00:00 2001 From: george Date: Wed, 29 Jul 2009 15:01:56 +0000 Subject: [PATCH] New C++ disklyap_fast Lyapunov doublig solver based on previous work of J. Pearlman. git-svn-id: https://www.dynare.org/svn/dynare/trunk@2856 ac1d8469-bf42-47a9-8791-bf33cf982152 --- mex/sources/kalman/cc/disclyap_fast.cpp | 102 ++++++++++++++++++ mex/sources/kalman/cc/disclyap_fast.h | 34 ++++++ mex/sources/kalman/matlab/Makefile | 4 + .../kalman/matlab/disclyap_fast_dll.cpp | 89 +++++++++++++++ mex/sources/kalman/qt/cc/qt.h | 74 ++++++------- 5 files changed, 266 insertions(+), 37 deletions(-) create mode 100644 mex/sources/kalman/cc/disclyap_fast.cpp create mode 100644 mex/sources/kalman/cc/disclyap_fast.h create mode 100644 mex/sources/kalman/matlab/disclyap_fast_dll.cpp diff --git a/mex/sources/kalman/cc/disclyap_fast.cpp b/mex/sources/kalman/cc/disclyap_fast.cpp new file mode 100644 index 000000000..3b66c295c --- /dev/null +++ b/mex/sources/kalman/cc/disclyap_fast.cpp @@ -0,0 +1,102 @@ +/* +* Copyright (C) 2008-2009 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 . +*/ +/**************************************************************** +% function X=disclyap_fast(G,V,tol,ch) +% +% Solve the discrete Lyapunov Equation +% X=G*X*G'+V +% Using the Doubling Algorithm +% +% INPUT: +% G and V - square General matrices of same size +% tol - double tollerance level +% flag_ch - integer flag: if 1 check if the result is positive +% definite and generate an error message if it is not +% OUTPUT: +% on completion V - square General matrice contains solution +% +% based on work of Joe Pearlman and Alejandro Justiniano +% 3/5/2005 +% C++ version 28/07/09 by Dynare team +****************************************************************/ +#include "ts_exception.h" +#include "cppblas.h" +#include "GeneralMatrix.h" +//#include "Vector.h" +#include "SylvException.h" +#include "utils.h" +#include "mex.h" + +void disclyap_fast(const GeneralMatrix &G, const GeneralMatrix & V, GeneralMatrix &X, double tol = 1e-16, int flag_ch=0) + { + /** + if nargin == 2 | isempty( ch ) == 1 + flag_ch = 0; + else + flag_ch = 1; + end + **/ + //P0=V; + GeneralMatrix P0(V); + //A0=G; + GeneralMatrix A0(G); + + //n=size(G,1); + int n= A0.numCols(); + const double alpha=1.0; + const double half=0.5; + const double neg_alpha=-1.0; + const double omega=0.0; + + GeneralMatrix A1(n,n); + GeneralMatrix Ptmp(n,n); + GeneralMatrix P1(P0); + GeneralMatrix I(n,n); + I.unit(); + bool matd=true; + while (matd ) // matrix diff > tol + { + //P1=P0+A0*P0*A0'; + // first step Ptmp=P0*A0'; + // DGEMM: C := alpha*op( A )*op( B ) + beta*C, + BLAS_dgemm("N", "T", &n, &n, &n, &alpha, P0.base(), &n, + A0.base(), &n, &omega, Ptmp.base(), &n); + // P1=P0+A0*Ptmp; + BLAS_dgemm("N", "N", &n, &n, &n, &alpha, A0.base(), &n, + Ptmp.base(), &n, &alpha, P1.base(), &n); + // A1=A0*A0; + // A0=A1 (=A0*A0); + // A0.multRight(A0); + BLAS_dgemm("N", "N", &n, &n, &n, &alpha, A0.base(), &n, + A0.base(), &n, &omega, A1.base(), &n); + + // check if max( max( abs( P1 - P0 ) ) )>tol + matd=P0.isDiffSym(P1, tol); + P0=P1; + A0=A1; + }//end while + + // X=P0=(P0+P0')/2; + BLAS_dgemm("T", "N", &n, &n, &n, &half, P1.base(), &n, + I.base(), &n, &half, P0.base(), &n); + X=P0; + // Check that X is positive definite + if (flag_ch==1) + NormCholesky chol(P0); + } diff --git a/mex/sources/kalman/cc/disclyap_fast.h b/mex/sources/kalman/cc/disclyap_fast.h new file mode 100644 index 000000000..abf7bca34 --- /dev/null +++ b/mex/sources/kalman/cc/disclyap_fast.h @@ -0,0 +1,34 @@ +/* +* Copyright (C) 2008-2009 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 . +*/ +/**************************************************************** +% function X=disclyap_fast(G,V,ch) +% +% Solve the discrete Lyapunov Equation +% X=G*X*G'+V +% Using the Doubling Algorithm +% +% If ch is defined then the code will check if the resulting X +% is positive definite and generate an error message if it is not +% +% based on work of Joe Pearlman and Alejandro Justiniano +% 3/5/2005 +% C++ version 28/07/09 by Dynare team +****************************************************************/ +#include "GeneralMatrix.h" +void disclyap_fast(const GeneralMatrix &G, const GeneralMatrix & V, GeneralMatrix &X, double tol, int ch); diff --git a/mex/sources/kalman/matlab/Makefile b/mex/sources/kalman/matlab/Makefile index 665bd581a..bfd5947a1 100644 --- a/mex/sources/kalman/matlab/Makefile +++ b/mex/sources/kalman/matlab/Makefile @@ -101,6 +101,10 @@ qtmvm.dll: qtmvm.o $(KALMANLIB) # $(hsource) $(cppsource) gcc -shared $(CC_FLAGS) -o qtmvm.dll qtmvm.o \ kalmanlib.a $(LD_LIBS) +disclyap_fast_dll.dll: disclyap_fast_dll.o $(KALMANLIB) # $(hsource) $(cppsource) + gcc -shared $(CC_FLAGS) -o disclyap_fast_dll.dll disclyap_fast_dll.o \ + $(KALMANLIB) $(LD_LIBS) kalmanlib.def + kalman_filter_dll.dll: kalman_filters.o $(KALMANLIB) # $(hsource) $(cppsource) gcc -shared $(CC_FLAGS) -o kalman_filter_dll.dll kalman_filters.o \ $(KALMANLIB) $(LD_LIBS) diff --git a/mex/sources/kalman/matlab/disclyap_fast_dll.cpp b/mex/sources/kalman/matlab/disclyap_fast_dll.cpp new file mode 100644 index 000000000..7fe1ac832 --- /dev/null +++ b/mex/sources/kalman/matlab/disclyap_fast_dll.cpp @@ -0,0 +1,89 @@ +/* +* Copyright (C) 2008-2009 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 . +*/ +/**************************************************************** +% entry Matlab DLL for function X=disclyap_fast(G,V,ch) +% +% which solve the discrete Lyapunov Equation +% X=G*X*G'+V +% Using the Doubling Algorithm +% +% If ch is defined then the code will check if the resulting X +% is positive definite and generate an error message if it is not +% +****************************************************************/ +#include "ts_exception.h" + +#include "GeneralMatrix.h" +#include "SylvException.h" +#include "mex.h" +#include "disclyap_fast.h" + +//void disclyap_fast(GeneralMatrix &G,GeneralMatrix & V, double tol= 1e-16, int ch=0); + + +extern "C" { + void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) + { + + if (nrhs < 2 || nrhs > 4) + mexErrMsgTxt("Must have 2, 3 or 4 input parameters.\n"); + if (nlhs != 1 ) + mexErrMsgTxt("Must have 1 output parameters.\n"); + int cholCheck = 0; + double LyapTol=1e-06; + try + { + // make input matrices + int s = mxGetM(prhs[0]); + GeneralMatrix G(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0])); + GeneralMatrix V(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1])); + + // create output + plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]),mxGetN(prhs[0]), mxREAL); + GeneralMatrix X(mxGetPr(plhs[0]),mxGetM(plhs[0]),mxGetN(plhs[0])); + if (nrhs > 2) + LyapTol = (double)mxGetScalar(prhs[2]); + if (nrhs > 3) + cholCheck = (int)mxGetScalar(prhs[3]); + +#ifdef TIMING_LOOP + for (int tt=0;tt<1000;++tt) + { +#endif + disclyap_fast(G, V, X, LyapTol, cholCheck); +#ifdef TIMING_LOOP + } + mexPrintf("disclyap_fast: finished 1000 loops"); +#endif + + } + catch (const TSException& e) + { + mexErrMsgTxt(e.getMessage()); + } + catch (SylvException& e) + { + char mes[300]; + e.printMessage(mes, 299); + mexErrMsgTxt(mes); + } + + } // mexFunction + }; // extern 'C' diff --git a/mex/sources/kalman/qt/cc/qt.h b/mex/sources/kalman/qt/cc/qt.h index f0413922a..1e40ccfd6 100644 --- a/mex/sources/kalman/qt/cc/qt.h +++ b/mex/sources/kalman/qt/cc/qt.h @@ -1,37 +1,37 @@ - -#ifndef QT_H -#define QT_H - -#define C_CHAR const char* -#define BLINT int* -#define C_BLINT const int* -#define C_BLDOU const double* -#define BLDOU double* - - -extern "C"{ -void ldm_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT); -void ldsld_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT); -void ldv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void mtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void qt2ld_(BLDOU,C_BLDOU, C_BLINT); -void qt2t_(BLDOU, C_BLDOU, C_BLINT); -void s2d_(BLDOU,C_BLDOU, C_BLINT); -void s2u_(BLDOU,C_BLDOU, C_BLINT); -void td_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void tm_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void tstt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void tt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void tu_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void tut_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void tv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -void qtv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -#ifdef WINDOWS -void qtv_1__(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -#else -void qtv_1_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); -#endif -void qtsqtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); - -}; -#endif + +#ifndef QT_H +#define QT_H + +#define C_CHAR const char* +#define BLINT int* +#define C_BLINT const int* +#define C_BLDOU const double* +#define BLDOU double* + + +extern "C"{ +void ldm_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT); +void ldsld_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT); +void ldv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void mtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void qt2ld_(BLDOU,C_BLDOU, C_BLINT); +void qt2t_(BLDOU, C_BLDOU, C_BLINT); +void s2d_(BLDOU,C_BLDOU, C_BLINT); +void s2u_(BLDOU,C_BLDOU, C_BLINT); +void td_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void tm_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void tstt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void tt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void tu_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void tut_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void tv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +void qtv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +#ifdef WINDOWS +void qtv_1__(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +#else +void qtv_1_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); +#endif +void qtsqtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT); + +}; +#endif