31#include <oomph-lib-config.h>
62 template<
typename MATRIX>
71 throw OomphLibError(
"No matrix was stored -- cannot re-solve",
72 OOMPH_CURRENT_FUNCTION,
73 OOMPH_EXCEPTION_LOCATION);
78 this->solve(Matrix_pt, rhs, result);
89 template<
typename MATRIX>
96 clock_t t_start = clock();
107 Matrix_pt =
new MATRIX;
113 Matrix_can_be_deleted =
true;
118 Jacobian_setup_time = t_end - t_start;
120 clock_t t_end = clock();
121 Jacobian_setup_time = double(t_end - t_start) / CLOCKS_PER_SEC;
126 oomph_info <<
"Time for setup of Jacobian [sec]: " << Jacobian_setup_time
134 this->build_distribution(
136 ->distribution_pt());
145 if ((result.
built()) &&
149 result.
build(this->distribution_pt(), 0.0);
150 this->solve_helper(Matrix_pt, f, result);
155 this->solve_helper(Matrix_pt, f, result);
170 template<
typename MATRIX>
179 std::ostringstream error_message_stream;
180 error_message_stream <<
"The vectors rhs must be setup";
182 OOMPH_CURRENT_FUNCTION,
183 OOMPH_EXCEPTION_LOCATION);
187 if (matrix_pt->
nrow() != matrix_pt->
ncol())
189 std::ostringstream error_message_stream;
190 error_message_stream <<
"The matrix at matrix_pt must be square.";
192 OOMPH_CURRENT_FUNCTION,
193 OOMPH_EXCEPTION_LOCATION);
197 if (matrix_pt->
nrow() != rhs.
nrow())
199 std::ostringstream error_message_stream;
201 <<
"The matrix and the rhs vector must have the same number of rows.";
203 OOMPH_CURRENT_FUNCTION,
204 OOMPH_EXCEPTION_LOCATION);
211 if (dist_matrix_pt != 0)
215 std::ostringstream error_message_stream;
217 <<
"The matrix matrix_pt must have the same communicator as the "
219 <<
" rhs and result must have the same communicator";
221 OOMPH_CURRENT_FUNCTION,
222 OOMPH_EXCEPTION_LOCATION);
231 std::ostringstream error_message_stream;
233 <<
"The matrix (matrix_pt) is not distributable and therefore the rhs"
234 <<
" vector must not be distributed";
236 OOMPH_CURRENT_FUNCTION,
237 OOMPH_EXCEPTION_LOCATION);
242 if (solution.
built())
246 std::ostringstream error_message_stream;
247 error_message_stream <<
"The solution vector distribution has been "
248 "setup; it must have the "
249 <<
"same distribution as the rhs vector.";
251 OOMPH_CURRENT_FUNCTION,
252 OOMPH_EXCEPTION_LOCATION);
260 solution.
build(this->distribution_pt(), 0.0);
270 unsigned nrow_local = this->nrow_local();
276 clock_t t_start = clock();
282 double residual_norm = residual.
norm();
283 double rhs_norm = residual_norm;
284 if (rhs_norm == 0.0) rhs_norm = 1.0;
291 double normalised_residual_norm = residual_norm / rhs_norm;
295 if (Doc_convergence_history)
297 if (!Output_file_stream.is_open())
299 oomph_info << 0 <<
" " << normalised_residual_norm << std::endl;
303 Output_file_stream << 0 <<
" " << normalised_residual_norm << std::endl;
308 if (normalised_residual_norm < Tolerance)
312 oomph_info <<
"BiCGStab converged immediately" << std::endl;
318 Solution_time = t_end - t_start;
322 oomph_info <<
"Time for solve with BiCGStab [sec]: " << Solution_time
332 if (Setup_preconditioner_before_solve)
337 preconditioner_pt()->setup(matrix_pt);
341 Preconditioner_setup_time = t_end_prec - t_start_prec;
345 oomph_info <<
"Time for setup of preconditioner [sec]: "
346 << Preconditioner_setup_time << std::endl;
354 oomph_info <<
"Setup of preconditioner is bypassed in resolve mode"
360 double rho_prev = 0.0;
363 double rho, beta, dot_prod, dot_prod_tt, dot_prod_ts;
364 double s_norm, r_norm;
368 p_hat(this->distribution_pt(), 0.0), v(this->distribution_pt(), 0.0),
369 z(this->distribution_pt(), 0.0),
t(this->distribution_pt(), 0.0),
370 s(this->distribution_pt(), 0.0);
373 for (
unsigned iter = 1; iter <=
Max_iter; iter++)
376 rho = r_hat.
dot(residual);
381 oomph_info <<
"BiCGStab has broken down after " << iter <<
" iterations"
383 oomph_info <<
"Returning with current normalised residual of "
384 << normalised_residual_norm << std::endl;
394 beta = (rho / rho_prev) * (alpha / omega);
395 for (
unsigned i = 0;
i < nrow_local;
i++)
397 p[
i] = residual[
i] + beta * (p[
i] - omega * v[
i]);
402 preconditioner_pt()->preconditioner_solve(p, p_hat);
406 dot_prod = r_hat.
dot(v);
407 alpha = rho / dot_prod;
408 for (
unsigned i = 0;
i < nrow_local;
i++)
410 s[
i] = residual[
i] - alpha * v[
i];
415 normalised_residual_norm = s_norm / rhs_norm;
419 if (Doc_convergence_history)
421 if (!Output_file_stream.is_open())
423 oomph_info << double(iter - 0.5) <<
" " << normalised_residual_norm
428 Output_file_stream << double(iter - 0.5) <<
" "
429 << normalised_residual_norm << std::endl;
434 if (normalised_residual_norm < Tolerance)
436 for (
unsigned i = 0;
i < nrow_local;
i++)
438 solution[
i] = x[
i] + alpha * p_hat[
i];
444 oomph_info <<
"BiCGStab converged. Normalised residual norm: "
445 << normalised_residual_norm << std::endl;
446 oomph_info <<
"Number of iterations to convergence: " << iter
456 Solution_time = t_end - t_start;
460 oomph_info <<
"Time for solve with BiCGStab [sec]: " << Solution_time
468 preconditioner_pt()->preconditioner_solve(
s, z);
472 dot_prod_ts =
t.dot(
s);
473 dot_prod_tt =
t.dot(
t);
474 omega = dot_prod_ts / dot_prod_tt;
475 for (
unsigned i = 0;
i < nrow_local;
i++)
477 x[
i] += alpha * p_hat[
i] + omega * z[
i];
478 residual[
i] =
s[
i] - omega *
t[
i];
480 r_norm = residual.
norm();
484 normalised_residual_norm = r_norm / rhs_norm;
488 if (Doc_convergence_history)
490 if (!Output_file_stream.is_open())
492 oomph_info << iter <<
" " << normalised_residual_norm << std::endl;
496 Output_file_stream << iter <<
" " << normalised_residual_norm
502 if (normalised_residual_norm < Tolerance)
507 oomph_info <<
"BiCGStab converged. Normalised residual norm: "
508 << normalised_residual_norm << std::endl;
509 oomph_info <<
"Number of iterations to convergence: " << iter
520 Solution_time = t_end - t_start;
524 oomph_info <<
"Time for solve with BiCGStab [sec]: " << Solution_time
535 oomph_info <<
"BiCGStab breakdown with omega=0.0. "
536 <<
"Normalised residual norm: " << normalised_residual_norm
538 oomph_info <<
"Number of iterations so far: " << iter << std::endl;
547 Solution_time = t_end - t_start;
551 oomph_info <<
"Time for solve with BiCGStab [sec]: " << Solution_time
563 oomph_info <<
"BiCGStab did not converge to required tolerance! "
565 oomph_info <<
"Returning with normalised residual norm: "
566 << normalised_residual_norm << std::endl;
574 Solution_time = t_end - t_start;
578 oomph_info <<
"Time for solve with BiCGStab [sec]: " << Solution_time
582 if (Throw_error_after_max_iter)
584 std::string err =
"Solver failed to converge and you requested an error";
585 err +=
" on convergence failures.";
587 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
604 template<
typename MATRIX>
613 std::ostringstream error_message_stream;
614 error_message_stream <<
"The vectors rhs must be setup";
616 OOMPH_CURRENT_FUNCTION,
617 OOMPH_EXCEPTION_LOCATION);
621 if (matrix_pt->
nrow() != matrix_pt->
ncol())
623 std::ostringstream error_message_stream;
624 error_message_stream <<
"The matrix at matrix_pt must be square.";
626 OOMPH_CURRENT_FUNCTION,
627 OOMPH_EXCEPTION_LOCATION);
631 if (matrix_pt->
nrow() != rhs.
nrow())
633 std::ostringstream error_message_stream;
635 <<
"The matrix and the rhs vector must have the same number of rows.";
637 OOMPH_CURRENT_FUNCTION,
638 OOMPH_EXCEPTION_LOCATION);
645 if (dist_matrix_pt != 0)
649 std::ostringstream error_message_stream;
651 <<
"The matrix matrix_pt must have the same communicator as the "
653 <<
" rhs and result must have the same communicator";
655 OOMPH_CURRENT_FUNCTION,
656 OOMPH_EXCEPTION_LOCATION);
665 std::ostringstream error_message_stream;
667 <<
"The matrix (matrix_pt) is not distributable and therefore the rhs"
668 <<
" vector must not be distributed";
670 OOMPH_CURRENT_FUNCTION,
671 OOMPH_EXCEPTION_LOCATION);
676 if (solution.
built())
680 std::ostringstream error_message_stream;
681 error_message_stream <<
"The solution vector distribution has been "
682 "setup; it must have the "
683 <<
"same distribution as the rhs vector.";
685 OOMPH_CURRENT_FUNCTION,
686 OOMPH_EXCEPTION_LOCATION);
694 solution.
build(this->distribution_pt(), 0.0);
704 unsigned nrow_local = this->nrow_local();
707 unsigned counter = 0;
716 double residual_norm = residual.
norm();
717 double rhs_norm = residual_norm;
718 if (rhs_norm == 0.0) rhs_norm = 1.0;
721 double normalised_residual_norm = residual_norm / rhs_norm;
725 if (Doc_convergence_history)
727 if (!Output_file_stream.is_open())
729 oomph_info << 0 <<
" " << normalised_residual_norm << std::endl;
733 Output_file_stream << 0 <<
" " << normalised_residual_norm << std::endl;
738 if (normalised_residual_norm < Tolerance)
742 oomph_info <<
"CG converged immediately" << std::endl;
748 Solution_time = t_end - t_start;
752 oomph_info <<
"Time for solve with CG [sec]: " << Solution_time
763 if (Setup_preconditioner_before_solve)
768 preconditioner_pt()->setup(matrix_pt);
772 Preconditioner_setup_time = t_end_prec - t_start_prec;
776 oomph_info <<
"Time for setup of preconditioner [sec]: "
777 << Preconditioner_setup_time << std::endl;
785 oomph_info <<
"Setup of preconditioner is bypassed in resolve mode"
794 p(this->distribution_pt(), 0.0),
795 jacobian_times_p(this->distribution_pt(), 0.0);
798 double alpha, beta, rz;
799 double prev_rz = 0.0;
802 while ((normalised_residual_norm > Tolerance) && (counter !=
Max_iter))
805 preconditioner_pt()->preconditioner_solve(residual, z);
811 rz = residual.
dot(z);
816 rz = residual.
dot(z);
818 for (
unsigned i = 0;
i < nrow_local;
i++)
820 p[
i] = z[
i] + beta * p[
i];
826 matrix_pt->
multiply(p, jacobian_times_p);
827 double pq = p.dot(jacobian_times_p);
832 for (
unsigned i = 0;
i < nrow_local;
i++)
834 x[
i] += alpha * p[
i];
835 residual[
i] -= alpha * jacobian_times_p[
i];
840 residual_norm = residual.
norm();
843 normalised_residual_norm = residual_norm / rhs_norm;
848 if (Doc_convergence_history)
850 if (!Output_file_stream.is_open())
852 oomph_info << counter <<
" " << normalised_residual_norm << std::endl;
856 Output_file_stream << counter <<
" " << normalised_residual_norm
861 counter = counter + 1;
869 oomph_info <<
"CG did not converge to required tolerance! " << std::endl;
870 oomph_info <<
"Returning with normalised residual norm: "
871 << normalised_residual_norm << std::endl;
872 oomph_info <<
"after " << counter <<
" iterations." << std::endl;
880 oomph_info <<
"CG converged. Normalised residual norm: "
881 << normalised_residual_norm << std::endl;
882 oomph_info <<
"Number of iterations to convergence: " << counter
890 Iterations = counter;
897 Solution_time = t_end - t_start;
901 oomph_info <<
"Time for solve with CG [sec]: " << Solution_time
905 if ((counter >=
Max_iter) && (Throw_error_after_max_iter))
907 std::string err =
"Solver failed to converge and you requested an error";
908 err +=
" on convergence failures.";
910 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
921 template<
typename MATRIX>
930 throw OomphLibError(
"No matrix was stored -- cannot re-solve",
931 OOMPH_CURRENT_FUNCTION,
932 OOMPH_EXCEPTION_LOCATION);
937 this->solve(Matrix_pt, rhs, result);
949 template<
typename MATRIX>
963 Matrix_pt =
new MATRIX;
968 Matrix_can_be_deleted =
true;
972 Jacobian_setup_time = t_end - t_start;
976 oomph_info <<
"Time for setup of Jacobian [sec]: " << Jacobian_setup_time
984 this->build_distribution(
986 ->distribution_pt());
997 result.
build(this->distribution_pt(), 0.0);
1004 result.
build(this->distribution_pt(), 0.0);
1005 this->solve_helper(Matrix_pt, f, result);
1010 this->solve_helper(Matrix_pt, f, result);
1028 template<
typename MATRIX>
1032 const double& n_dof)
1041 std::ostringstream error_message_stream;
1042 error_message_stream <<
"The matrix must not be distributed.";
1044 OOMPH_CURRENT_FUNCTION,
1045 OOMPH_EXCEPTION_LOCATION);
1051 std::ostringstream error_message_stream;
1052 error_message_stream <<
"The rhs vector distribution must be setup.";
1054 OOMPH_CURRENT_FUNCTION,
1055 OOMPH_EXCEPTION_LOCATION);
1058 if (rhs.
nrow() != n_dof)
1060 std::ostringstream error_message_stream;
1061 error_message_stream <<
"RHS does not have the same dimension as "
1062 <<
"the linear system";
1064 OOMPH_CURRENT_FUNCTION,
1065 OOMPH_EXCEPTION_LOCATION);
1070 std::ostringstream error_message_stream;
1071 error_message_stream <<
"The rhs vector must not be distributed.";
1073 OOMPH_CURRENT_FUNCTION,
1074 OOMPH_EXCEPTION_LOCATION);
1078 if (solution.
built())
1082 std::ostringstream error_message_stream;
1083 error_message_stream <<
"If the result distribution is setup then it "
1084 <<
"must be the same as the rhs distribution";
1086 OOMPH_CURRENT_FUNCTION,
1087 OOMPH_EXCEPTION_LOCATION);
1103 template<
typename MATRIX>
1108 Use_as_smoother =
false;
1111 unsigned n_dof = problem_pt->
ndof();
1126 this->build_distribution(dist);
1130 Matrix_pt =
new MATRIX;
1139 Matrix_can_be_deleted =
true;
1143 Jacobian_setup_time = t_end - t_start;
1148 oomph_info <<
"Time for setup of Jacobian [sec]: " << Jacobian_setup_time
1153 this->solve_helper(Matrix_pt, f, result);
1163 template<
typename MATRIX>
1169 unsigned n_dof = rhs.
nrow();
1173 MATRIX* tmp_matrix_pt =
dynamic_cast<MATRIX*
>(matrix_pt);
1176 this->check_validity_of_solve_helper_inputs<MATRIX>(
1177 tmp_matrix_pt, rhs, solution, n_dof);
1188 solution.
build(this->distribution_pt(), 0.0);
1197 if (!Use_as_smoother)
1216 double norm_res = 0.0;
1220 double norm_f = 0.0;
1226 if (!Use_as_smoother)
1229 current_residual.
build(this->distribution_pt(), 0.0);
1232 matrix_pt->
residual(x, rhs, current_residual);
1235 norm_res = current_residual.
norm();
1242 if (Doc_convergence_history)
1245 if (!Output_file_stream.is_open())
1248 oomph_info << Iterations <<
" " << norm_res << std::endl;
1254 Output_file_stream << Iterations <<
" " << norm_res << std::endl;
1260 for (
unsigned iter_num = 0; iter_num <
Max_iter; iter_num++)
1263 for (
unsigned i = 0;
i < n_dof;
i++)
1265 double dummy = rhs[
i];
1266 for (
unsigned j = 0; j <
i; j++)
1268 dummy -= (*(matrix_pt))(
i, j) * x[j];
1270 for (
unsigned j = (
i + 1); j < n_dof; j++)
1272 dummy -= (*(matrix_pt))(
i, j) * x[j];
1274 x[
i] = dummy / (*(matrix_pt))(
i,
i);
1281 if (!Use_as_smoother)
1284 matrix_pt->
residual(x, rhs, current_residual);
1288 norm_res = current_residual.
norm() / norm_f;
1292 if (Doc_convergence_history)
1294 if (!Output_file_stream.is_open())
1296 oomph_info << Iterations <<
" " << norm_res << std::endl;
1300 Output_file_stream << Iterations <<
" " << norm_res << std::endl;
1305 if (norm_res < Tolerance)
1314 if (!Use_as_smoother)
1319 oomph_info <<
"\nGS converged. Residual norm: " << norm_res
1320 <<
"\nNumber of iterations to convergence: " << Iterations
1331 Solution_time = t_end - t_start;
1334 oomph_info <<
"Time for solve with GS [sec]: " << Solution_time
1340 if ((Iterations >
Max_iter - 1) && (Throw_error_after_max_iter))
1343 "Solver failed to converge and you requested ";
1344 error_message +=
"an error on convergence failures.";
1346 error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1370 unsigned n_dof = problem_pt->
ndof();
1453 Index_of_diagonal_entries =
Matrix_pt->get_index_of_diagonal_entries();
1458 bool is_matrix_sorted =
true;
1461 is_matrix_sorted =
Matrix_pt->entries_are_sorted();
1464 if (!is_matrix_sorted)
1469 OOMPH_CURRENT_FUNCTION,
1470 OOMPH_EXCEPTION_LOCATION);
1485 unsigned n_dof = rhs.
nrow();
1492 this->check_validity_of_solve_helper_inputs<CRDoubleMatrix>(
1493 cr_matrix_pt, rhs, solution, n_dof);
1531 double norm_res = 0.0;
1535 double norm_f = 0.0;
1547 matrix_pt->
residual(x, rhs, current_residual);
1550 norm_res = current_residual.
norm();
1584 const double* value_pt = tmp_matrix_pt->
value();
1585 const int* row_start_pt = tmp_matrix_pt->
row_start();
1586 const int* column_index_pt = tmp_matrix_pt->
column_index();
1594 for (
unsigned i = 0;
i < n_dof;
i++)
1597 unsigned diag_index = Index_of_diagonal_entries[
i];
1599 if (
unsigned(*(column_index_pt + diag_index)) !=
i)
1601 std::string err_strng =
"Gauss-Seidel cannot operate on a matrix with ";
1602 err_strng +=
"zero diagonal entries.";
1604 err_strng, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1615 constant_term = rhs;
1618 for (
unsigned i = 0;
i < n_dof;
i++)
1621 unsigned diag_index = Index_of_diagonal_entries[
i];
1624 unsigned row_i_start = *(row_start_pt +
i);
1631 if (
unsigned(*(column_index_pt + row_i_start)) !=
i)
1634 unsigned column_index = 0;
1637 for (
unsigned j = row_i_start; j < diag_index; j++)
1640 column_index = *(column_index_pt + j);
1643 constant_term[
i] -= (*(value_pt + j)) * constant_term[column_index];
1648 constant_term[
i] /= *(value_pt + diag_index);
1655 for (
unsigned iter_num = 0; iter_num <
Max_iter; iter_num++)
1672 unsigned upper_tri_start = 0;
1676 unsigned next_row_start = 0;
1684 for (
unsigned i = 0;
i < n_dof - 1;
i++)
1689 upper_tri_start = Index_of_diagonal_entries[
i] + 1;
1692 next_row_start = *(row_start_pt +
i + 1);
1695 unsigned column_index = 0;
1698 for (
unsigned j = upper_tri_start; j < next_row_start; j++)
1701 column_index = *(column_index_pt + j);
1704 temp_vec[
i] += (*(value_pt + j)) * x[column_index];
1724 for (
unsigned i = 0;
i < n_dof;
i++)
1727 unsigned diag_index = Index_of_diagonal_entries[
i];
1730 unsigned row_i_start = *(row_start_pt +
i);
1735 if (
unsigned(*(column_index_pt + row_i_start)) !=
i)
1738 unsigned column_index = 0;
1741 for (
unsigned j = row_i_start; j < diag_index; j++)
1744 column_index = *(column_index_pt + j);
1747 temp_vec[
i] -= (*(value_pt + j)) * temp_vec[column_index];
1752 temp_vec[
i] /= *(value_pt + diag_index);
1778 matrix_pt->
residual(x, rhs, current_residual);
1782 norm_res = current_residual.
norm() / norm_f;
1790 oomph_info << iter_num <<
" " << norm_res << std::endl;
1813 oomph_info <<
"\nGS converged. Residual norm: " << norm_res
1814 <<
"\nNumber of iterations to convergence: " <<
Iterations
1839 "Solver failed to converge and you requested ";
1840 error_message +=
"an error on convergence failures.";
1842 error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1857 template<
typename MATRIX>
1863 Use_as_smoother =
false;
1866 unsigned n_dof = problem_pt->
ndof();
1881 this->build_distribution(dist);
1885 Matrix_pt =
new MATRIX;
1894 extract_diagonal_entries(Matrix_pt);
1897 Matrix_can_be_deleted =
true;
1901 Jacobian_setup_time = t_end - t_start;
1906 oomph_info <<
"Time for setup of Jacobian [sec]: " << Jacobian_setup_time
1911 solve_helper(Matrix_pt, f, result);
1921 template<
typename MATRIX>
1927 unsigned n_dof = rhs.
nrow();
1931 MATRIX* tmp_matrix_pt =
dynamic_cast<MATRIX*
>(matrix_pt);
1934 this->check_validity_of_solve_helper_inputs<MATRIX>(
1935 tmp_matrix_pt, rhs, solution, n_dof);
1947 solution.
build(this->distribution_pt(), 0.0);
1956 if (!Use_as_smoother)
1971 DoubleVector constant_term(this->distribution_pt(), 0.0);
1974 for (
unsigned i = 0;
i < n_dof;
i++)
1977 constant_term[
i] = Omega * Matrix_diagonal[
i] * rhs[
i];
1986 double norm_res = 0.0;
1990 double norm_f = 0.0;
1996 if (!Use_as_smoother)
1999 local_residual.
build(this->distribution_pt(), 0.0);
2002 matrix_pt->
residual(x, rhs, local_residual);
2005 norm_res = local_residual.
norm();
2012 if (Doc_convergence_history)
2014 if (!Output_file_stream.is_open())
2016 oomph_info << Iterations <<
" " << norm_res << std::endl;
2020 Output_file_stream << Iterations <<
" " << norm_res << std::endl;
2031 for (
unsigned iter_num = 0; iter_num <
Max_iter; iter_num++)
2038 for (
unsigned idof = 0; idof < n_dof; idof++)
2042 temp_vec[idof] *= Omega * Matrix_diagonal[idof];
2053 if (!Use_as_smoother)
2056 matrix_pt->
residual(x, rhs, local_residual);
2059 norm_res = local_residual.
norm() / norm_f;
2063 if (Doc_convergence_history)
2065 if (!Output_file_stream.is_open())
2067 oomph_info << Iterations <<
" " << norm_res << std::endl;
2071 Output_file_stream << Iterations <<
" " << norm_res << std::endl;
2076 if (norm_res < Tolerance)
2085 if (!Use_as_smoother)
2090 oomph_info <<
"\nDamped Jacobi converged. Residual norm: " << norm_res
2091 <<
"\nNumber of iterations to convergence: " << Iterations
2102 Solution_time = t_end - t_start;
2105 oomph_info <<
"Time for solve with damped Jacobi [sec]: " << Solution_time
2111 if ((Iterations >
Max_iter - 1) && (Throw_error_after_max_iter))
2114 "Solver failed to converge and you requested ";
2115 error_message +=
"an error on convergence failures.";
2117 error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
2132 template<
typename MATRIX>
2141 throw OomphLibError(
"No matrix was stored -- cannot re-solve",
2142 OOMPH_CURRENT_FUNCTION,
2143 OOMPH_EXCEPTION_LOCATION);
2148 this->solve(Matrix_pt, rhs, result);
2160 template<
typename MATRIX>
2164 unsigned n_dof = problem_pt->
ndof();
2177 this->build_distribution(dist);
2181 Matrix_pt =
new MATRIX;
2188 this->distribution_pt());
2189 f.
build(this->distribution_pt(), 0.0);
2195 Matrix_can_be_deleted =
true;
2199 Jacobian_setup_time = t_end - t_start;
2203 oomph_info <<
"Time for setup of Jacobian [sec]: " << Jacobian_setup_time
2209 if (Compute_gradient)
2212 Matrix_pt->multiply_transpose(f, Gradient_for_glob_conv_newton_solve);
2214 Gradient_has_been_computed =
true;
2221 if ((result.
built()) &&
2225 result.
build(this->distribution_pt(), 0.0);
2226 this->solve_helper(Matrix_pt, f, result);
2232 this->solve_helper(Matrix_pt, f, result);
2249 template<
typename MATRIX>
2255 unsigned n_dof = rhs.
nrow();
2265 std::ostringstream error_message_stream;
2266 error_message_stream <<
"The matrix must not be distributed.";
2268 OOMPH_CURRENT_FUNCTION,
2269 OOMPH_EXCEPTION_LOCATION);
2275 std::ostringstream error_message_stream;
2276 error_message_stream <<
"The rhs vector distribution must be setup.";
2278 OOMPH_CURRENT_FUNCTION,
2279 OOMPH_EXCEPTION_LOCATION);
2282 if (rhs.
nrow() != n_dof)
2285 "RHS does not have the same dimension as the linear system",
2286 OOMPH_CURRENT_FUNCTION,
2287 OOMPH_EXCEPTION_LOCATION);
2292 std::ostringstream error_message_stream;
2293 error_message_stream <<
"The rhs vector must not be distributed.";
2295 OOMPH_CURRENT_FUNCTION,
2296 OOMPH_EXCEPTION_LOCATION);
2300 if (solution.
built())
2304 std::ostringstream error_message_stream;
2305 error_message_stream <<
"If the result distribution is setup then it "
2306 "must be the same as the "
2307 <<
"rhs distribution";
2309 OOMPH_CURRENT_FUNCTION,
2310 OOMPH_EXCEPTION_LOCATION);
2316 Preconditioner_application_time = 0.0;
2319 if (!solution.
built())
2321 solution.
build(this->distribution_pt(), 0.0);
2339 if (!Iteration_restart)
2354 if (Setup_preconditioner_before_solve)
2360 preconditioner_pt()->setup(matrix_pt);
2364 Preconditioner_setup_time = t_end_prec - t_start_prec;
2368 oomph_info <<
"Time for setup of preconditioner [sec]: "
2369 << Preconditioner_setup_time << std::endl;
2377 oomph_info <<
"Setup of preconditioner is bypassed in resolve mode"
2384 if (Preconditioner_LHS)
2390 preconditioner_pt()->preconditioner_solve(rhs, r);
2393 Preconditioner_application_time +=
2404 for (
unsigned i = 0;
i < n_dof;
i++)
2406 normb += r_pt[
i] * r_pt[
i];
2408 normb = sqrt(normb);
2411 double beta = normb;
2414 if (normb == 0.0) normb = 1;
2415 resid = beta / normb;
2419 if (Doc_convergence_history)
2421 if (!Output_file_stream.is_open())
2423 oomph_info << 0 <<
" " << resid << std::endl;
2427 Output_file_stream << 0 <<
" " << resid << std::endl;
2432 if (resid <= Tolerance)
2436 oomph_info <<
"GMRES converged immediately. Normalised residual norm: "
2437 << resid << std::endl;
2441 Solution_time = t_end - t_start;
2446 oomph_info <<
"Time for all preconditioner applications [sec]: "
2447 << Preconditioner_application_time
2448 <<
"\n\nTime for solve with GMRES [sec]: " << Solution_time
2460 v.resize(Restart + 1);
2467 v[0].build(this->distribution_pt(), 0.0);
2468 double* v0_pt = v[0].values_pt();
2469 for (
unsigned i = 0;
i < n_dof;
i++)
2471 v0_pt[
i] = r_pt[
i] / beta;
2478 unsigned iter_restart;
2481 for (iter_restart = 0; iter_restart < Restart && iter <=
Max_iter;
2482 iter_restart++, iter++)
2485 H[iter_restart].resize(iter_restart + 2);
2490 if (Preconditioner_LHS)
2493 matrix_pt->
multiply(v[iter_restart], temp);
2499 preconditioner_pt()->preconditioner_solve(temp, w);
2502 Preconditioner_application_time +=
2511 preconditioner_pt()->preconditioner_solve(v[iter_restart], temp);
2514 Preconditioner_application_time +=
2524 for (
unsigned k = 0; k <= iter_restart; k++)
2527 H[iter_restart][k] = 0.0;
2528 double* vk_pt = v[k].values_pt();
2529 for (
unsigned i = 0;
i < n_dof;
i++)
2531 H[iter_restart][k] += w_pt[
i] * vk_pt[
i];
2535 for (
unsigned i = 0;
i < n_dof;
i++)
2537 w_pt[
i] -= H[iter_restart][k] * vk_pt[
i];
2543 double temp_norm_w = 0.0;
2544 for (
unsigned i = 0;
i < n_dof;
i++)
2546 temp_norm_w += w_pt[
i] * w_pt[
i];
2548 temp_norm_w = sqrt(temp_norm_w);
2549 H[iter_restart][iter_restart + 1] = temp_norm_w;
2553 v[iter_restart + 1].build(this->distribution_pt(), 0.0);
2554 double* v_pt = v[iter_restart + 1].values_pt();
2555 for (
unsigned i = 0;
i < n_dof;
i++)
2557 v_pt[
i] = w_pt[
i] / H[iter_restart][iter_restart + 1];
2561 for (
unsigned k = 0; k < iter_restart; k++)
2563 apply_plane_rotation(
2564 H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
2566 generate_plane_rotation(H[iter_restart][iter_restart],
2567 H[iter_restart][iter_restart + 1],
2570 apply_plane_rotation(H[iter_restart][iter_restart],
2571 H[iter_restart][iter_restart + 1],
2574 apply_plane_rotation(
s[iter_restart],
2575 s[iter_restart + 1],
2580 beta = std::fabs(
s[iter_restart + 1]);
2583 resid = beta / normb;
2587 if (Doc_convergence_history)
2589 if (!Output_file_stream.is_open())
2591 oomph_info << iter <<
" " << resid << std::endl;
2595 Output_file_stream << iter <<
" " << resid << std::endl;
2600 if (resid < Tolerance)
2603 update(iter_restart, H,
s, v, solution);
2609 oomph_info <<
"GMRES converged (1). Normalised residual norm: "
2610 << resid << std::endl;
2611 oomph_info <<
"Number of iterations to convergence: " << iter
2618 Solution_time = t_end - t_start;
2625 oomph_info <<
"Time for all preconditioner applications [sec]: "
2626 << Preconditioner_application_time
2627 <<
"\n\nTime for solve with GMRES [sec]: "
2628 << Solution_time << std::endl;
2635 if (iter_restart > 0) update((iter_restart - 1), H,
s, v, solution);
2640 matrix_pt->
multiply(solution, temp);
2643 for (
unsigned i = 0;
i < n_dof;
i++)
2645 temp_pt[
i] = rhs_pt[
i] - temp_pt[
i];
2648 if (Preconditioner_LHS)
2653 preconditioner_pt()->preconditioner_solve(temp, r);
2656 Preconditioner_application_time +=
2664 for (
unsigned i = 0;
i < n_dof;
i++)
2666 beta += r_pt[
i] * r_pt[
i];
2671 resid = beta / normb;
2672 if (resid < Tolerance)
2677 oomph_info <<
"GMRES converged (2). Normalised residual norm: "
2678 << resid << std::endl;
2679 oomph_info <<
"Number of iterations to convergence: " << iter
2686 Solution_time = t_end - t_start;
2693 oomph_info <<
"Time for all preconditioner applications [sec]: "
2694 << Preconditioner_application_time
2695 <<
"\n\nTime for solve with GMRES [sec]: "
2696 << Solution_time << std::endl;
2708 matrix_pt->
multiply(solution, temp);
2711 resid = temp.
norm() / normb;
2716 oomph_info <<
"GMRES did not converge to required tolerance! " << std::endl;
2717 oomph_info <<
"Returning with normalised residual norm: " << resid
2723 if (Throw_error_after_max_iter)
2725 std::string err =
"Solver failed to converge and you requested an error";
2726 err +=
" on convergence failures.";
2728 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
2750 unsigned n_dof = problem_pt->
ndof();
2807 result.
build(dist, 0.0);
2840 throw OomphLibError(
"No matrix was stored -- cannot re-solve",
2841 OOMPH_CURRENT_FUNCTION,
2842 OOMPH_EXCEPTION_LOCATION);
2868 unsigned n_dof = input_matrix_pt->
nrow();
2871 unsigned n_aug_dof = n_dof + 1;
2883 if ((rhs.
nrow() != n_dof) || (lhs.
nrow() != n_dof))
2886 std::ostringstream error_message_stream;
2889 error_message_stream <<
"RHS vector has " << rhs.
nrow()
2891 <<
"LHS vector has " << rhs.
nrow()
2892 <<
" rows but\nthey "
2893 <<
"should both have " << n_dof <<
" rows!"
2898 OOMPH_CURRENT_FUNCTION,
2899 OOMPH_EXCEPTION_LOCATION);
2910 std::ostringstream error_message_stream;
2911 error_message_stream <<
"The matrix must not be distributed.";
2913 OOMPH_CURRENT_FUNCTION,
2914 OOMPH_EXCEPTION_LOCATION);
2920 std::ostringstream error_message_stream;
2921 error_message_stream <<
"The rhs vector distribution must be setup.";
2923 OOMPH_CURRENT_FUNCTION,
2924 OOMPH_EXCEPTION_LOCATION);
2929 std::ostringstream error_message_stream;
2930 error_message_stream <<
"The rhs vector must not be distributed.";
2932 OOMPH_CURRENT_FUNCTION,
2933 OOMPH_EXCEPTION_LOCATION);
2941 std::ostringstream error_message_stream;
2942 error_message_stream <<
"If the result distribution is setup then it "
2943 "must be the same as the "
2944 <<
"rhs distribution";
2946 OOMPH_CURRENT_FUNCTION,
2947 OOMPH_EXCEPTION_LOCATION);
2957 lhs.
build(dist, 0.0);
3008 oomph_info <<
"Time for setup of preconditioner [sec]: "
3018 oomph_info <<
"Setup of preconditioner is bypassed in resolve mode"
3027 double t_prec_application_total = 0.0;
3044 throw OomphLibError(
"Bordering vectors have not been created!",
3045 OOMPH_CURRENT_FUNCTION,
3046 OOMPH_EXCEPTION_LOCATION);
3052 throw OomphLibError(
"Bordering vectors do not have the right size!",
3053 OOMPH_CURRENT_FUNCTION,
3054 OOMPH_EXCEPTION_LOCATION);
3082 for (
unsigned i = 0;
i < n_dof;
i++)
3085 rhs_aug[
i] = rhs[
i];
3089 rhs_aug[n_dof] = (*Rhs_pt);
3107 t_prec_application_total +=
3108 (t_prec_application_end - t_prec_application_start);
3118 double normb = r.
norm();
3121 double beta = normb;
3131 resid = beta / normb;
3139 oomph_info << 0 <<
" " << resid << std::endl;
3152 oomph_info <<
"AugmentedProblemGMRES converged immediately. Normalised "
3153 <<
"residual norm: " << resid << std::endl;
3163 <<
"Time for all preconditioner applications [sec]: "
3164 << t_prec_application_total
3165 <<
"\nTotal time for solve with AugmentedProblemGMRES [sec]: "
3194 unsigned iter_restart = 0;
3198 iter_restart++, iter++)
3201 H[iter_restart].resize(iter_restart + 2);
3225 t_prec_application_total +=
3226 (t_prec_application_end - t_prec_application_start);
3240 t_prec_application_total +=
3241 (t_prec_application_end - t_prec_application_start);
3252 for (
unsigned k = 0; k <= iter_restart; k++)
3256 H[iter_restart][k] = w.
dot(v[k]);
3259 double* vk_pt = v[k].values_pt();
3262 for (
unsigned i = 0;
i < n_aug_dof;
i++)
3265 w_pt[
i] -= H[iter_restart][k] * vk_pt[
i];
3270 H[iter_restart][iter_restart + 1] = w.
norm();
3273 v[iter_restart + 1] = w;
3276 v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
3279 for (
unsigned k = 0; k < iter_restart; k++)
3283 H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
3287 H[iter_restart][iter_restart + 1],
3291 H[iter_restart][iter_restart + 1],
3295 s[iter_restart + 1],
3300 beta = std::fabs(
s[iter_restart + 1]);
3303 resid = beta / normb;
3315 oomph_info << iter <<
" " << resid << std::endl;
3327 update(iter_restart, H,
s, v, solution);
3330 for (
unsigned i = 0;
i < n_dof;
i++)
3333 lhs[
i] = solution[
i];
3337 (*X_pt) -= solution[n_dof];
3343 oomph_info <<
"\nAugmentedProblemGMRES converged (1). Normalised "
3344 <<
"residual norm: " << resid
3345 <<
"\nNumber of iterations "
3346 <<
"to convergence: " << iter <<
"\n"
3360 <<
"Time for all preconditioner applications [sec]: "
3361 << t_prec_application_total
3362 <<
"\nTotal time for solve with AugmentedProblemGMRES [sec]: "
3375 if (std::fabs(H[iter_restart][iter_restart + 1]) < 1.0e-12)
3378 std::ostringstream warning_message_stream;
3381 warning_message_stream <<
"Hessenberg subdiagonal entry (="
3382 << H[iter_restart][iter_restart + 1]
3383 <<
") almost zero at iteration "
3384 << iter_restart <<
"\nbut residual "
3385 <<
"is only " << resid <<
" which does "
3386 <<
"not meet the tolerance (" <<
Tolerance
3387 <<
")!" << std::endl;
3391 OOMPH_CURRENT_FUNCTION,
3392 OOMPH_EXCEPTION_LOCATION);
3399 if (iter_restart > 0)
3402 update((iter_restart - 1), H,
s, v, solution);
3411 for (
unsigned i = 0;
i < n_dof;
i++)
3414 temp_pt[
i] = rhs_pt[
i] - temp_pt[
i];
3431 t_prec_application_total +=
3432 (t_prec_application_end - t_prec_application_start);
3440 resid = beta / normb;
3446 oomph_info <<
"\nAugmentedProblemGMRES converged (2). Normalised "
3447 <<
"residual norm: " << resid <<
"\nNumber of iterations "
3448 <<
"to convergence: " << iter <<
"\n"
3453 for (
unsigned i = 0;
i < n_dof;
i++)
3456 lhs[
i] = solution[
i];
3460 (*X_pt) -= solution[n_dof];
3472 <<
"Time for all preconditioner applications [sec]: "
3473 << t_prec_application_total
3474 <<
"\nTotal time for solve with AugmentedProblemGMRES [sec]: "
3488 resid = temp.
norm() / normb;
3492 for (
unsigned i = 0;
i < n_dof;
i++)
3495 lhs[
i] = solution[
i];
3499 (*X_pt) -= solution[n_dof];
3502 oomph_info <<
"\nAugmentedProblemGMRES did not converge to required "
3503 <<
"tolerance!\nReturning with normalised residual norm: "
3504 << resid <<
"\nafter " <<
Max_iter <<
" iterations.\n"
3509 std::string error_message_string =
"Solver failed to converge and you ";
3510 error_message_string +=
"requested an error on convergence failures.";
3512 error_message_string, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
void apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update:
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
DoubleVector * C_pt
Pointer to the row vector in the bordered system.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
CRDoubleMatrix * Matrix_pt
Pointer to matrix.
bool Iteration_restart
boolean indicating if iteration restarting is used
unsigned Restart
The number of iterations before the iteration proceedure is restarted if iteration restart is used.
void apply_schur_complement_preconditioner(const DoubleVector &rhs, DoubleVector &soln)
Apply the block-diagonal Schur complement preconditioner to compute the LHS which has size N+1 (the a...
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void generate_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Generate a plane rotation. This is done by finding the values of (i....
double Schur_complement_scalar
The scalar component of the Schur complement preconditioner.
void augmented_matrix_multiply(CRDoubleMatrix *matrix_pt, const DoubleVector &x, DoubleVector &soln)
Multiply the vector x by the augmented system matrix.
unsigned Iterations
Number of iterations taken.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
DoubleVector * B_pt
Pointer to the column vector in the bordered system.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
void update(const unsigned &k, const Vector< Vector< double > > &H, const Vector< double > &s, const Vector< DoubleVector > &v, DoubleVector &x)
Helper function to update the result vector using the result, x=x_0+V_m*y.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
A class for compressed row matrices. This is a distributable object.
int * row_start()
Access to C-style row_start array.
int * column_index()
Access to C-style column index array.
double * value()
Access to C-style value array.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
This is where the actual work is done – different implementations for different matrix types.
void solve(Problem *const &problem_pt, DoubleVector &result)
Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the r...
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
bool distributed() const
distribution is serial or distributed
unsigned nrow() const
access function to the number of global rows.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
virtual void multiply(const DoubleVector &x, DoubleVector &soln) const =0
Multiply the matrix by the vector x: soln=Ax.
virtual void residual(const DoubleVector &x, const DoubleVector &b, DoubleVector &residual_)
Find the residual, i.e. r=b-Ax the residual.
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
A vector in the mathematical sense, initially developed for linear algebra type applications....
void initialise(const double &v)
initialise the whole vector with value v
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
double norm() const
compute the 2 norm of this vector
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
double * values_pt()
access function to the underlying values
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
unsigned Iterations
Number of iterations taken.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
MATRIX * Matrix_pt
System matrix pointer in the format specified by the template argument.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
double Tolerance
Convergence tolerance.
bool Throw_error_after_max_iter
Should we throw an error instead of just returning when we hit the max iterations?
double Preconditioner_setup_time
Preconditioner setup time.
double Jacobian_setup_time
Jacobian setup time.
static IdentityPreconditioner Default_preconditioner
Default preconditioner: The base class for preconditioners is a fully functional (if trivial!...
unsigned Max_iter
Maximum number of iterations.
double Solution_time
linear solver solution time
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
std::ofstream Output_file_stream
Output file stream for convergence history.
bool Doc_convergence_history
Flag indicating if the convergence history is to be documented.
bool Setup_preconditioner_before_solve
indicates whether the preconditioner should be setup before solve. Default = true;
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
bool Doc_time
Boolean flag that indicates whether the time taken.
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
unsigned long ndof() const
Return the number of dofs.
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
void check_validity_of_solve_helper_inputs(MATRIX *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution, const double &n_dof)
Self-test to check that all the dimensions of the inputs to solve helper are consistent and everythin...
bool Use_as_smoother
When a derived class object is being used as a smoother in the MG solver (or elsewhere) the residual ...
unsigned Max_iter
Max. # of Newton iterations.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void clean_up_memory()
Clean up function that deletes anything dynamically allocated in this namespace.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...