diff --git a/dynare++/sylv/matlab/gensylv.cpp b/dynare++/sylv/matlab/gensylv.cpp index 9cfad4300..0a0d7e391 100644 --- a/dynare++/sylv/matlab/gensylv.cpp +++ b/dynare++/sylv/matlab/gensylv.cpp @@ -2,7 +2,7 @@ /* Tag $Name: $ */ - +#include "dynmex.h" #include "mex.h" #include "GeneralSylvester.h" @@ -13,14 +13,8 @@ void gen_sylv_solve(int order, int n, int m, int zero_cols, const double* A, const double* B, const double* C, double* X) { - try { GeneralSylvester s(order, n, m, zero_cols, A, B, C, X, false); s.solve(); - } catch (const SylvException& e) { - char mes[1000]; - e.printMessage(mes, 999); - mexErrMsgTxt(mes); - } } void gen_sylv_solve_and_check(int order, int n, int m, int zero_cols, @@ -28,44 +22,18 @@ void gen_sylv_solve_and_check(int order, int n, int m, int zero_cols, const double* C, const double* D, double* X, mxArray*& p) { - try { GeneralSylvester s(order, n, m, zero_cols, A, B, C, X, true); s.solve(); s.check(D); p = s.getParams().createStructArray(); - } catch (const SylvException& e) { - char mes[1000]; - e.printMessage(mes, 999); - mexErrMsgTxt(mes); - } -} - -void checkDimensions(const mwSize* const Adims, const mwSize* const Bdims, - const mwSize* const Cdims, const mwSize* const Ddims, - int order) -{ - if (Adims[0] != Adims[1]) - mexErrMsgTxt("Matrix A must be a square matrix."); - if (Adims[0] != Bdims[0]) - mexErrMsgTxt("Matrix A and matrix B must have the same number of rows."); - if (Adims[0] != Ddims[0]) - mexErrMsgTxt("Matrix A and matrix B must have the same number of rows."); - if (Cdims[0] != Cdims[1]) - mexErrMsgTxt("Matrix C must be square."); - if (Bdims[0] < Bdims[1]) - mexErrMsgTxt("Matrix B must not have more columns than rows."); - if (Ddims[1] != power(Cdims[0], order)) - mexErrMsgTxt("Matrix D has wrong number of columns."); } extern "C" { - void mexFunction(int nhls, mxArray* plhs[], + void mexFunction(int nlhs, mxArray* plhs[], int nhrs, const mxArray* prhs[]) { - if (nhrs != 5) - mexErrMsgTxt("Must have exactly 5 input args."); - if (nhls !=1 && nhls != 2) - mexErrMsgTxt("Must have 1 or 2 output args."); + if (nhrs != 5 || nlhs > 3 || nlhs < 2 ) + DYN_MEX_FUNC_ERR_MSG_TXT("Gensylv: Must have exactly 5 input args and either 2 or 3 output args."); int order = (int)mxGetScalar(prhs[0]); const mxArray* const A = prhs[1]; @@ -76,7 +44,20 @@ extern "C" { const mwSize* const Bdims = mxGetDimensions(B); const mwSize* const Cdims = mxGetDimensions(C); const mwSize* const Ddims = mxGetDimensions(D); - checkDimensions(Adims, Bdims, Cdims, Ddims, order); + + if (Adims[0] != Adims[1]) + DYN_MEX_FUNC_ERR_MSG_TXT("Matrix A must be a square matrix."); + if (Adims[0] != Bdims[0]) + DYN_MEX_FUNC_ERR_MSG_TXT("Matrix A and matrix B must have the same number of rows."); + if (Adims[0] != Ddims[0]) + DYN_MEX_FUNC_ERR_MSG_TXT("Matrix A and matrix B must have the same number of rows."); + if (Cdims[0] != Cdims[1]) + DYN_MEX_FUNC_ERR_MSG_TXT("Matrix C must be square."); + if (Bdims[0] < Bdims[1]) + DYN_MEX_FUNC_ERR_MSG_TXT("Matrix B must not have more columns than rows."); + if (Ddims[1] != power(Cdims[0], order)) + DYN_MEX_FUNC_ERR_MSG_TXT("Matrix D has wrong number of columns."); + int n = Adims[0]; int m = Cdims[0]; int zero_cols = Bdims[0] - Bdims[1]; @@ -86,15 +67,25 @@ extern "C" { ConstVector Dvec((double*)mxGetPr(D), power(m, order)*n); Xvec = Dvec; // solve (or solve and check) - if (nhls == 1) { - gen_sylv_solve(order, n, m, zero_cols, + try + { + if (nlhs == 2) { + gen_sylv_solve(order, n, m, zero_cols, mxGetPr(A), mxGetPr(B), mxGetPr(C), mxGetPr(X)); - } else if (nhls == 2) { - gen_sylv_solve_and_check(order, n, m, zero_cols, + } else if (nlhs == 3) { + gen_sylv_solve_and_check(order, n, m, zero_cols, mxGetPr(A), mxGetPr(B), mxGetPr(C), - mxGetPr(D), mxGetPr(X), plhs[1]); - } - plhs[0] = X; + mxGetPr(D), mxGetPr(X), plhs[2]); + } + } + catch (const SylvException& e) + { + char mes[1000]; + e.printMessage(mes, 999); + DYN_MEX_FUNC_ERR_MSG_TXT(mes); + } + plhs[1] = X; + plhs[0] = mxCreateDoubleScalar(0); } }; diff --git a/dynare++/sylv/matlab/gensylv.m b/dynare++/sylv/matlab/gensylv.m index b6cb0dccb..56f5f0e73 100644 --- a/dynare++/sylv/matlab/gensylv.m +++ b/dynare++/sylv/matlab/gensylv.m @@ -12,10 +12,11 @@ % ...... Kronecker power of C of order 'order' % D ..... rectangular (n, m^order) matrix. % -% X = gensylv(order, A, B, C, D) -% returns X as the solution, doesn't perform any checks +% [err, X] = gensylv(order, A, B, C, D) +% returns err a scalar where 1 indicates failure, 0 indicates success +% and X as the solution, doesn't perform any checks % -% [X, par] = gensylv(order, A, B, C, D) +% [err, X, par] = gensylv(order, A, B, C, D) % solves the system, and performs checks. 'par' is struct % containing information about solution and error norms % returned by the check. This is a list of the struct @@ -52,7 +53,7 @@ % $Header: /var/lib/cvs/dynare_cpp/sylv/matlab/gensylv.m,v 1.1.1.1 2004/06/04 13:01:13 kamenik Exp $ % Tag $Name: $ -function [X, varargout] = gensylv(order, A, B, C, D) +function [err, X, varargout] = gensylv(order, A, B, C, D) % in Windows, ensure that aa_gensylv.dll is loaded, this prevents % clearing the function and a successive Matlab crash @@ -63,12 +64,13 @@ if strcmp('PCWIN', computer) end % launch aa_gensylv -if nargout == 1 +if nargout == 2 X = aa_gensylv(order, A, B, C, D); -elseif nargout == 2 +elseif nargout == 3 [X, par] = aa_gensylv(order, A, B, C, D); varargout(1) = {par}; else - error('Must have 1 or 2 output arguments.'); + error('Must have 2 or 3 output arguments.'); end +err = 0; \ No newline at end of file diff --git a/matlab/dr1.m b/matlab/dr1.m index 1323e6559..7996169d9 100644 --- a/matlab/dr1.m +++ b/matlab/dr1.m @@ -578,7 +578,8 @@ C = hx; D = [rhs; zeros(n-M_.endo_nbr,size(rhs,2))]; -dr.ghxx = gensylv(2,A,B,C,D); +[err, dr.ghxx] = gensylv(2,A,B,C,D); +mexErrCheck('gensylv', err); %ghxu %rhs diff --git a/matlab/dr11_sparse.m b/matlab/dr11_sparse.m index ad48a362e..41a5c2fda 100644 --- a/matlab/dr11_sparse.m +++ b/matlab/dr11_sparse.m @@ -332,7 +332,8 @@ C = hx; D = [rhs; zeros(n-M_.endo_nbr,size(rhs,2))]; -dr.ghxx = gensylv(2,A,B,C,D); +[err, dr.ghxx] = gensylv(2,A,B,C,D); +mexErrCheck('gensylv', err); %ghxu %rhs diff --git a/matlab/gensylv/gensylv.m b/matlab/gensylv/gensylv.m index 5ccacb271..16e0d559b 100644 --- a/matlab/gensylv/gensylv.m +++ b/matlab/gensylv/gensylv.m @@ -1,5 +1,5 @@ -function E = gensylv(fake,A,B,C,D) -%function E = gensylv(fake,A,B,C,D) +function [err, E] = gensylv(fake,A,B,C,D) +%function [err, E] = gensylv(fake,A,B,C,D) % Solves a Sylvester equation. % % INPUTS @@ -10,6 +10,7 @@ function E = gensylv(fake,A,B,C,D) % D % % OUTPUTS +% err [double] scalar: 1 indicates failure, 0 indicates success % E % % ALGORITHM @@ -18,7 +19,7 @@ function E = gensylv(fake,A,B,C,D) % SPECIAL REQUIREMENTS % none. -% Copyright (C) 1996-2008 Dynare Team +% Copyright (C) 1996-2010 Dynare Team % % This file is part of Dynare. % @@ -38,3 +39,4 @@ function E = gensylv(fake,A,B,C,D) C = kron(C,C); x0 = sylvester3(A,B,C,D); E = sylvester3a(x0,A,B,C,D); +err = 0; \ No newline at end of file