Changed the order of the output arguments in kronecker product related mex routines.
Stéphane Adjemian (Charybdis)
Stéphane Adjemian (Scylla)
parent
9b9fb7df3a
commit
560800bcc5
16
matlab/dr1.m
16
matlab/dr1.m
|
|
@ -555,7 +555,7 @@ zx=[zx; zeros(M_.exo_nbr,np);zeros(M_.exo_det_nbr,np)];
|
|||
zu=[zu; eye(M_.exo_nbr);zeros(M_.exo_det_nbr,M_.exo_nbr)];
|
||||
[nrzx,nczx] = size(zx);
|
||||
|
||||
[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zx,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zx,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
|
||||
rhs = -rhs;
|
||||
|
||||
|
|
@ -615,7 +615,7 @@ hu = dr.ghu(nstatic+1:nstatic+npred,:);
|
|||
%kk = kk(1:npred,1:npred);
|
||||
%rhs = -hessian*kron(zx,zu)-f1*dr.ghxx(end-nyf+1:end,kk(:))*kron(hx(1:npred,:),hu(1:npred,:));
|
||||
|
||||
[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zx,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zx,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
|
||||
|
||||
nyf1 = sum(kstate(:,2) == M_.maximum_endo_lag+2);
|
||||
|
|
@ -624,7 +624,7 @@ hu1 = [hu;zeros(np-npred,M_.exo_nbr)];
|
|||
[nrhx,nchx] = size(hx);
|
||||
[nrhu1,nchu1] = size(hu1);
|
||||
|
||||
[err, abcOut] = A_times_B_kronecker_C(dr.ghxx,hx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[abcOut,err] = A_times_B_kronecker_C(dr.ghxx,hx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
B1 = B*abcOut;
|
||||
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
|
||||
|
|
@ -638,10 +638,10 @@ dr.ghxu = A\rhs;
|
|||
kk = reshape([1:np*np],np,np);
|
||||
kk = kk(1:npred,1:npred);
|
||||
|
||||
[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
|
||||
|
||||
[err, B1] = A_times_B_kronecker_C(B*dr.ghxx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[B1, err] = A_times_B_kronecker_C(B*dr.ghxx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
|
||||
|
||||
|
|
@ -686,7 +686,7 @@ for i=1:M_.maximum_endo_lead
|
|||
[junk,k3a,k3] = ...
|
||||
find(M_.lead_lag_incidenceordered(M_.maximum_endo_lag+j+1,:));
|
||||
nk3a = length(k3a);
|
||||
[err, B1] = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k3,k3)),gu(k3a,:),options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
[B1, err] = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k3,k3)),gu(k3a,:),options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
|
||||
RHS = RHS + jacobia_(:,k2)*guu(k2a,:)+B1;
|
||||
end
|
||||
|
|
@ -701,9 +701,9 @@ 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);
|
||||
[err, G1] = A_times_B_kronecker_C(dr.ghxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[G1, err] = A_times_B_kronecker_C(dr.ghxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err, G2] = A_times_B_kronecker_C(hxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[G2, err] = A_times_B_kronecker_C(hxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
guu = dr.ghx*Guu+G1;
|
||||
Gu = hx*Gu;
|
||||
|
|
|
|||
|
|
@ -77,11 +77,11 @@ u = oo.exo_simul(1,:)';
|
|||
|
||||
Uyy = full(Uyy);
|
||||
|
||||
[err,Uyygygy] = A_times_B_kronecker_C(Uyy,gy,gy,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Uyygygy, err] = A_times_B_kronecker_C(Uyy,gy,gy,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Uyygugu] = A_times_B_kronecker_C(Uyy,gu,gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Uyygugu, err] = A_times_B_kronecker_C(Uyy,gu,gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Uyygygu] = A_times_B_kronecker_C(Uyy,gy,gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Uyygygu, err] = A_times_B_kronecker_C(Uyy,gy,gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
|
||||
Wbar =U/(1-beta);
|
||||
|
|
@ -92,9 +92,9 @@ if LQ
|
|||
else
|
||||
Wyy = (Uy*gyy+Uyygygy+beta*Wy*Gyy)/(eye(npred*npred)-beta*kron(Gy,Gy));
|
||||
end
|
||||
[err,Wyygugu] = A_times_B_kronecker_C(Wyy,Gu,Gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Wyygugu, err] = A_times_B_kronecker_C(Wyy,Gu,Gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Wyygygu] = A_times_B_kronecker_C(Wyy,Gy,Gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Wyygygu,err] = A_times_B_kronecker_C(Wyy,Gy,Gu,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
if LQ
|
||||
Wuu = Uyygugu+beta*Wyygugu;
|
||||
|
|
@ -114,11 +114,11 @@ end
|
|||
yhat = yhat(dr.order_var(nstatic+(1:npred)),1)-dr.ys(dr.order_var(nstatic+(1:npred)));
|
||||
u = oo.exo_simul(1,:)';
|
||||
|
||||
[err,Wyyyhatyhat] = A_times_B_kronecker_C(Wyy,yhat,yhat,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Wyyyhatyhat, err] = A_times_B_kronecker_C(Wyy,yhat,yhat,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Wuuuu] = A_times_B_kronecker_C(Wuu,u,u,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Wuuuu, err] = A_times_B_kronecker_C(Wuu,u,u,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err,Wyuyhatu] = A_times_B_kronecker_C(Wyu,yhat,u,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
[Wyuyhatu, err] = A_times_B_kronecker_C(Wyu,yhat,u,options.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
planner_objective_value(1) = Wbar+Wy*yhat+Wu*u+Wyuyhatu ...
|
||||
+ 0.5*(Wyyyhatyhat + Wuuuu+Wss);
|
||||
|
|
|
|||
|
|
@ -1,4 +1,39 @@
|
|||
function [err, D] = A_times_B_kronecker_C(A,B,C,fake)
|
||||
function [D, err] = A_times_B_kronecker_C(A,B,C,fake)
|
||||
|
||||
%@info:
|
||||
%! @deftypefn {Function File} {@var{D}, @var{err} =} A_timesB_kronecker_C (@var{A},@var{B},@var{C},@var{fake})
|
||||
%! @anchor{kronecker/A_times_B_kronecker_C}
|
||||
%! @sp 1
|
||||
%! Computes A*kron(B,C).
|
||||
%! @sp 2
|
||||
%! @strong{Inputs}
|
||||
%! @sp 1
|
||||
%! @table @ @var
|
||||
%! @item A
|
||||
%! mA*nA matrix of doubles.
|
||||
%! @item B
|
||||
%! mB*nB matrix of doubles.
|
||||
%! @item C
|
||||
%! mC*nC matrix of doubles.
|
||||
%! @item fake
|
||||
%! Anything you want, just a fake parameter (because the mex version admits a last argument specifying the number of threads to be used in parallel mode).
|
||||
%! @end table
|
||||
%! @sp 2
|
||||
%! @strong{Outputs}
|
||||
%! @sp 1
|
||||
%! @table @ @var
|
||||
%! @item c
|
||||
%! Integer scalar, the number of years, quarters, months or weeks between @var{a} and @var{B}.
|
||||
%! @end table
|
||||
%! @sp 2
|
||||
%! @strong{This function is called by:}
|
||||
%! @sp 2
|
||||
%! @strong{This function calls:}
|
||||
%! @ref{@@dynDates/eq},@ref{@@dynDates/lt}
|
||||
%!
|
||||
%! @end deftypefn
|
||||
%@eod:
|
||||
|
||||
%function [err, D] = A_times_B_kronecker_C(A,B,C)
|
||||
% Computes A * kron(B,C).
|
||||
%
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
function [err, D] = sparse_hessian_times_B_kronecker_C(varargin)
|
||||
function [D, err] = sparse_hessian_times_B_kronecker_C(varargin)
|
||||
%function [err, D] = sparse_hessian_times_B_kronecker_C(A,B,C, fake)
|
||||
% Computes A * kron(B,C) where A is a sparse matrix.
|
||||
%
|
||||
|
|
@ -43,9 +43,9 @@ end
|
|||
|
||||
switch nargin
|
||||
case 4
|
||||
[fake,D] = A_times_B_kronecker_C(A,B,C,fake);
|
||||
[D, fake] = A_times_B_kronecker_C(A,B,C,fake);
|
||||
case 3
|
||||
[fake,D] = A_times_B_kronecker_C(A,B,C);
|
||||
[D, fake] = A_times_B_kronecker_C(A,B,C);
|
||||
otherwise
|
||||
error('Two or Three input arguments required!')
|
||||
end
|
||||
|
|
|
|||
|
|
@ -105,11 +105,11 @@ else
|
|||
yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
|
||||
yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
|
||||
epsilon = ex_(i-1,:)';
|
||||
[err, abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err, abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err, abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
|
||||
+ abcOut1 + abcOut2 + abcOut3;
|
||||
|
|
@ -119,11 +119,11 @@ else
|
|||
for i = 2:iter+M_.maximum_lag
|
||||
yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
|
||||
epsilon = ex_(i-1,:)';
|
||||
[err, abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err, abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
[err, abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
[abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
|
||||
mexErrCheck('A_times_B_kronecker_C', err);
|
||||
y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
|
||||
+ abcOut1 + abcOut2 + abcOut3;
|
||||
|
|
|
|||
|
|
@ -126,8 +126,8 @@ void
|
|||
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
||||
{
|
||||
// Check input and output:
|
||||
if (nrhs > 4 || nrhs < 3 || nlhs != 2)
|
||||
DYN_MEX_FUNC_ERR_MSG_TXT("A_times_B_kronecker_C takes 3 or 4 input arguments and provides exactly 2 output arguments.");
|
||||
if (nrhs > 4 || nrhs < 3)
|
||||
DYN_MEX_FUNC_ERR_MSG_TXT("A_times_B_kronecker_C takes 3 or 4 input arguments and provides 2 output arguments.");
|
||||
|
||||
// Get & Check dimensions (columns and rows):
|
||||
mwSize mA, nA, mB, nB, mC, nC;
|
||||
|
|
@ -164,13 +164,13 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
double *D;
|
||||
if (nrhs == 4)
|
||||
{
|
||||
plhs[1] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
}
|
||||
else
|
||||
{
|
||||
plhs[1] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
}
|
||||
D = mxGetPr(plhs[1]);
|
||||
D = mxGetPr(plhs[0]);
|
||||
// Computational part:
|
||||
if (nrhs == 3)
|
||||
{
|
||||
|
|
@ -180,5 +180,5 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
{
|
||||
full_A_times_kronecker_B_C(A, B, C, &D[0], mA, nA, mB, nB, mC, nC, numthreads);
|
||||
}
|
||||
plhs[0] = mxCreateDoubleScalar(0);
|
||||
plhs[1] = mxCreateDoubleScalar(0);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -143,8 +143,8 @@ void
|
|||
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
||||
{
|
||||
// Check input and output:
|
||||
if ((nrhs > 4) || (nrhs < 3) || nlhs != 2)
|
||||
DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C takes 3 or 4 input arguments and provides exactly 2 output argument.");
|
||||
if ((nrhs > 4) || (nrhs < 3) )
|
||||
DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C takes 3 or 4 input arguments and provides 2 output arguments.");
|
||||
|
||||
if (!mxIsSparse(prhs[0]))
|
||||
DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C: First input must be a sparse (dynare) hessian matrix.");
|
||||
|
|
@ -185,13 +185,13 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
double *D;
|
||||
if (nrhs == 4)
|
||||
{
|
||||
plhs[1] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
}
|
||||
else
|
||||
{
|
||||
plhs[1] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
}
|
||||
D = mxGetPr(plhs[1]);
|
||||
D = mxGetPr(plhs[0]);
|
||||
// Computational part:
|
||||
if (nrhs == 3)
|
||||
{
|
||||
|
|
@ -201,5 +201,5 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
{
|
||||
sparse_hessian_times_B_kronecker_C(isparseA, jsparseA, vsparseA, B, C, D, mA, nA, mB, nB, mC, nC, numthreads);
|
||||
}
|
||||
plhs[0] = mxCreateDoubleScalar(0);
|
||||
plhs[1] = mxCreateDoubleScalar(0);
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue