- extends steady algorithms to solve_algo=5, 6, 7 and 8 for a block decomposed model without bytcode
- solve_algo is reordered = * 0: fsolve * 1: solve1 * 2, 4: solve1 + block decomposition * 3: csolve * 5: bytecode own solver (use Gaussian elimination + sparse matrix) * 6: LU decomposition with UMFPack (method handling sparse matrix in Matlab) * 7: GMRES * 8: BiCGStabtime-shift
parent
fc31afb356
commit
39718c8645
|
@ -97,7 +97,7 @@ for it_=start:incr:finish
|
||||||
if(is_dynamic)
|
if(is_dynamic)
|
||||||
[r, y, 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, y, g1, g2, g3] = feval(fname, y, x, params, 0);
|
[r, y, g1] = feval(fname, y, x, params);
|
||||||
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;
|
||||||
|
@ -219,6 +219,9 @@ for it_=start:incr:finish
|
||||||
ya_save=ya;
|
ya_save=ya;
|
||||||
g1a=g1;
|
g1a=g1;
|
||||||
if(~is_dynamic & options_.solve_algo == 0)
|
if(~is_dynamic & options_.solve_algo == 0)
|
||||||
|
if (verbose == 1)
|
||||||
|
disp("steady: fsolve");
|
||||||
|
end
|
||||||
if exist('OCTAVE_VERSION') || isempty(ver('optim'))
|
if exist('OCTAVE_VERSION') || isempty(ver('optim'))
|
||||||
% Note that fsolve() exists under Octave, but has a different syntax
|
% Note that fsolve() exists under Octave, but has a different syntax
|
||||||
% So we fail for the moment under Octave, until we add the corresponding code
|
% So we fail for the moment under Octave, until we add the corresponding code
|
||||||
|
@ -238,6 +241,9 @@ for it_=start:incr:finish
|
||||||
info = -Block_Num*10;
|
info = -Block_Num*10;
|
||||||
end
|
end
|
||||||
elseif((~is_dynamic & options_.solve_algo==2) || (is_dynamic & stack_solve_algo==4))
|
elseif((~is_dynamic & options_.solve_algo==2) || (is_dynamic & stack_solve_algo==4))
|
||||||
|
if (verbose == 1 & ~is_dynamic)
|
||||||
|
disp("steady: LU + lnsrch1");
|
||||||
|
end
|
||||||
lambda=1;
|
lambda=1;
|
||||||
stpmx = 100 ;
|
stpmx = 100 ;
|
||||||
if (is_dynamic)
|
if (is_dynamic)
|
||||||
|
@ -261,10 +267,16 @@ for it_=start:incr:finish
|
||||||
y = ya';
|
y = ya';
|
||||||
end;
|
end;
|
||||||
elseif(~is_dynamic & options_.solve_algo==3)
|
elseif(~is_dynamic & options_.solve_algo==3)
|
||||||
|
if (verbose == 1)
|
||||||
|
disp("steady: csolve");
|
||||||
|
end
|
||||||
[yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1);
|
[yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1);
|
||||||
dx = ya - yn;
|
dx = ya - yn;
|
||||||
y(y_index_eq) = yn;
|
y(y_index_eq) = yn;
|
||||||
elseif((stack_solve_algo==1 & is_dynamic) | (stack_solve_algo==0 & is_dynamic) | (~is_dynamic & options_.solve_algo==1)),
|
elseif((stack_solve_algo==1 & is_dynamic) | (stack_solve_algo==0 & is_dynamic) | (~is_dynamic & (options_.solve_algo==1 | options_.solve_algo==6))),
|
||||||
|
if (verbose == 1 & ~is_dynamic)
|
||||||
|
disp("steady: Sparse LU ");
|
||||||
|
end
|
||||||
dx = g1\r;
|
dx = g1\r;
|
||||||
ya = ya - lambda*dx;
|
ya = ya - lambda*dx;
|
||||||
if(is_dynamic)
|
if(is_dynamic)
|
||||||
|
@ -272,8 +284,14 @@ for it_=start:incr:finish
|
||||||
else
|
else
|
||||||
y(y_index_eq) = ya;
|
y(y_index_eq) = ya;
|
||||||
end;
|
end;
|
||||||
elseif(stack_solve_algo==2 & is_dynamic),
|
elseif((stack_solve_algo==2 & is_dynamic) | (options_.solve_algo==7 & ~is_dynamic)),
|
||||||
flag1=1;
|
flag1=1;
|
||||||
|
if exist('OCTAVE_VERSION')
|
||||||
|
error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=7 since GMRES is not implemented in Octave')
|
||||||
|
end
|
||||||
|
if (verbose == 1 & ~is_dynamic)
|
||||||
|
disp("steady: GMRES ");
|
||||||
|
end
|
||||||
while(flag1>0)
|
while(flag1>0)
|
||||||
[L1, U1]=luinc(g1,luinc_tol);
|
[L1, U1]=luinc(g1,luinc_tol);
|
||||||
[dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
|
[dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
|
||||||
|
@ -296,8 +314,11 @@ for it_=start:incr:finish
|
||||||
end;
|
end;
|
||||||
end;
|
end;
|
||||||
end;
|
end;
|
||||||
elseif(stack_solve_algo==3 & is_dynamic),
|
elseif((stack_solve_algo==3 & is_dynamic) | (options_.solve_algo==8 & ~is_dynamic)),
|
||||||
flag1=1;
|
flag1=1;
|
||||||
|
if (verbose == 1 & ~is_dynamic)
|
||||||
|
disp("steady: BiCGStab");
|
||||||
|
end
|
||||||
while(flag1>0)
|
while(flag1>0)
|
||||||
[L1, U1]=luinc(g1,luinc_tol);
|
[L1, U1]=luinc(g1,luinc_tol);
|
||||||
[dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
|
[dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
|
||||||
|
@ -352,14 +373,16 @@ for it_=start:incr:finish
|
||||||
return;
|
return;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
info = 1;
|
|
||||||
if(is_dynamic)
|
if(is_dynamic)
|
||||||
|
info = 1;
|
||||||
oo_.deterministic_simulation.status = 1;
|
oo_.deterministic_simulation.status = 1;
|
||||||
oo_.deterministic_simulation.error = max_res;
|
oo_.deterministic_simulation.error = max_res;
|
||||||
oo_.deterministic_simulation.iterations = iter;
|
oo_.deterministic_simulation.iterations = iter;
|
||||||
oo_.deterministic_simulation.block(Block_Num).status = 1;
|
oo_.deterministic_simulation.block(Block_Num).status = 1;
|
||||||
oo_.deterministic_simulation.block(Block_Num).error = max_res;
|
oo_.deterministic_simulation.block(Block_Num).error = max_res;
|
||||||
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
|
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
|
||||||
|
else
|
||||||
|
info = 0;
|
||||||
end;
|
end;
|
||||||
return;
|
return;
|
||||||
|
|
||||||
|
|
|
@ -115,13 +115,20 @@ elseif options_.block && ~options_.bytecode
|
||||||
n = size(M_.blocksMFS{b}, 1);
|
n = size(M_.blocksMFS{b}, 1);
|
||||||
ss = oo_.steady_state;
|
ss = oo_.steady_state;
|
||||||
if n ~= 0
|
if n ~= 0
|
||||||
[y, check] = dynare_solve('block_mfs_steadystate', ...
|
if options_.solve_algo <= 4
|
||||||
ss(M_.blocksMFS{b}), ...
|
[y, check] = dynare_solve('block_mfs_steadystate', ...
|
||||||
options_.jacobian_flag, b);
|
ss(M_.blocksMFS{b}), ...
|
||||||
if check ~= 0
|
options_.jacobian_flag, b);
|
||||||
error(['STEADY: convergence problems in block ' int2str(b)])
|
if check ~= 0
|
||||||
|
error(['STEADY: convergence problems in block ' int2str(b)])
|
||||||
|
end
|
||||||
|
ss(M_.blocksMFS{b}) = y;
|
||||||
|
else
|
||||||
|
[ss, check] = solve_one_boundary([M_.fname '_static_' int2str(b)], ss, [oo_.exo_steady_state; ...
|
||||||
|
oo_.exo_det_steady_state], M_.params, M_.blocksMFS{b}, n, 1, 0, b, 0, options_.maxit_, ...
|
||||||
|
options_.solve_tolf, options_.slowc, 0, options_.solve_algo, 1, 0, 0);
|
||||||
|
|
||||||
end
|
end
|
||||||
ss(M_.blocksMFS{b}) = y;
|
|
||||||
end
|
end
|
||||||
[r, g1, oo_.steady_state] = feval([M_.fname '_static'], b, ss, ...
|
[r, g1, oo_.steady_state] = feval([M_.fname '_static'], b, ss, ...
|
||||||
[oo_.exo_steady_state; ...
|
[oo_.exo_steady_state; ...
|
||||||
|
|
|
@ -329,7 +329,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 8 && steady_state))
|
if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 5 && steady_state))
|
||||||
{
|
{
|
||||||
for (i = 0; i < u_count_init; i++)
|
for (i = 0; i < u_count_init; i++)
|
||||||
{
|
{
|
||||||
|
@ -340,7 +340,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
|
||||||
IM_i[make_pair(make_pair(eq, var), lag)] = j;
|
IM_i[make_pair(make_pair(eq, var), lag)] = j;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
else if ( ((stack_solve_algo >= 0 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 5 || solve_algo <= 7) && steady_state) )
|
else if ( ((stack_solve_algo >= 0 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 6 || solve_algo <= 8) && steady_state) )
|
||||||
{
|
{
|
||||||
for (i = 0; i < u_count_init; i++)
|
for (i = 0; i < u_count_init; i++)
|
||||||
{
|
{
|
||||||
|
@ -1949,7 +1949,7 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double sl
|
||||||
#ifdef OCTAVE_MEX_FILE
|
#ifdef OCTAVE_MEX_FILE
|
||||||
ostringstream tmp;
|
ostringstream tmp;
|
||||||
if (steady_state)
|
if (steady_state)
|
||||||
tmp << " GMRES method is not implemented in Octave. You cannot use solve_algo=6, change solve_algo.\n";
|
tmp << " GMRES method is not implemented in Octave. You cannot use solve_algo=7, change solve_algo.\n";
|
||||||
else
|
else
|
||||||
tmp << " GMRES method is not implemented in Octave. You cannot use stack_solve_algo=2, change stack_solve_algo.\n";
|
tmp << " GMRES method is not implemented in Octave. You cannot use stack_solve_algo=2, change stack_solve_algo.\n";
|
||||||
throw FatalExceptionHandling(tmp.str());
|
throw FatalExceptionHandling(tmp.str());
|
||||||
|
@ -2956,19 +2956,19 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
|
||||||
mexPrintf("MODEL STEADY STATE: MATLAB csolve\n");
|
mexPrintf("MODEL STEADY STATE: MATLAB csolve\n");
|
||||||
break;
|
break;
|
||||||
case 5:
|
case 5:
|
||||||
mexPrintf("MODEL STEADY STATE: Sparse LU\n");
|
mexPrintf("MODEL STEADY STATE: (method=ByteCode own solver)\n");
|
||||||
break;
|
break;
|
||||||
case 6:
|
case 6:
|
||||||
mexPrintf("MODEL SIMULATION: (method=GMRES)\n");
|
mexPrintf("MODEL STEADY STATE: Sparse LU\n");
|
||||||
break;
|
break;
|
||||||
case 7:
|
case 7:
|
||||||
mexPrintf("MODEL SIMULATION: (method=BiCGStab)\n");
|
mexPrintf("MODEL STEADY STATE: (method=GMRES)\n");
|
||||||
break;
|
break;
|
||||||
case 8:
|
case 8:
|
||||||
mexPrintf("MODEL SIMULATION: (method=ByteCode own solver)\n");
|
mexPrintf("MODEL STEADY STATE: (method=BiCGStab)\n");
|
||||||
break;
|
break;
|
||||||
default:
|
default:
|
||||||
mexPrintf("MODEL SIMULATION: (method=Unknown - %d - )\n", stack_solve_algo);
|
mexPrintf("MODEL STEADY STATE: (method=Unknown - %d - )\n", stack_solve_algo);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -2980,7 +2980,7 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
|
||||||
mexPrintf("-----------------------------------\n");
|
mexPrintf("-----------------------------------\n");
|
||||||
}
|
}
|
||||||
bool zero_solution;
|
bool zero_solution;
|
||||||
if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 5 && !steady_state))
|
if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state))
|
||||||
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i, zero_solution);
|
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i, zero_solution);
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
|
@ -3012,13 +3012,13 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 5 && !steady_state))
|
if ((solve_algo == 5 && steady_state) || (stack_solve_algo == 5 && !steady_state))
|
||||||
Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
|
Solve_ByteCode_Sparse_GaussianElimination(Size, blck, steady_state, it_);
|
||||||
else if ((solve_algo == 6 && steady_state) || (stack_solve_algo == 2 && !steady_state))
|
else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 2 && !steady_state))
|
||||||
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_, steady_state);
|
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck, false, it_, steady_state);
|
||||||
else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 3 && !steady_state))
|
else if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 3 && !steady_state))
|
||||||
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_);
|
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck, false, it_);
|
||||||
else if ((solve_algo == 5 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1) && !steady_state))
|
else if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1) && !steady_state))
|
||||||
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
|
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc, false, it_);
|
||||||
}
|
}
|
||||||
return;
|
return;
|
||||||
|
|
Loading…
Reference in New Issue