Commit 0ed63ac0 by george

### 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`
parent a2448f63
 /* * 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); }
 /* * 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);
 ... ... @@ -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) ... ...
 /* * 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'
 #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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment