Commit 6479edf5 by ferhat

### - Correction of several bugs

```- normalize an equation linear in its endogenous variable
- Chained rule derivatives (necessary to reduce a block to the feedback equations and variables)

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2726 ac1d8469-bf42-47a9-8791-bf33cf982152```
parent d67f1ecd
 ... @@ -90,15 +90,16 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe ... @@ -90,15 +90,16 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe g1=spalloc( Blck_size, Blck_size, nze); g1=spalloc( Blck_size, Blck_size, nze); while ~(cvg==1 | iter>maxit_), while ~(cvg==1 | iter>maxit_), if(is_dynamic) if(is_dynamic) [r, g1, g2, g3] = feval(fname, y, x, params, it_, 0); [r, y, g1, g2, g3] = feval(fname, y, x, params, it_, 0); else else [r, g1, g2, g3] = feval(fname, y, x, params, 0); [r, y, g1, g2, g3] = feval(fname, y, x, params, 0); end; end; if(~isreal(r)) if(~isreal(r)) max_res=(-(max(max(abs(r))))^2)^0.5; max_res=(-(max(max(abs(r))))^2)^0.5; else else max_res=max(max(abs(r))); max_res=max(max(abs(r))); end; end; %['max_res=' num2str(max_res) ' Block_Num=' int2str(Block_Num) ' it_=' int2str(it_)] if(verbose==1) if(verbose==1) disp(['iteration : ' int2str(iter) ' => ' num2str(max_res) ' time = ' int2str(it_)]); disp(['iteration : ' int2str(iter) ' => ' num2str(max_res) ' time = ' int2str(it_)]); if(is_dynamic) if(is_dynamic) ... @@ -315,7 +316,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe ... @@ -315,7 +316,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, y_index_eq, nze, pe function [err, G]=local_fname(yl, x, params, y, y_index_eq, fname, is_csolve) function [err, G]=local_fname(yl, x, params, y, y_index_eq, fname, is_csolve) y(y_index_eq) = yl; y(y_index_eq) = yl; [err, G] = feval(fname, y, x, params, 0); [err, y, G] = feval(fname, y, x, params, 0); if(is_csolve) if(is_csolve) G = full(G); G = full(G); end; end;
 ... @@ -58,7 +58,7 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ ... @@ -58,7 +58,7 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ % You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . % along with Dynare. If not, see . global oo_; global oo_ M_; cvg=0; cvg=0; iter=0; iter=0; Per_u_=0; Per_u_=0; ... @@ -72,17 +72,18 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ ... @@ -72,17 +72,18 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods); g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods); reduced = 0; reduced = 0; while ~(cvg==1 | iter>maxit_), while ~(cvg==1 | iter>maxit_), [r, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, y_kmin, Blck_size); [r, y, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, y_kmin, Blck_size); g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size); g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size); %disp(['size(g1)=' int2str(size(g1))]); b = b -g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)'; %disp(['g1(:,' int2str(1:y_kmin_l*Blck_size) ')']); [max_res, max_indx]=max(max(abs(r'))); %disp(['g1(:,' int2str((periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size) ')']); b = b' -g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)'; if(~isreal(r)) if(~isreal(r)) max_res=(-(max(max(abs(r))))^2)^0.5; max_res = (-max_res^2)^0.5; else max_res=max(max(abs(r))); end; end; % if(~isreal(r)) % max_res=(-(max(max(abs(r))))^2)^0.5; % else % max_res=max(max(abs(r))); % end; if(~isreal(max_res) | isnan(max_res)) if(~isreal(max_res) | isnan(max_res)) cvg = 0; cvg = 0; elseif(is_linear & iter>0) elseif(is_linear & iter>0) ... @@ -93,6 +94,9 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ ... @@ -93,6 +94,9 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ if(~cvg) if(~cvg) if(iter>0) if(iter>0) if(~isreal(max_res) | isnan(max_res) | (max_resa1)) if(~isreal(max_res) | isnan(max_res) | (max_resa1)) if(~isreal(max_res)) disp(['Variable ' M_.endo_names(max_indx,:) ' (' int2str(max_indx) ') returns an undefined value']); end; if(isnan(max_res)) if(isnan(max_res)) detJ=det(g1aa); detJ=det(g1aa); if(abs(detJ)<1e-7) if(abs(detJ)<1e-7) ... ...
 ... @@ -16,9 +16,8 @@ ... @@ -16,9 +16,8 @@ * You should have received a copy of the GNU General Public License * You should have received a copy of the GNU General Public License * along with Dynare. If not, see . * along with Dynare. If not, see . */ */ //#define DEBUGC #include #include #include "Interpreter.hh" #include "Interpreter.hh" Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *direction_arg, int y_size_arg, Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *direction_arg, int y_size_arg, ... @@ -77,16 +76,22 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ ... @@ -77,16 +76,22 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ case FLDV : case FLDV : //load a variable in the processor //load a variable in the processor #ifdef DEBUGC #ifdef DEBUGC mexPrintf("FLDV"); if(Block_Count==2) mexEvalString("drawnow;"); { mexPrintf("FLDV\n"); mexEvalString("drawnow;"); } #endif #endif switch (get_code_char) switch (get_code_char) { { case eParameter : case eParameter : var=get_code_int var=get_code_int #ifdef DEBUGC #ifdef DEBUGC mexPrintf(" params[%d]=%f\n",var,params[var]); if(Block_Count==2) mexEvalString("drawnow;"); { mexPrintf(" params[%d]=%f\n",var,params[var]); mexEvalString("drawnow;"); } #endif #endif Stack.push(params[var]); Stack.push(params[var]); break; break; ... @@ -94,27 +99,32 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ ... @@ -94,27 +99,32 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ var=get_code_int var=get_code_int lag=get_code_int lag=get_code_int #ifdef DEBUGC #ifdef DEBUGC if(var==153) if(Block_Count==2) { { mexPrintf(" FLD y[var=%d,time=%d,lag=%d,%d]=%f\n",var,it_,lag,(it_+lag)*y_size+var,y[(it_+lag)*y_size+var]); mexPrintf("y[%d, %d]=%f\n",it_+lag, var+1, y[(it_+lag)*y_size+var]); mexEvalString("drawnow;"); } } #endif #endif Stack.push(y[(it_+lag)*y_size+var]); Stack.push(y[(it_+lag)*y_size+var]); break; break; case eExogenous : case eExogenous : #ifdef DEBUGC mexPrintf("Exogenous\n"); #endif var=get_code_int var=get_code_int lag=get_code_int lag=get_code_int #ifdef DEBUGC #ifdef DEBUGC if(var==650 or var==643 or var==628) if(Block_Count==2) { { mexPrintf(" FLD x[%d, time=%d, var=%d, lag=%d]=%f\n",it_+lag+var*nb_row_x,it_,var,lag,x[it_+lag+var*nb_row_x]); mexPrintf("x[%d, %d]\n",it_+lag, var+1); mexEvalString("drawnow;"); } } #endif #endif Stack.push(x[it_+lag+var*nb_row_x]); Stack.push(x[it_+lag+var*nb_row_x]); break; break; case eExogenousDet : case eExogenousDet : #ifdef DEBUGC mexPrintf("ExogenousDet\n"); #endif var=get_code_int var=get_code_int lag=get_code_int lag=get_code_int #ifdef DEBUGC #ifdef DEBUGC ... @@ -123,6 +133,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ ... @@ -123,6 +133,8 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ #endif #endif Stack.push(x[it_+lag+var*nb_row_xd]); Stack.push(x[it_+lag+var*nb_row_xd]); break; break; default: mexPrintf("Unknown vraibale type\n"); } } break; break; case FLDT : case FLDT : ... @@ -183,39 +195,53 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ ... @@ -183,39 +195,53 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/ case eParameter : case eParameter : var=get_code_int var=get_code_int params[var] = Stack.top(); params[var] = Stack.top(); #ifdef DEBUGC mexPrintf("FSTP params[%d]=%f\n", var, params[var]); mexEvalString("drawnow;"); #endif Stack.pop(); Stack.pop(); break; break; case eEndogenous : case eEndogenous : #ifdef DEBUGC mexPrintf("FSTP Endogenous\n"); #endif var=get_code_int var=get_code_int lag=get_code_int lag=get_code_int #ifdef DEBUGC #ifdef DEBUGC mexPrintf("y[%d(it_=%d, lag=%d, y_size=%d, var=%d)](%d)=",(it_+lag)*y_size+var,it_, lag, y_size, var, Stack.size()); //mexPrintf("y[%d(it_=%d, lag=%d, y_size=%d, var=%d)](%d)=",(it_+lag)*y_size+var,it_, lag, y_size, var, Stack.size()); mexPrintf("y[it_=%d, lag=%d, y_size=%d, var=%d, block=%d)]=",it_, lag, y_size, var+1, Block_Count+1); mexEvalString("drawnow;"); mexEvalString("drawnow;"); #endif #endif y[(it_+lag)*y_size+var] = Stack.top(); y[(it_+lag)*y_size+var] = Stack.top(); #ifdef DEBUGC #ifdef DEBUGC if(var==557 || var==558) mexPrintf("%f\n",y[(it_+lag)*y_size+var]); { mexEvalString("drawnow;"); mexPrintf(" FSTP y[var=%d,time=%d,lag=%d,%d]=%f\n",var,it_,lag,(it_+lag)*y_size+var,y[(it_+lag)*y_size+var]); mexEvalString("drawnow;"); } /*mexPrintf("%f\n",y[(it_+lag)*y_size+var]); mexEvalString("drawnow;");*/ #endif #endif Stack.pop(); Stack.pop(); break; break; case eExogenous : case eExogenous : #ifdef DEBUGC mexPrintf("Exogenous\n"); #endif var=get_code_int var=get_code_int var=get_code_int lag=get_code_int lag=get_code_int x[it_+lag+var*nb_row_x] = Stack.top(); x[it_+lag+var*nb_row_x] = Stack.top(); Stack.pop(); Stack.pop(); break; break; case eExogenousDet : case eExogenousDet : #ifdef DEBUGC mexPrintf("ExogenousDet\n"); #endif var=get_code_int var=get_code_int var=get_code_int lag=get_code_int lag=get_code_int x[it_+lag+var*nb_row_xd] = Stack.top(); x[it_+lag+var*nb_row_xd] = Stack.top(); Stack.pop(); Stack.pop(); break; break; default: mexPrintf("Unknown vraibale type\n"); } } break; break; case FSTPT : case FSTPT : ... @@ -584,7 +610,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba ... @@ -584,7 +610,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba Mat_DP a; Mat_DP a; Vec_INT indx; Vec_INT indx; #endif #endif //SparseMatrix sparse_matrix; //SparseMatrix sparse_matrix; //int nb_endo, u_count_init; //int nb_endo, u_count_init; ... @@ -599,19 +625,22 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba ... @@ -599,19 +625,22 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba switch (type) switch (type) { { case EVALUATE_FORWARD : case EVALUATE_FORWARD : case EVALUATE_FORWARD_R : //case EVALUATE_FORWARD_R : #ifdef DEBUGC #ifdef DEBUGC mexPrintf("EVALUATE_FORWARD\n"); mexPrintf("EVALUATE_FORWARD\n"); #endif #endif begining=get_code_pointer; begining=get_code_pointer; for (it_=y_kmin;it_y_kmin;it_--) for (it_=periods+y_kmin-1;it_>=y_kmin;it_--) { { set_code_pointer(begining); set_code_pointer(begining); Per_y_=it_*y_size; Per_y_=it_*y_size; ... @@ -824,6 +853,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba ... @@ -824,6 +853,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba if (!is_linear) if (!is_linear) { { max_res_idx=0; for (it_=y_kmin;it_ in SOLVE_FORWARD_COMPLETE before compute_block_time OK\n"); /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE before compute_block_time OK\n"); mexEvalString("drawnow;");*/ mexEvalString("drawnow;");*/ compute_block_time(0); compute_block_time(0); ... @@ -916,6 +949,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba ... @@ -916,6 +949,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba begining=get_code_pointer; begining=get_code_pointer; if (!is_linear) if (!is_linear) { { max_res_idx=0; for (it_=periods+y_kmin;it_>y_kmin;it_--) for (it_=periods+y_kmin;it_>y_kmin;it_--) { { cvg=false; cvg=false; ... @@ -937,7 +971,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba ... @@ -937,7 +971,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba double rr; double rr; rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]); if (max_res
 ... @@ -453,7 +453,7 @@ void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::m ... @@ -453,7 +453,7 @@ void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::m it4=IM.begin(); it4=IM.begin(); eq=-1; eq=-1; double tmp_b[Size]; double tmp_b[Size]; #pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ///#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) for(i=0; i< Size;i++) for(i=0; i< Size;i++) { { tmp_b[i]=0;//u[i]; tmp_b[i]=0;//u[i]; ... @@ -509,7 +509,7 @@ void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::m ... @@ -509,7 +509,7 @@ void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::m } } it4++; it4++; } } #pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ///#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) for(i=0;i