### - Correction of several bugs in sparse_dll

```- Pound expressions accepted with sparse option

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2329 ac1d8469-bf42-47a9-8791-bf33cf982152```
parent 7e23ff1e
 ... ... @@ -66,7 +66,7 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i % along with Dynare. If not, see . global oo_ M_; global oo_ M_ options_; Blck_size=size(y_index_eq,2); g2 = []; g3 = []; ... ... @@ -90,9 +90,9 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i g1=spalloc( Blck_size, Blck_size, nze); while ~(cvg==1 | iter>maxit_), if(is_dynamic) [r, g1] = feval(fname, y, x, params, it_, 0, g1, g2, g3); [r, g1, g2, g3, b] = feval(fname, y, x, params, it_, 0, g1, g2, g3); else [r, g1] = feval(fname, y, x, params, 0); [r, g1, g2, g3, b] = feval(fname, y, x, params, 0); end; if(~isreal(r)) max_res=(-(max(max(abs(r))))^2)^0.5; ... ... @@ -100,10 +100,16 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i max_res=max(max(abs(r))); end; if(verbose==1) disp(['iteration : ' int2str(iter) ' => ' num2str(max_res)]); disp([M_.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]); disp(['iteration : ' int2str(iter) ' => ' num2str(max_res) ' time = ' int2str(it_)]); if(is_dynamic) disp([M_.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])]); else disp([M_.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]); end; end; if(is_linear) if(~isreal(max_res) | isnan(max_res)) cvg = 0; elseif(is_linear & iter>0) cvg = 1; else cvg=(max_res 0 info = 0; else info = 1; end elseif(~is_dynamic & options_.solve_algo==2) lambda=1; stpmx = 100 ; stpmax = stpmx*max([sqrt(y'*y);size(y_index_eq,2)]); nn=1:size(y_index_eq,2); g = (r'*g1)'; f = 0.5*r'*r; p = -g1\r ; [y,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, 0); elseif(~is_dynamic & options_.solve_algo==3) [yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1); y(y_index_eq) = yn; elseif((simulation_method==0 & is_dynamic) | (~is_dynamic & options_.solve_algo==1)), dx = g1\r; ya = ya - lambda*dx; if(is_dynamic) y(it_,y_index_eq) = ya'; else y(y_index_eq) = ya; end; elseif(simulation_method==2), elseif(simulation_method==2 & is_dynamic), flag1=1; while(flag1>0) [L1, U1]=luinc(g1,luinc_tol); [za,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1); %[za,flag1] = gmres(-g1,b',Blck_size,1e-6,Blck_size,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]); ... ... @@ -204,11 +241,12 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i end; end; end; elseif(simulation_method==3), elseif(simulation_method==3 & is_dynamic), flag1=1; while(flag1>0) [L1, U1]=luinc(g1,luinc_tol); [za,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1); %[za,flag1] = bicgstab(-g1,b',1e-7,Blck_size,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]); ... ... @@ -229,6 +267,14 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i end; end; end; else disp('unknown option : '); if(is_dynamic) disp(['options_.simulation_method = ' num2str(simulation_method) ' not implemented']); else disp(['options_.solve_algo = ' num2str(options_.solve_algo) ' not implemented']); end; return; end; iter=iter+1; end ... ... @@ -254,8 +300,15 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i oo_.deterministic_simulation.status = 1; oo_.deterministic_simulation.error = max_res; oo_.deterministic_simulation.iterations = iter; oo_.deterministic_simulation.block(Block_Num).status = 1;% Convergency failed. oo_.deterministic_simulation.block(Block_Num).status = 1; oo_.deterministic_simulation.block(Block_Num).error = max_res; oo_.deterministic_simulation.block(Block_Num).iterations = iter; end; return; \ No newline at end of file return; function [err, G]=local_fname(yl, x, params, y, y_index_eq, fname, is_csolve) y(y_index_eq) = yl; [err, G] = feval(fname, y, x, params, 0); if(is_csolve) G = full(G); end;
 ... ... @@ -32,7 +32,7 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ % - 3 BicGStab % % OUTPUTS % y [matrix] All endogenous variables of the model % y [matrix] All endogenous variables of the model % % ALGORITHM % Newton with LU or GMRES or BicGstab ... ... @@ -72,7 +72,6 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods); reduced = 0; while ~(cvg==1 | iter>maxit_), %fname [r, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, g1, g2, g3, y_kmin, Blck_size); g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size); b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+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)'; ... ... @@ -81,113 +80,117 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ else max_res=max(max(abs(r))); end; if(iter>0) if(~isreal(max_res) | isnan(max_res) | (max_resa1)) if(isnan(max_res)) detJ=det(g1aa); if(abs(detJ)<1e-7) max_factor=max(max(abs(g1aa))); ze_elem=sum(diag(g1aa)0) cvg = 1; else cvg=(max_res0) if(~isreal(max_res) | isnan(max_res) | (max_resa1)) if(isnan(max_res)) detJ=det(g1aa); if(abs(detJ)<1e-7) max_factor=max(max(abs(g1aa))); ze_elem=sum(diag(g1aa)1e-8) lambda=lambda/2; reduced = 1; disp(['reducing the path length: lambda=' num2str(lambda,'%f')]); y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)'; continue; else if(cutoff == 0) fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_".\n',Block_Num, iter); else disp('The singularity of the jacobian matrix could not be corrected'); return; fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_" or set "cutoff=0" in model options.\n',Block_Num, iter); end; oo_.deterministic_simulation.status = 0; oo_.deterministic_simulation.error = max_res; oo_.deterministic_simulation.iterations = iter; oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed. oo_.deterministic_simulation.block(Block_Num).error = max_res; oo_.deterministic_simulation.block(Block_Num).iterations = iter; return; end; elseif(lambda>1e-8) lambda=lambda/2; reduced = 1; disp(['reducing the path length: lambda=' num2str(lambda,'%f')]); y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)'; continue; else if(cutoff == 0) fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_".\n',Block_Num, iter); else fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_" or set "cutoff=0" in model options.\n',Block_Num, iter); if(lambda<1) lambda=max(lambda*2, 1); end; oo_.deterministic_simulation.status = 0; oo_.deterministic_simulation.error = max_res; oo_.deterministic_simulation.iterations = iter; oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed. oo_.deterministic_simulation.block(Block_Num).error = max_res; oo_.deterministic_simulation.block(Block_Num).iterations = iter; return; end; else if(lambda<1) lambda=max(lambda*2, 1); end; end; end; ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)'; ya_save=ya; g1aa=g1a; ba=b; max_resa=max_res; if(simulation_method==0), dx = g1a\b- ya; ya = ya + lambda*dx; y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; elseif(simulation_method==2), flag1=1; while(flag1>0) [L1, U1]=luinc(g1a,luinc_tol); [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); elseif(flag1==2) disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)'; ya_save=ya; g1aa=g1a; ba=b; max_resa=max_res; if(simulation_method==0), dx = g1a\b- ya; ya = ya + lambda*dx; y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; elseif(simulation_method==2), flag1=1; while(flag1>0) [L1, U1]=luinc(g1a,luinc_tol); [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); elseif(flag1==2) disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); end; luinc_tol = luinc_tol/10; reduced = 0; else dx = za - ya; ya = ya + lambda*dx; y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; end; luinc_tol = luinc_tol/10; reduced = 0; else dx = za - ya; ya = ya + lambda*dx; y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; end; end; elseif(simulation_method==3), flag1=1; while(flag1>0) [L1, U1]=luinc(g1a,luinc_tol); [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); elseif(flag1==2) disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); elseif(simulation_method==3), flag1=1; while(flag1>0) [L1, U1]=luinc(g1a,luinc_tol); [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); elseif(flag1==2) disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); end; luinc_tol = luinc_tol/10; reduced = 0; else dx = za - ya; ya = ya + lambda*dx; y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; end; luinc_tol = luinc_tol/10; reduced = 0; else dx = za - ya; ya = ya + lambda*dx; y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; end; end; end; if(is_linear) cvg = 1; else cvg=(max_resmaxit_) disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')]); oo_.deterministic_simulation.status = 0; ... ...
