2010-11-19 14:39:54 +01:00
/* kalman_steady_state.cc
* *
2019-04-16 11:40:38 +02:00
* * Copyright © 2009 - 2013 Dynare Team .
2010-11-19 14:39:54 +01:00
* *
* * 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/>.
* *
* * This mex file calls fortran routines from the Slicot library .
*/
/*
2011-02-04 16:53:12 +01:00
+ + INPUTS
+ + = = = = = =
+ +
+ +
+ + [ 0 ] T ( double ) n - by - n transition matrix .
+ +
+ + [ 1 ] QQ ( double ) n - by - n matrix ( = R * Q * R ' , where Q is the covariance matrix of the structural innovations ) .
+ +
+ + [ 2 ] Z ( double ) n - by - p selection matrix .
+ +
+ + [ 3 ] H ( double ) p - by - p covariance matrix of the measurement errors .
+ +
+ +
+ +
+ +
+ + OUTPUTS
+ + = = = = = = =
+ +
+ +
+ + [ 0 ] P ( double ) n - by - n covariance matrix of the state vector .
+ +
+ +
+ + NOTES
+ + = = = = =
+ +
+ + [ 1 ] T = transpose ( dynare transition matrix ) and Z = transpose ( dynare selection matrix ) .
*/
2010-11-19 14:39:54 +01:00
# include <string.h>
# include <stdlib.h>
# include <dynmex.h>
2010-11-25 14:26:24 +01:00
# include <dynlapack.h>
2010-11-19 14:39:54 +01:00
2010-11-22 12:22:22 +01:00
# if !defined(MATLAB_MEX_FILE) || !defined(_WIN32)
2011-02-04 16:53:12 +01:00
# define sb02od sb02od_
2010-11-19 14:39:54 +01:00
# endif
extern " C "
{
2011-02-04 16:53:12 +01:00
int sb02od ( char * , char * , char * , char * , char * , char * , mwSize * , mwSize * , mwSize * , double * , mwSize * , double * , mwSize * , double * , mwSize * , double * , mwSize * , double * , mwSize * , double * , double * , mwSize * , double * , double * , double * , double * , mwSize * , double * , mwSize * , double * , mwSize * , double * , lapack_int * , double * , mwSize * , lapack_int * , lapack_int * ) ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
template < typename T >
T
max ( T x , T y )
2010-11-19 14:39:54 +01:00
{
2011-02-04 16:53:12 +01:00
return x < y ? y : x ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
template < typename T >
T
max ( T x , T y , T z )
2010-11-19 14:39:54 +01:00
{
return max ( x , max ( y , z ) ) ;
}
2011-02-04 16:53:12 +01:00
void
mexFunction ( int nlhs , mxArray * plhs [ ] , int nrhs , const mxArray * prhs [ ] )
2010-11-19 14:39:54 +01:00
{
// Check the number of arguments and set some flags.
int measurement_error_flag = 1 ;
2011-02-11 12:18:48 +01:00
if ( nrhs < 3 | | 4 < nrhs )
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state accepts either 3 or 4 input arguments! " ) ;
if ( nlhs < 1 | | 2 < nlhs )
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state requires at least 1, but no more than 2, output arguments! " ) ;
if ( nrhs = = 3 )
measurement_error_flag = 0 ;
2010-11-19 14:39:54 +01:00
// Check the type of the input arguments and get the size of the matrices.
mwSize n = mxGetM ( prhs [ 0 ] ) ;
2019-04-23 12:58:38 +02:00
if ( static_cast < size_t > ( n ) ! = mxGetN ( prhs [ 0 ] ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The first input argument (T) must be a square matrix! " ) ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
if ( ( mxIsNumeric ( prhs [ 0 ] ) = = 0 ) | | ( mxIsComplex ( prhs [ 0 ] ) = = 1 ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The first input argument (T) must be a real matrix! " ) ;
2010-11-19 14:39:54 +01:00
}
mwSize q = mxGetM ( prhs [ 1 ] ) ;
2019-04-23 12:58:38 +02:00
if ( static_cast < size_t > ( q ) ! = mxGetN ( prhs [ 1 ] ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The second input argument (QQ) must be a square matrix! " ) ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
if ( ( mxIsNumeric ( prhs [ 1 ] ) = = 0 ) | | ( mxIsComplex ( prhs [ 1 ] ) = = 1 ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The second input argument (QQ) must be a real matrix! " ) ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
if ( q ! = n )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The size of the second input argument (QQ) must match the size of the first argument (T)! " ) ;
2010-11-19 14:39:54 +01:00
}
mwSize p = mxGetN ( prhs [ 2 ] ) ;
2019-04-23 12:58:38 +02:00
if ( mxGetM ( prhs [ 2 ] ) ! = static_cast < size_t > ( n ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The number of rows of the third argument (Z) must match the number of rows of the first argument (T)! " ) ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
if ( ( mxIsNumeric ( prhs [ 2 ] ) = = 0 ) | | ( mxIsComplex ( prhs [ 2 ] ) = = 1 ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The third input argument (Z) must be a real matrix! " ) ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
if ( measurement_error_flag )
2010-11-19 14:39:54 +01:00
{
2011-02-04 16:53:12 +01:00
if ( mxGetM ( prhs [ 3 ] ) ! = mxGetN ( prhs [ 3 ] ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The fourth input argument (H) must be a square matrix! " ) ;
2010-11-19 14:39:54 +01:00
}
2019-04-23 12:58:38 +02:00
if ( mxGetM ( prhs [ 3 ] ) ! = static_cast < size_t > ( p ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The number of rows of the fourth input argument (H) must match the number of rows of the third input argument! " ) ;
2010-11-19 14:39:54 +01:00
}
2011-02-04 16:53:12 +01:00
if ( ( mxIsNumeric ( prhs [ 3 ] ) = = 0 ) | | ( mxIsComplex ( prhs [ 3 ] ) = = 1 ) )
2010-11-19 14:39:54 +01:00
{
2011-02-11 12:18:48 +01:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " kalman_steady_state: The fifth input argument (H) must be a real matrix! " ) ;
2010-11-19 14:39:54 +01:00
}
}
// Get input matrices.
2011-02-04 16:53:12 +01:00
double * T , * QQ , * Z , * H , * L ; // Remark. L will not be used.
2019-04-23 12:58:38 +02:00
T = static_cast < double * > ( mxCalloc ( n * n , sizeof ( double ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( T , mxGetPr ( prhs [ 0 ] ) , n * n * sizeof ( double ) ) ;
2019-04-23 12:58:38 +02:00
QQ = static_cast < double * > ( mxCalloc ( n * n , sizeof ( double ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( QQ , mxGetPr ( prhs [ 1 ] ) , n * n * sizeof ( double ) ) ;
2019-04-23 12:58:38 +02:00
Z = static_cast < double * > ( mxCalloc ( n * p , sizeof ( double ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( Z , mxGetPr ( prhs [ 2 ] ) , n * p * sizeof ( double ) ) ;
2019-04-23 12:58:38 +02:00
H = static_cast < double * > ( mxCalloc ( p * p , sizeof ( double ) ) ) ;
2010-11-19 14:39:54 +01:00
if ( measurement_error_flag )
{
2011-02-04 16:53:12 +01:00
memcpy ( H , mxGetPr ( prhs [ 3 ] ) , p * p * sizeof ( double ) ) ;
2010-11-19 14:39:54 +01:00
}
2019-04-23 12:58:38 +02:00
L = static_cast < double * > ( mxCalloc ( n * p , sizeof ( double ) ) ) ;
2011-02-04 16:53:12 +01:00
char * DICO , * JOBB , * FACT , * UPLO , * JOBL , * SORT ;
2019-04-23 12:58:38 +02:00
DICO = static_cast < char * > ( mxCalloc ( 2 , sizeof ( char ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( DICO , " D " , 2 * sizeof ( char ) ) ; // We want to solve a discrete Riccati equation.
2019-04-23 12:58:38 +02:00
JOBB = static_cast < char * > ( mxCalloc ( 2 , sizeof ( char ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( JOBB , " B " , 2 * sizeof ( char ) ) ; // Matrices Z and H are given.
2019-04-23 12:58:38 +02:00
FACT = static_cast < char * > ( mxCalloc ( 2 , sizeof ( char ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( FACT , " N " , 2 * sizeof ( char ) ) ; // Given matrices H and QQ are not factored.
2019-04-23 12:58:38 +02:00
UPLO = static_cast < char * > ( mxCalloc ( 2 , sizeof ( char ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( UPLO , " U " , 2 * sizeof ( char ) ) ; // Upper triangle of matrix H is stored.
2019-04-23 12:58:38 +02:00
JOBL = static_cast < char * > ( mxCalloc ( 2 , sizeof ( char ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( JOBL , " Z " , 2 * sizeof ( char ) ) ; // L matrix is zero.
2019-04-23 12:58:38 +02:00
SORT = static_cast < char * > ( mxCalloc ( 2 , sizeof ( char ) ) ) ;
2011-02-04 16:53:12 +01:00
memcpy ( SORT , " S " , 2 * sizeof ( char ) ) ; // Stable eigenvalues come first.
2010-11-19 14:39:54 +01:00
mwSize nn = 2 * n ;
2019-04-23 12:58:38 +02:00
mwSize LDA = max ( static_cast < mwSize > ( 1 ) , n ) ;
2010-11-19 14:39:54 +01:00
mwSize LDQ = LDA ;
2019-04-23 12:58:38 +02:00
mwSize LDU = max ( static_cast < mwSize > ( 1 ) , nn ) ;
mwSize LDS = max ( static_cast < mwSize > ( 1 ) , nn + p ) ;
mwSize LIWORK = max ( static_cast < mwSize > ( 1 ) , p , nn ) ;
mwSize LDR = max ( static_cast < mwSize > ( 1 ) , p ) ;
2010-11-19 14:39:54 +01:00
mwSize LDB = LDA ;
mwSize LDL = LDA ;
mwSize LDT = LDS ;
mwSize LDX = LDA ;
2019-04-23 12:58:38 +02:00
mwSize LDWORK = max ( static_cast < mwSize > ( 7 ) * ( static_cast < mwSize > ( 2 ) * n + static_cast < mwSize > ( 1 ) ) + static_cast < mwSize > ( 16 ) , static_cast < mwSize > ( 16 ) * n ) ;
LDWORK = max ( LDWORK , static_cast < mwSize > ( 2 ) * n + p , static_cast < mwSize > ( 3 ) * p ) ;
2011-06-28 14:48:36 +02:00
double tolerance = .0000000000000001 ;
2010-11-25 14:26:24 +01:00
lapack_int INFO ;
2010-11-19 14:39:54 +01:00
// Outputs of subroutine sb02OD
double rcond ;
double * WR , * WI , * BETA , * S , * TT , * UU ;
2019-04-23 12:58:38 +02:00
WR = static_cast < double * > ( mxCalloc ( nn , sizeof ( double ) ) ) ;
WI = static_cast < double * > ( mxCalloc ( nn , sizeof ( double ) ) ) ;
BETA = static_cast < double * > ( mxCalloc ( nn , sizeof ( double ) ) ) ;
S = static_cast < double * > ( mxCalloc ( LDS * ( nn + p ) , sizeof ( double ) ) ) ;
TT = static_cast < double * > ( mxCalloc ( LDT * nn , sizeof ( double ) ) ) ;
UU = static_cast < double * > ( mxCalloc ( LDU * nn , sizeof ( double ) ) ) ;
2010-11-19 14:39:54 +01:00
// Working arrays
2010-11-25 14:26:24 +01:00
lapack_int * IWORK ;
2019-04-23 12:58:38 +02:00
IWORK = static_cast < lapack_int * > ( mxCalloc ( LIWORK , sizeof ( lapack_int ) ) ) ;
2010-11-19 14:39:54 +01:00
double * DWORK ;
2019-04-23 12:58:38 +02:00
DWORK = static_cast < double * > ( mxCalloc ( LDWORK , sizeof ( double ) ) ) ;
2010-11-25 14:26:24 +01:00
lapack_int * BWORK ;
2019-04-23 12:58:38 +02:00
BWORK = static_cast < lapack_int * > ( mxCalloc ( nn , sizeof ( lapack_int ) ) ) ;
2010-11-19 14:39:54 +01:00
// Initialize the output of the mex file
double * P ;
2011-02-11 12:18:48 +01:00
plhs [ 1 ] = mxCreateDoubleMatrix ( n , n , mxREAL ) ;
P = mxGetPr ( plhs [ 1 ] ) ;
2010-11-19 14:39:54 +01:00
// Call the slicot routine
sb02od ( DICO , JOBB , FACT , UPLO , JOBL , SORT , & n , & p , & p , & T [ 0 ] , & LDA , & Z [ 0 ] , & LDB , & QQ [ 0 ] , & LDQ , & H [ 0 ] , & LDR , & L [ 0 ] , & LDL , & rcond , & P [ 0 ] , & LDX , & WR [ 0 ] , & WI [ 0 ] , & BETA [ 0 ] , & S [ 0 ] , & LDS , & TT [ 0 ] , & LDT , & UU [ 0 ] , & LDU , & tolerance , & IWORK [ 0 ] , & DWORK [ 0 ] , & LDWORK , & BWORK [ 0 ] , & INFO ) ;
mxFree ( T ) ;
mxFree ( QQ ) ;
mxFree ( Z ) ;
mxFree ( H ) ;
mxFree ( L ) ;
2011-02-04 16:53:12 +01:00
mxFree ( DICO ) ;
2010-11-19 14:39:54 +01:00
mxFree ( JOBB ) ;
mxFree ( FACT ) ;
mxFree ( UPLO ) ;
mxFree ( JOBL ) ;
mxFree ( SORT ) ;
mxFree ( WR ) ;
mxFree ( WI ) ;
mxFree ( BETA ) ;
mxFree ( S ) ;
mxFree ( TT ) ;
mxFree ( UU ) ;
mxFree ( IWORK ) ;
mxFree ( DWORK ) ;
mxFree ( BWORK ) ;
2011-02-04 16:53:12 +01:00
if ( INFO ! = 0 )
2010-11-19 14:39:54 +01:00
{
2011-02-04 16:53:12 +01:00
switch ( INFO )
2010-11-19 14:39:54 +01:00
{
case 1 :
{
2011-06-24 16:32:46 +02:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " The computed extended matrix pencil is singular, possibly due to rounding errors. \n " ) ;
2010-11-19 14:39:54 +01:00
break ;
}
case 2 :
{
2011-06-24 16:32:46 +02:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " The QZ (or QR) algorithm failed! \n " ) ;
2010-11-19 14:39:54 +01:00
break ;
}
2011-02-04 16:53:12 +01:00
case 3 :
2010-11-19 14:39:54 +01:00
{
2011-06-24 16:32:46 +02:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " The reordering of the (generalized) eigenvalues failed! \n " ) ;
2010-11-19 14:39:54 +01:00
break ;
}
case 4 :
{
2011-06-24 16:32:46 +02:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " After reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues \n in the (generalized) Schur form no longer satisfy the stability condition; this could also be caused due to scaling. " ) ;
2010-11-19 14:39:54 +01:00
break ;
}
case 5 :
{
2011-06-24 16:32:46 +02:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " The computed dimension of the solution does not equal n! \n " ) ;
2010-11-19 14:39:54 +01:00
break ;
}
case 6 :
{
2011-06-24 16:32:46 +02:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " A singular matrix was encountered during the computation of the solution matrix P! \n " ) ;
2010-11-19 14:39:54 +01:00
break ;
2011-02-04 16:53:12 +01:00
}
2010-11-22 12:22:45 +01:00
default :
2010-11-19 14:39:54 +01:00
{
2011-06-24 16:32:46 +02:00
DYN_MEX_FUNC_ERR_MSG_TXT ( " Unknown problem! \n " ) ;
2010-11-19 14:39:54 +01:00
}
}
}
2011-02-11 12:18:48 +01:00
plhs [ 0 ] = mxCreateDoubleScalar ( 0 ) ;
2010-11-19 14:39:54 +01:00
}