Commit ca24c92e authored by Houtan Bastani's avatar Houtan Bastani
Browse files

A_times_B_kronecker_C: remove instances of mexErrMsgTxt

parent 7b0d6da9
......@@ -593,7 +593,9 @@ hu1 = [hu;zeros(np-npred,M_.exo_nbr)];
[nrhx,nchx] = size(hx);
[nrhu1,nchu1] = size(hu1);
B1 = B*A_times_B_kronecker_C(dr.ghxx,hx,hu1);
[err, abcOut] = A_times_B_kronecker_C(dr.ghxx,hx,hu1);
mexErrCheck('A_times_B_kronecker_C', err);
B1 = B*abcOut;
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
......@@ -608,7 +610,8 @@ kk = kk(1:npred,1:npred);
rhs = sparse_hessian_times_B_kronecker_C(hessian,zu);
B1 = A_times_B_kronecker_C(B*dr.ghxx,hu1);
[err, B1] = A_times_B_kronecker_C(B*dr.ghxx,hu1);
mexErrCheck('A_times_B_kronecker_C', err);
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
%lhs
......@@ -666,8 +669,10 @@ for i=1:M_.maximum_endo_lead
kk = find(kstate(:,2) == M_.maximum_endo_lag+i+1);
gu = dr.ghx*Gu;
[nrGu,ncGu] = size(Gu);
G1 = A_times_B_kronecker_C(dr.ghxx,Gu);
G2 = A_times_B_kronecker_C(hxx,Gu);
[err, G1] = A_times_B_kronecker_C(dr.ghxx,Gu);
mexErrCheck('A_times_B_kronecker_C', err);
[err, G2] = A_times_B_kronecker_C(hxx,Gu);
mexErrCheck('A_times_B_kronecker_C', err);
guu = dr.ghx*Guu+G1;
Gu = hx*Gu;
Guu = hx*Guu;
......
......@@ -347,7 +347,9 @@ hu1 = [hu;zeros(np-npred,M_.exo_nbr)];
[nrhx,nchx] = size(hx);
[nrhu1,nchu1] = size(hu1);
B1 = B*A_times_B_kronecker_C(dr.ghxx,hx,hu1);
[err, abcOut] = A_times_B_kronecker_C(dr.ghxx,hx,hu1);
mexErrCheck('A_times_B_kronecker_C', err);
B1 = B*abcOut;
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
......@@ -362,7 +364,8 @@ kk = kk(1:npred,1:npred);
rhs = sparse_hessian_times_B_kronecker_C(hessian,zu);
B1 = A_times_B_kronecker_C(B*dr.ghxx,hu1);
[err, B1] = A_times_B_kronecker_C(B*dr.ghxx,hu1);
mexErrCheck('A_times_B_kronecker_C', err);
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
%lhs
......@@ -420,8 +423,10 @@ for i=1:M_.maximum_endo_lead
kk = find(kstate(:,2) == M_.maximum_endo_lag+i+1);
gu = dr.ghx*Gu;
[nrGu,ncGu] = size(Gu);
G1 = A_times_B_kronecker_C(dr.ghxx,Gu);
G2 = A_times_B_kronecker_C(hxx,Gu);
[err, G1] = A_times_B_kronecker_C(dr.ghxx,Gu);
mexErrCheck('A_times_B_kronecker_C', err);
[err, G2] = A_times_B_kronecker_C(hxx,Gu);
mexErrCheck('A_times_B_kronecker_C', err);
guu = dr.ghx*Guu+G1;
Gu = hx*Gu;
Guu = hx*Guu;
......
function D = A_times_B_kronecker_C(A,B,C)
%function D = A_times_B_kronecker_C(A,B,C)
function [err, D] = A_times_B_kronecker_C(A,B,C)
%function [err, D] = A_times_B_kronecker_C(A,B,C)
% Computes A * kron(B,C).
%
% INPUTS
% A [double] mA*nA matrix.
% B [double] mB*nB matrix.
% C [double] mC*nC matrix.
%
% A [double] mA*nA matrix.
% B [double] mB*nB matrix.
% C [double] mC*nC matrix.
%
% OUTPUTS
% D [double] mA*(nC*nB) or mA*(nB*nB) matrix.
%
% err [double] scalar: 1 indicates failure, 0 indicates success
% D [double] mA*(nC*nB) or mA*(nB*nB) matrix.
%
% ALGORITHM
% none.
%
......@@ -34,11 +35,8 @@ function D = A_times_B_kronecker_C(A,B,C)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Chek number of inputs and outputs.
if nargin>3 || nargin<2
error('Two or Three input arguments required!')
end
if nargout>1
error('Too many output arguments!')
if nargin>3 || nargin<2 || nargout~=2
error('A_times_B_kronecker_C takes 2 or 3 input arguments and provides exactly 2 output arguments.')
end
% Get & check dimensions. Initialization of the output matrix.
[mA,nA] = size(A);
......@@ -82,4 +80,5 @@ else
else
D = A * kron(B,B);
end
end
\ No newline at end of file
end
err = 0;
\ No newline at end of file
......@@ -84,20 +84,32 @@ else
yhat1 = y__(dr.order_var(k2))-dr.ys(dr.order_var(k2));
yhat2 = y_(dr.order_var(k2),i-1)-dr.ys(dr.order_var(k2));
epsilon = ex_(i-1,:)';
[err, abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1);
mexErrCheck('A_times_B_kronecker_C', err);
[err, abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
[err, abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
y_(dr.order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
+ A_times_B_kronecker_C(.5*dr.ghxx,yhat1) ...
+ A_times_B_kronecker_C(.5*dr.ghuu,epsilon) ...
+ A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
+ abcOut1 + abcOut2 + abcOut3;
y__(dr.order_var) = dr.ys(dr.order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
end
else
for i = 2:iter+M_.maximum_lag
yhat = y_(dr.order_var(k2),i-1)-dr.ys(dr.order_var(k2));
epsilon = ex_(i-1,:)';
[err, abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat);
mexErrCheck('A_times_B_kronecker_C', err);
[err, abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
[err, abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
mexErrCheck('A_times_B_kronecker_C', err);
y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
+ A_times_B_kronecker_C(.5*dr.ghxx,yhat) ...
+ A_times_B_kronecker_C(.5*dr.ghuu,epsilon) ...
+ A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
+ abcOut1 + abcOut2 + abcOut3;
end
end
end
......
/*
* Copyright (C) 2007-2009 Dynare Team
* Copyright (C) 2007-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -124,14 +124,9 @@ void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
// Check input and output:
if ((nrhs > 3) || (nrhs < 2))
{
mexErrMsgTxt("Two or Three input arguments required.");
}
if (nlhs > 1)
{
mexErrMsgTxt("Too many output arguments.");
}
if (nrhs > 3 || nrhs < 2 || nlhs != 2)
DYN_MEX_FUNC_ERR_MSG_TXT("A_times_B_kronecker_C takes 2 or 3 input arguments and provides exactly 2 output arguments.");
// Get & Check dimensions (columns and rows):
mwSize mA, nA, mB, nB, mC, nC;
mA = mxGetM(prhs[0]);
......@@ -143,16 +138,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mC = mxGetM(prhs[2]);
nC = mxGetN(prhs[2]);
if (mB*mC != nA)
{
mexErrMsgTxt("Input dimension error!");
}
DYN_MEX_FUNC_ERR_MSG_TXT( "Input dimension error!");
}
else // A*kron(B,B) is to be computed.
{
if (mB*mB != nA)
{
mexErrMsgTxt("Input dimension error!");
}
DYN_MEX_FUNC_ERR_MSG_TXT( "Input dimension error!");
}
// Get input matrices:
double *B, *C, *A;
......@@ -166,13 +157,13 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
double *D;
if (nrhs == 3)
{
plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
plhs[1] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
}
else
{
plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
plhs[1] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
}
D = mxGetPr(plhs[0]);
D = mxGetPr(plhs[1]);
// Computational part:
if (nrhs == 2)
{
......@@ -182,4 +173,5 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
full_A_times_kronecker_B_C(A, B, C, &D[0], (int) mA, (int) nA, (int) mB, (int) nB, (int) mC, (int) nC);
}
plhs[0] = mxCreateDoubleScalar(0);
}
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