No preview for this file type
No preview for this file type
 ... ... @@ -106,7 +106,7 @@ Interpreter::compute_block_time() /*throw(EvalException)*/ var=get_code_int lag=get_code_int #ifdef DEBUGC if(var==6) if(var==650 or var==643 or var==628) { 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]); mexEvalString("drawnow;"); ... ... @@ -194,13 +194,13 @@ Interpreter::compute_block_time() /*throw(EvalException)*/ #endif y[(it_+lag)*y_size+var] = Stack.top(); #ifdef DEBUGC if(var==153) if(var==557 || var==558) { 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;"); /*mexPrintf("%f\n",y[(it_+lag)*y_size+var]); mexEvalString("drawnow;");*/ #endif Stack.pop(); break; ... ... @@ -686,6 +686,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba set_code_pointer(begining); Per_y_=it_*y_size; compute_block_time(); mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_SIMPLE : OK\n"); mexEvalString("drawnow;"); y[Per_y_+Block_Contain.Variable] += -r/g1; double rr; rr=r/(1+y[Per_y_+Block_Contain.Variable]); ... ... @@ -705,7 +707,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mxFree(g1); mxFree(r); break; case SOLVE_TWO_BOUNDARIES_SIMPLE : /*case SOLVE_TWO_BOUNDARIES_SIMPLE : #ifdef DEBUGC mexPrintf("SOLVE_TWO_BOUNDARIES_SIMPLE\n"); #endif ... ... @@ -777,17 +779,17 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba cvg=(max_res in SOLVE_FORWARD_COMPLETE : OK\n"); mexEvalString("drawnow;"); Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); res2=0; res1=0; ... ... @@ -833,8 +837,17 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba { set_code_pointer(begining); Per_y_=it_*y_size; iter = 0; res1=res2=max_res=0; /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE before compute_block_time OK\n"); mexEvalString("drawnow;");*/ compute_block_time(); Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : %d OK\n",it_); mexEvalString("drawnow;");*/ //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, /*true*/false, cvg, iter); //simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter); } } mxFree(g1); ... ... @@ -866,6 +879,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba { set_code_pointer(begining); compute_block_time(); mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_COMPLETE : OK\n"); mexEvalString("drawnow;"); Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); res2=0; res1=0; ... ... @@ -897,13 +912,14 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba set_code_pointer(begining); Per_y_=it_*y_size; compute_block_time(); Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 1, false, iter); } } mxFree(g1); mxFree(r); mxFree(u); break; case SOLVE_TWO_BOUNDARIES_SIMPLE : case SOLVE_TWO_BOUNDARIES_COMPLETE: #ifdef DEBUGC mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n"); ... ... @@ -948,11 +964,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mexPrintf("u=%x\n",u); #endif Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax); //mexPrintf("aft reading_sparse_matrix\n"); //mexEvalString("drawnow;"); u_count=u_count_int*(periods+y_kmax+y_kmin); g1=(double*)mxMalloc(size*size*sizeof(double)); //g1=(double*)mxMalloc(size*size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double)); y_save=(double*)mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); #ifdef DEBUGC ... ... @@ -1093,8 +1108,9 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba //mexErrMsgTxt("End of simulate"); #endif mxFree(g1); //mxFree(g1); mxFree(r); mxFree(y_save); mxFree(u); mxFree(index_vara); memset(direction,0,size_of_direction); ... ... @@ -1113,8 +1129,6 @@ Interpreter::compute_blocks(string file_name, string bin_basename) { ifstream CompiledCode; int Code_Size, var; /*mexPrintf("compute_blocks %s\n",file_name.c_str()); mexEvalString("drawnow;");*/ //First read and store inn memory the code CompiledCode.open((file_name + ".cod").c_str(),std::ios::in | std::ios::binary| std::ios::ate); if (!CompiledCode.is_open()) ... ... @@ -1132,8 +1146,6 @@ Interpreter::compute_blocks(string file_name, string bin_basename) mexEvalString("drawnow;"); #endif CompiledCode.seekg(std::ios::beg); /*Code_Size=CompiledCode.tellg(); mexPrintf("Code_Size=%d\n",Code_Size);*/ Code=(char*)mxMalloc(Code_Size); CompiledCode.seekg(0); CompiledCode.read(reinterpret_cast(Code), Code_Size); ... ...
 ... ... @@ -61,25 +61,32 @@ NonZeroElem* Mem_Mngr::mxMalloc_NZE() { int i; if (!Chunk_Stack.empty()) if (!Chunk_Stack.empty()) /*An unused block of memory available inside the heap*/ { NonZeroElem* p1 = Chunk_Stack.back(); Chunk_Stack.pop_back(); return(p1); } else if (CHUNK_heap_pos