From 39718c86456246f0950080201e6ca910f47ba176 Mon Sep 17 00:00:00 2001 From: Ferhat Mihoubi Date: Fri, 22 Oct 2010 16:49:47 +0200 Subject: [PATCH] - 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: BiCGStab --- matlab/solve_one_boundary.m | 33 +++++++++++++++++++++++----- matlab/steady_.m | 19 +++++++++++----- mex/sources/bytecode/SparseMatrix.cc | 26 +++++++++++----------- 3 files changed, 54 insertions(+), 24 deletions(-) diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index 48e1f579b..1fc231e05 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -97,7 +97,7 @@ for it_=start:incr:finish if(is_dynamic) [r, y, g1, g2, g3] = feval(fname, y, x, params, it_, 0); else - [r, y, g1, g2, g3] = feval(fname, y, x, params, 0); + [r, y, g1] = feval(fname, y, x, params); end; if(~isreal(r)) max_res=(-(max(max(abs(r))))^2)^0.5; @@ -219,6 +219,9 @@ for it_=start:incr:finish ya_save=ya; g1a=g1; if(~is_dynamic & options_.solve_algo == 0) + if (verbose == 1) + disp("steady: fsolve"); + end if exist('OCTAVE_VERSION') || isempty(ver('optim')) % 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 @@ -238,6 +241,9 @@ for it_=start:incr:finish info = -Block_Num*10; end 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; stpmx = 100 ; if (is_dynamic) @@ -261,10 +267,16 @@ for it_=start:incr:finish y = ya'; end; 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); dx = ya - 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; ya = ya - lambda*dx; if(is_dynamic) @@ -272,8 +284,14 @@ for it_=start:incr:finish else y(y_index_eq) = ya; end; - elseif(stack_solve_algo==2 & is_dynamic), + elseif((stack_solve_algo==2 & is_dynamic) | (options_.solve_algo==7 & ~is_dynamic)), 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) [L1, U1]=luinc(g1,luinc_tol); [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; - elseif(stack_solve_algo==3 & is_dynamic), + elseif((stack_solve_algo==3 & is_dynamic) | (options_.solve_algo==8 & ~is_dynamic)), flag1=1; + if (verbose == 1 & ~is_dynamic) + disp("steady: BiCGStab"); + end while(flag1>0) [L1, U1]=luinc(g1,luinc_tol); [dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1); @@ -352,14 +373,16 @@ for it_=start:incr:finish return; end end -info = 1; if(is_dynamic) + info = 1; oo_.deterministic_simulation.status = 1; oo_.deterministic_simulation.error = max_res; oo_.deterministic_simulation.iterations = iter; 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; +else + info = 0; end; return; diff --git a/matlab/steady_.m b/matlab/steady_.m index 4b7261122..27e2780ca 100644 --- a/matlab/steady_.m +++ b/matlab/steady_.m @@ -115,13 +115,20 @@ elseif options_.block && ~options_.bytecode n = size(M_.blocksMFS{b}, 1); ss = oo_.steady_state; if n ~= 0 - [y, check] = dynare_solve('block_mfs_steadystate', ... - ss(M_.blocksMFS{b}), ... - options_.jacobian_flag, b); - if check ~= 0 - error(['STEADY: convergence problems in block ' int2str(b)]) + if options_.solve_algo <= 4 + [y, check] = dynare_solve('block_mfs_steadystate', ... + ss(M_.blocksMFS{b}), ... + options_.jacobian_flag, 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 - ss(M_.blocksMFS{b}) = y; end [r, g1, oo_.steady_state] = feval([M_.fname '_static'], b, ss, ... [oo_.exo_steady_state; ... diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 256e36ae1..c1d2ee6a0 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -329,7 +329,7 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i } 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++) { @@ -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; } } - 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++) { @@ -1949,7 +1949,7 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double sl #ifdef OCTAVE_MEX_FILE ostringstream tmp; 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 tmp << " GMRES method is not implemented in Octave. You cannot use stack_solve_algo=2, change stack_solve_algo.\n"; 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"); break; case 5: - mexPrintf("MODEL STEADY STATE: Sparse LU\n"); + mexPrintf("MODEL STEADY STATE: (method=ByteCode own solver)\n"); break; case 6: - mexPrintf("MODEL SIMULATION: (method=GMRES)\n"); + mexPrintf("MODEL STEADY STATE: Sparse LU\n"); break; case 7: - mexPrintf("MODEL SIMULATION: (method=BiCGStab)\n"); + mexPrintf("MODEL STEADY STATE: (method=GMRES)\n"); break; case 8: - mexPrintf("MODEL SIMULATION: (method=ByteCode own solver)\n"); + mexPrintf("MODEL STEADY STATE: (method=BiCGStab)\n"); break; 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"); } 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); else { @@ -3012,13 +3012,13 @@ SparseMatrix::Simulate_Newton_One_Boundary(int blck, int y_size, int it_, int y_ } 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_); - 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); - 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_); - 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_); } return;