00001
00010 #ifndef BS_JACOBIAN_IMPL_H_
00011 #define BS_JACOBIAN_IMPL_H_
00012 #include "solve_helper.h"
00013
00014 namespace blue_sky {
00015
00021 template <typename strategy_t>
00022 struct jacobian_impl
00023 {
00024 typedef typename strategy_t::item_t item_t;
00025 typedef typename strategy_t::index_t index_t;
00026 typedef typename strategy_t::matrix_t matrix_t;
00027 typedef typename strategy_t::csr_matrix_t bcsr_matrix_t;
00028 typedef typename strategy_t::item_array_t item_array_t;
00029 typedef typename strategy_t::rhs_item_array_t rhs_item_array_t;
00030
00031 typedef jacobian <strategy_t> jacobian_t;
00032 typedef jacobian_matrix <strategy_t> jmatrix_t;
00033 typedef linear_solver_base <strategy_t> linear_solver_t;
00034
00035 typedef smart_ptr <jacobian_t, true> sp_jacobian_t;
00036 typedef smart_ptr <jmatrix_t, true> sp_jmatrix_t;
00037 typedef smart_ptr <linear_solver_t, true> sp_linear_solver_t;
00038 typedef smart_ptr <matrix_t> sp_matrix_t;
00039 typedef smart_ptr <bcsr_matrix_t, true> sp_bcsr_matrix_t;
00040
00041 typedef sp_linear_solver_t sp_solver_t;
00042
00043 public:
00049 jacobian_impl (sp_jacobian_t &jacobian, sp_jmatrix_t &jmatrix)
00050 : jacobian_ (jacobian),
00051 jmatrix_ (jmatrix),
00052 solver_ (jacobian_->get_solver ()),
00053 preconditioner_ (jacobian_->get_prec ()),
00054 rhs_ (jmatrix_->get_rhs ()),
00055 sol_ (jmatrix_->get_solution ()),
00056 regular_matrix_ (jmatrix_->get_regular_matrix ()),
00057 irregular_matrix_ (jmatrix_->get_irregular_matrix ())
00058 {
00059 if (!jmatrix_)
00060 {
00061 bs_throw_exception ("Jacobian matrix is not inited!");
00062 }
00063 if (!solver_)
00064 {
00065 bs_throw_exception ("Solver is not inited!");
00066 }
00067 if (!preconditioner_)
00068 {
00069 bs_throw_exception ("Preconditioner is not inited!");
00070 }
00071 }
00072
00076 void
00077 setup_jacobian ()
00078 {
00079 full_matrix_print ();
00080
00081 jmatrix_->prepare_matrix ();
00082
00083 OMP_TIME_MEASURE_START(gmres_setup_timer);
00084 if (0 != solver_->setup (jmatrix_))
00085 {
00086 bs_throw_exception ("Internal error: cannot build preconditioner for linear solver.");
00087 }
00088 OMP_TIME_MEASURE_END (gmres_setup_timer);
00089 }
00090
00098 item_t
00099 solve_jacobian (index_t &n_lin_iters)
00100 {
00101 #ifdef _MPI
00102 BS_ASSERT (false && "MPI: NOT IMPL YET");
00104 #else //_MPI
00105 sp_matrix_t working_matrix;
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117 working_matrix = jmatrix_->get_prepared_matrix ();
00118 const sp_matrix_t &locked_working_mx (working_matrix);
00119
00120 debug::print_memory_info ("-> jacobian_solver_solve");
00121 index_t ret_code = solve_helper (&*solver_, &(*locked_working_mx), rhs_, sol_);
00122 debug::print_memory_info ("<- jacobian_solver_solve");
00123 #endif //_MPI
00124
00125 jmatrix_->restore_sec_solution ();
00126
00127 static size_t iter_counter = 0;
00128 ++iter_counter;
00131
00132
00133
00134
00135
00136 if (ret_code < 0)
00137 {
00138 BOSERR (section::solvers, level::error) << "Linear solver failed with retcode = " << ret_code << bs_end;
00139 return 10.0f;
00140 }
00141
00142 #ifdef __FULL_MATRIX_PRINT__
00143 ++cur_iter;
00144 #endif //__FULL_MATRIX_PRINT__
00145
00146 n_lin_iters = solver_->get_prop()->get_iters();
00147 item_t tol = solver_->get_prop()->get_final_resid();
00148 if (tol > solver_->get_prop()->get_tolerance())
00149 {
00150 BOSERR (section::solvers, level::error)
00151 << "Linear solver failed with tolerance "
00152 << tol
00153 << bs_end;
00154 }
00155 else
00156 {
00157 BOSOUT (section::solvers, level::medium)
00158 << "Linear solver iterations "
00159 << n_lin_iters
00160 << ", tol = " << tol
00161 << ", ret_code = " << ret_code
00162 << bs_end;
00163 }
00164
00165 return tol;
00166 }
00167
00172 void
00173 full_matrix_print ()
00174 {
00175 #ifdef __FULL_MATRIX_PRINT__
00176 if (jacobian_->cur_iter % jacobian_->sav_iters == 0)
00177 {
00178 char reg_filename [50];
00179 char irr_filename [50];
00180
00181 sprintf (reg_filename, "regular_matrix%d.bs.b", jacobian_->cur_iter);
00182 sprintf (irr_filename, "irregular_matrix%d.bs.csr", jacobian_->cur_iter);
00183
00184 regular_matrix_->ascii_write_in_csr_format(reg_filename);
00185 irregular_matrix_->ascii_write_in_csr_format(irr_filename);
00186 }
00187 #endif
00188 }
00189
00190 private:
00191
00192 sp_jacobian_t &jacobian_;
00193 sp_jmatrix_t &jmatrix_;
00194 sp_solver_t solver_;
00195 sp_solver_t preconditioner_;
00196
00197 rhs_item_array_t &rhs_;
00198 item_array_t &sol_;
00199
00200 sp_matrix_t regular_matrix_;
00201 sp_matrix_t irregular_matrix_;
00202 };
00203
00204
00205 }
00206
00207
00208 #endif // #ifndef BS_JACOBIAN_IMPL_H_
00209