27#ifndef OOMPH_COMPLEX_SMOOTHER_HEADER
28#define OOMPH_COMPLEX_SMOOTHER_HEADER
32#include <oomph-lib-config.h>
97 <<
"Can only deal with two input vectors. You have " << x.
size()
98 <<
" vectors." << std::endl;
113 <<
"Can only deal with two output vectors. You have " <<
soln.
size()
114 <<
" output vectors." << std::endl;
149 matrices_pt[1]->multiply(x[1],
temp);
163 template<
typename MATRIX>
169 const double&
n_dof);
184 template<
typename MATRIX>
260 <<
"must be the same as the rhs distribution";
280 template<
typename MATRIX>
342 double omega = (12.0 - 4.0 *
pow(
kh, 2.0)) / (18.0 - 3.0 *
pow(
kh, 2.0));
360 else if ((
k >
pi) && (
kh > 2 *
cos(
pi * h / 2)))
454 for (
unsigned i = 0;
i <
n_dof;
i++)
531 template<
typename MATRIX>
536 unsigned n_dof = Matrix_real_pt->nrow();
559 if ((!
solution[0].distribution_pt()->built()) ||
560 (!
solution[1].distribution_pt()->built()))
612 for (
unsigned i = 0;
i <
n_dof;
i++)
665 if (!Use_as_smoother)
694 if (Doc_convergence_history)
696 if (!Output_file_stream.is_open())
702 Output_file_stream << Iterations <<
" " <<
norm_res << std::endl;
725 for (
unsigned i = 0;
i <
n_dof;
i++)
740 dummy_r *= Omega * Matrix_diagonal_inv_sq[
i];
741 dummy_c *= Omega * Matrix_diagonal_inv_sq[
i];
769 if (!Use_as_smoother)
790 if (Doc_convergence_history)
792 if (!Output_file_stream.is_open())
798 Output_file_stream << Iterations <<
" " <<
norm_res << std::endl;
812 if (!Use_as_smoother)
817 oomph_info <<
"\n(Complex) damped Jacobi converged. Residual norm: "
819 <<
"\nNumber of iterations to convergence: " << Iterations
834 oomph_info <<
"Time for solve with (complex) damped Jacobi [sec]: "
835 << Solution_time << std::endl;
840 if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
842 std::string error_message =
843 "Solver failed to converge and you requested ";
844 error_message +=
"an error on convergence failures.";
853 template<
typename MATRIX>
895 std::string error_message =
"Solve function for class\n\n";
896 error_message +=
"ComplexGMRES\n\n";
897 error_message +=
"is deliberately broken to avoid the accidental \n";
898 error_message +=
"use of the inappropriate C++ default. If you \n";
899 error_message +=
"really need one for this class, write it yourself...\n";
1031 const Vector<std::complex<double>>&
s,
1053 for (
int i =
int(
k);
i >= 0;
i--)
1059 for (
int j =
i - 1;
j >= 0;
j--)
1068 unsigned n_row = x[0].nrow();
1073 for (
unsigned j = 0;
j <=
k;
j++)
1076 const double*
vj_r_pt =
v[
j][0].values_pt();
1079 const double*
vj_c_pt =
v[
j][1].values_pt();
1082 for (
unsigned i = 0;
i <
n_row;
i++)
1107 std::complex<double>&
dy,
1108 std::complex<double>&
cs,
1109 std::complex<double>&
sn)
1125 else if (std::abs(
dy) > std::abs(
dx))
1128 std::complex<double>
temp =
dx /
dy;
1146 std::complex<double>
temp =
dy /
dx;
1171 <<
"and/or nonnegative. Value is: " <<
sn
1186 std::complex<double>&
dy,
1187 std::complex<double>&
cs,
1188 std::complex<double>&
sn)
1191 std::complex<double>
temp = std::conj(
cs) *
dx + std::conj(
sn) *
dy;
1218 template<
typename MATRIX>
1228 unsigned n_row = Matrices_storage_pt[0]->nrow();
1232 if (Max_iter >
n_row)
1239 <<
"the number of rows in the problem."
1240 <<
"\nMaximum number of iterations: " << Max_iter
1241 <<
"\nNumber of rows: " <<
n_row << std::endl;
1256 Matrices_storage_pt[
dof_type]) != 0)
1283 <<
" linear system";
1289 if (
rhs[
dof_type].distribution_pt()->distributed())
1306 <<
"must be the same as the rhs distribution";
1355 w[
dof_type].build(this->distribution_pt(), 0.0);
1366 r[
dof_type].build(this->distribution_pt(), 0.0);
1386 double beta =
normb;
1405 if (Doc_convergence_history)
1408 if (!Output_file_stream.is_open())
1417 Output_file_stream << 0 <<
" " <<
resid << std::endl;
1422 if (
resid <= Tolerance)
1428 oomph_info <<
"GMRES converged immediately. Normalised residual norm: "
1429 <<
resid << std::endl;
1462 v[0][
dof_type].build(this->distribution_pt(), 0.0);
1475 for (
unsigned i = 0;
i <
n_row;
i++)
1488 for (
unsigned j = 0;
j < Max_iter;
j++)
1491 hessenberg[
j].resize(
j + 2, std::complex<double>(0.0, 0.0));
1495 complex_matrix_multiplication(Matrices_storage_pt,
v[
j], w);
1498 double*
w_r_pt = w[0].values_pt();
1501 double*
w_c_pt = w[1].values_pt();
1507 for (
unsigned i = 0;
i <
j + 1;
i++)
1510 const double*
vi_r_pt =
v[
i][0].values_pt();
1513 const double*
vi_c_pt =
v[
i][1].values_pt();
1516 for (
unsigned k = 0;
k <
n_row;
k++)
1531 for (
unsigned k = 0;
k <
n_row;
k++)
1554 v[
j + 1][0].build(this->distribution_pt(), 0.0);
1557 v[
j + 1][1].build(this->distribution_pt(), 0.0);
1564 double*
v_r_pt =
v[
j + 1][0].values_pt();
1567 double*
v_c_pt =
v[
j + 1][1].values_pt();
1570 const double*
w_r_pt = w[0].values_pt();
1573 const double*
w_c_pt = w[1].values_pt();
1583 for (
unsigned k = 0;
k <
n_row;
k++)
1601 oomph_info <<
"Subdiagonal Hessenberg entry is zero."
1602 <<
"Do something here..." << std::endl;
1607 for (
unsigned k = 0;
k <
j;
k++)
1611 apply_plane_rotation(
1616 generate_plane_rotation(
1620 apply_plane_rotation(
1625 apply_plane_rotation(
s[
j],
s[
j + 1],
cs[
j],
sn[
j]);
1630 beta = std::abs(
s[
j + 1]);
1637 if (Doc_convergence_history)
1640 if (!Output_file_stream.is_open())
1649 Output_file_stream <<
j + 1 <<
" " <<
resid << std::endl;
1654 if (
resid < Tolerance)
1667 oomph_info <<
"\nGMRES converged (1). Normalised residual norm: "
1668 <<
resid << std::endl;
1671 oomph_info <<
"Number of iterations to convergence: " <<
j + 1 <<
"\n"
1685 oomph_info <<
"Time for solve with GMRES [sec]: " << Solution_time
1696 Iterations = Max_iter;
1716 oomph_info <<
"Time for solve with GMRES [sec]: " << Solution_time
1733 template<
typename MATRIX>
1777 <<
"HelmholtzGMRESMG is only meant to be used as "
1778 <<
"a linear solver.\n";
1795 <<
"only meant to be used as a linear solver.\n";
1814 OomphLibWarning(
"Can't guarantee the MG solver will work in parallel!",
1887 if (nblock_types != 2)
1890 std::stringstream
tmp;
1891 tmp <<
"There are supposed to be two block types.\nYours has "
1922 oomph_info <<
"\nTime for setup of block Jacobian [sec]: "
1930 if ((!(*
result.distribution_pt() ==
1972 <<
"This function is broken. The block preconditioner "
1973 <<
"needs access to the underlying mesh.\n";
2004 throw OomphLibError(
"No matrix was stored -- cannot re-solve",
2115 <<
"Can only deal with two input vectors. You have " << x.
size()
2116 <<
" vectors." << std::endl;
2131 <<
"Can only deal with two output vectors. You have " <<
soln.
size()
2132 <<
" output vectors." << std::endl;
2186 const Vector<std::complex<double>>&
s,
2209 for (
int i =
int(
k);
i >= 0;
i--)
2215 for (
int j =
i - 1;
j >= 0;
j--)
2224 unsigned n_row = x[0].nrow();
2249 for (
unsigned j = 0;
j <=
k;
j++)
2252 const double*
vj_r_pt =
v[
j][0].values_pt();
2255 const double*
vj_c_pt =
v[
j][1].values_pt();
2258 for (
unsigned i = 0;
i <
n_row;
i++)
2328 std::complex<double>&
dy,
2329 std::complex<double>&
cs,
2330 std::complex<double>&
sn)
2346 else if (std::abs(
dy) > std::abs(
dx))
2349 std::complex<double>
temp =
dx /
dy;
2367 std::complex<double>
temp =
dy /
dx;
2392 <<
"and/or nonnegative. Value is: " <<
sn
2407 std::complex<double>&
dy,
2408 std::complex<double>&
cs,
2409 std::complex<double>&
sn)
2412 std::complex<double>
temp = std::conj(
cs) *
dx + std::conj(
sn) *
dy;
2452 template<
typename MATRIX>
2463 if (n_dof_types != 2)
2482 unsigned n_row = Matrices_storage_pt[0]->nrow();
2486 if (Max_iter >
n_row)
2493 <<
"the number of rows in the problem."
2494 <<
"\nMaximum number of iterations: " << Max_iter
2495 <<
"\nNumber of rows: " <<
n_row << std::endl;
2510 Matrices_storage_pt[
dof_type]) != 0)
2538 <<
" linear system";
2544 if (
rhs.distribution_pt()->distributed())
2556 if (!(*
rhs.distribution_pt() == *
solution.distribution_pt()))
2560 <<
"must be the same as the rhs distribution";
2641 if (Setup_preconditioner_before_solve)
2647 preconditioner_pt()->setup(
dynamic_cast<MATRIX*
>(matrix_pt));
2657 oomph_info <<
"Time for setup of preconditioner [sec]: "
2658 << Preconditioner_setup_time << std::endl;
2668 oomph_info <<
"Setup of preconditioner is bypassed in resolve mode"
2686 if (Preconditioner_LHS)
2692 this->return_block_vectors(
block_r,
r);
2695 preconditioner_pt()->preconditioner_solve(
rhs,
r);
2698 this->get_block_vectors(
r,
block_r);
2720 double beta =
normb;
2739 if (Doc_convergence_history)
2742 if (!Output_file_stream.is_open())
2751 Output_file_stream << 0 <<
" " <<
resid << std::endl;
2756 if (
resid <= Tolerance)
2762 oomph_info <<
"GMRES converged immediately. Normalised residual norm: "
2763 <<
resid << std::endl;
2809 for (
unsigned i = 0;
i <
n_row;
i++)
2822 for (
unsigned j = 0;
j < Max_iter;
j++)
2825 hessenberg[
j].resize(
j + 2, std::complex<double>(0.0, 0.0));
2849 if (Preconditioner_LHS)
2853 complex_matrix_multiplication(
2860 this->return_block_vectors(
block_w, w);
2863 this->preconditioner_pt()->preconditioner_solve(
temp, w);
2866 this->get_block_vectors(w,
block_w);
2875 this->preconditioner_pt()->preconditioner_solve(
vj,
temp);
2882 complex_matrix_multiplication(
2897 for (
unsigned i = 0;
i <
j + 1;
i++)
2906 for (
unsigned k = 0;
k <
n_row;
k++)
2921 for (
unsigned k = 0;
k <
n_row;
k++)
2976 for (
unsigned k = 0;
k <
n_row;
k++)
2994 oomph_info <<
"Subdiagonal Hessenberg entry is zero. "
2995 <<
"Do something here..." << std::endl;
3000 for (
unsigned k = 0;
k <
j;
k++)
3004 apply_plane_rotation(
3009 generate_plane_rotation(
3013 apply_plane_rotation(
3018 apply_plane_rotation(
s[
j],
s[
j + 1],
cs[
j],
sn[
j]);
3023 beta = std::abs(
s[
j + 1]);
3030 if (Doc_convergence_history)
3033 if (!Output_file_stream.is_open())
3042 Output_file_stream <<
j + 1 <<
" " <<
resid << std::endl;
3047 if (
resid < Tolerance)
3063 oomph_info <<
"\nGMRES converged (1). Normalised residual norm: "
3064 <<
resid << std::endl;
3067 oomph_info <<
"Number of iterations to convergence: " <<
j + 1 <<
"\n"
3081 oomph_info <<
"Time for solve with GMRES [sec]: " << Solution_time
3092 Iterations = Max_iter;
3118 complex_matrix_multiplication(
3134 for (
unsigned i = 0;
i <
n_row;
i++)
3144 if (Preconditioner_LHS)
3156 this->return_block_vectors(
block_r,
r);
3159 preconditioner_pt()->preconditioner_solve(
temp,
r);
3179 if (
resid < Tolerance)
3185 oomph_info <<
"\nGMRES converged (2). Normalised residual norm: "
3187 <<
"\nNumber of iterations to convergence: " << Iterations
3201 oomph_info <<
"Time for solve with GMRES [sec]: " << Solution_time
3208 oomph_info <<
"\nGMRES did not converge to required tolerance! "
3209 <<
"\nReturning with normalised residual norm: " <<
resid
3210 <<
"\nafter " << Max_iter <<
" iterations.\n"
3214 if (Throw_error_after_max_iter)
3216 std::string
err =
"Solver failed to converge and you requested an error";
3217 err +=
" on convergence failures.";
3239 template<
typename MATRIX>
3272 <<
"only capable of using right preconditioning."
3292 OomphLibWarning(
"Can't guarantee the MG solver will work in parallel!",
3365 if (nblock_types != 2)
3368 std::stringstream
tmp;
3369 tmp <<
"There are supposed to be two block types.\nYours has "
3400 oomph_info <<
"\nTime for setup of block Jacobian [sec]: "
3408 if ((!(*
result.distribution_pt() ==
3448 const Vector<std::complex<double>>&
s,
3471 for (
int i =
int(
k);
i >= 0;
i--)
3477 for (
int j =
i - 1;
j >= 0;
j--)
3486 unsigned n_row = x[0].nrow();
3505 for (
unsigned j = 0;
j <=
k;
j++)
3514 for (
unsigned i = 0;
i <
n_row;
i++)
3550 template<
typename MATRIX>
3561 if (n_dof_types != 2)
3580 unsigned n_row = this->Matrices_storage_pt[0]->nrow();
3584 if (this->Max_iter >
n_row)
3591 <<
"the number of rows in the problem."
3592 <<
"\nMaximum number of iterations: "
3593 << this->Max_iter <<
"\nNumber of rows: " <<
n_row
3609 this->Matrices_storage_pt[
dof_type]) != 0)
3612 this->Matrices_storage_pt[
dof_type])
3637 <<
" linear system";
3643 if (
rhs.distribution_pt()->distributed())
3655 if (!(*
rhs.distribution_pt() == *
solution.distribution_pt()))
3659 <<
"must be the same as the rhs distribution";
3737 if (!(this->Resolving))
3740 if (this->Setup_preconditioner_before_solve)
3746 this->preconditioner_pt()->setup(
dynamic_cast<MATRIX*
>(matrix_pt));
3756 oomph_info <<
"Time for setup of preconditioner [sec]: "
3757 << this->Preconditioner_setup_time << std::endl;
3767 oomph_info <<
"Setup of preconditioner is bypassed in resolve mode"
3800 double beta =
normb;
3819 if (this->Doc_convergence_history)
3822 if (!(this->Output_file_stream.is_open()))
3831 this->Output_file_stream << 0 <<
" " <<
resid << std::endl;
3842 oomph_info <<
"FGMRES converged immediately. Normalised residual norm: "
3843 <<
resid << std::endl;
3864 for (
unsigned i = 0;
i <
n_row + 1;
i++)
3898 for (
unsigned i = 0;
i <
n_row;
i++)
3911 for (
unsigned j = 0;
j < this->Max_iter;
j++)
3914 hessenberg[
j].resize(
j + 2, std::complex<double>(0.0, 0.0));
3939 this->preconditioner_pt()->preconditioner_solve(
vj,
zj);
3946 this->complex_matrix_multiplication(
3960 for (
unsigned i = 0;
i <
j + 1;
i++)
3969 for (
unsigned k = 0;
k <
n_row;
k++)
3984 for (
unsigned k = 0;
k <
n_row;
k++)
4039 for (
unsigned k = 0;
k <
n_row;
k++)
4057 oomph_info <<
"Subdiagonal Hessenberg entry is zero. "
4058 <<
"Do something here..." << std::endl;
4063 for (
unsigned k = 0;
k <
j;
k++)
4067 this->apply_plane_rotation(
4072 this->generate_plane_rotation(
4076 this->apply_plane_rotation(
4081 this->apply_plane_rotation(
s[
j],
s[
j + 1],
cs[
j],
sn[
j]);
4086 beta = std::abs(
s[
j + 1]);
4093 if (this->Doc_convergence_history)
4096 if (!(this->Output_file_stream.is_open()))
4105 this->Output_file_stream <<
j + 1 <<
" " <<
resid << std::endl;
4113 this->Iterations =
j + 1;
4126 oomph_info <<
"\nFGMRES converged (1). Normalised residual norm: "
4127 <<
resid << std::endl;
4130 oomph_info <<
"Number of iterations to convergence: " <<
j + 1 <<
"\n"
4144 oomph_info <<
"Time for solve with FGMRES [sec]: "
4145 << this->Solution_time << std::endl;
4155 this->Iterations = this->Max_iter;
4158 if (this->Max_iter > 0)
4190 oomph_info <<
"\nFGMRES converged (2). Normalised residual norm: "
4191 <<
resid <<
"\nNumber of iterations to convergence: "
4192 << this->Iterations <<
"\n"
4205 oomph_info <<
"Time for solve with FGMRES [sec]: "
4206 << this->Solution_time << std::endl;
4212 oomph_info <<
"\nFGMRES did not converge to required tolerance! "
4213 <<
"\nReturning with normalised residual norm: " <<
resid
4214 <<
"\nafter " << this->Max_iter <<
" iterations.\n"
4218 if (this->Throw_error_after_max_iter)
4220 std::string
err =
"Solver failed to converge and you requested an error";
4221 err +=
" on convergence failures.";
Block Preconditioner base class. The block structure of the overall problem is determined from the Me...
void return_block_vectors(const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
Takes the vector of block vectors, s, and copies its entries into the naturally ordered vector,...
void get_block(const unsigned &i, const unsigned &j, MATRIX &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
unsigned nblock_types() const
Return the number of block types.
void get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Takes the naturally ordered vector and rearranges it into a vector of sub vectors corresponding to th...
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
MATRIX * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*,...
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Set the i-th mesh for this block preconditioner. Note: The method set_nmesh(...) must be called befor...
A class for compressed row matrices. This is a distributable object.
unsigned long nrow() const
Return the number of rows of the matrix.
Vector< double > diagonal_entries() const
returns a Vector of diagonal entries of this matrix. This only works with square matrices....
Damped Jacobi "solver" templated by matrix type. The "solver" exists in many different incarnations: ...
void complex_solve_helper(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
This is where the actual work is done.
~ComplexDampedJacobi()
Empty destructor.
void operator=(const ComplexDampedJacobi &)=delete
Broken assignment operator.
unsigned iterations() const
Number of iterations taken.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrices pointed to by Matrix_real_pt and Matrix_imag_pt can be delet...
void complex_smoother_solve(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
void complex_smoother_setup(Vector< CRDoubleMatrix * > helmholtz_matrix_pt)
Setup: Pass pointer to the matrix and store in cast form.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
Vector< double > Matrix_diagonal_inv_sq
Vector whose j'th entry is given by 1/(A_r(j,j)^2 + A_c(j,j)^2)
ComplexDampedJacobi(const double &omega=0.5)
Constructor (empty)
Vector< double > Matrix_diagonal_real
Vector containing the diagonal entries of A_r (real(A))
Vector< double > Matrix_diagonal_imag
Vector containing the diagonal entries of A_c (imag(A))
CRDoubleMatrix * Matrix_imag_pt
Pointer to the real part of the system matrix.
double Omega
Damping factor.
double & omega()
Get access to the value of Omega (lvalue)
CRDoubleMatrix * Matrix_real_pt
Pointer to the real part of the system matrix.
void solve(Problem *const &problem_pt, DoubleVector &result)
Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the r...
ComplexDampedJacobi(const ComplexDampedJacobi &)=delete
Broken copy constructor.
unsigned Iterations
Number of iterations taken.
void calculate_omega(const double &k, const double &h)
Function to calculate the value of Omega by passing in the value of k and h [see Elman et al....
The GMRES method rewritten for complex matrices.
void generate_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Generate a plane rotation. This is done by finding the value of (i....
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void update(const unsigned &k, const Vector< Vector< std::complex< double > > > &hessenberg, const Vector< std::complex< double > > &s, const Vector< Vector< DoubleVector > > &v, Vector< DoubleVector > &x)
Helper function to update the result vector.
unsigned iterations() const
Number of iterations taken.
void complex_smoother_setup(Vector< CRDoubleMatrix * > helmholtz_matrix_pt)
Setup: Pass pointer to the matrix and store in cast form.
unsigned Iterations
Number of iterations taken.
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
ComplexGMRES(const ComplexGMRES &)=delete
Broken copy constructor.
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 complex_solve_helper(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
This is where the actual work is done.
void complex_smoother_solve(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
void apply_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Apply plane rotation. This is done using the update:
ComplexGMRES()
Constructor.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
bool Matrix_can_be_deleted
Boolean flag to indicate if the real and imaginary system matrices can be deleted.
void operator=(const ComplexGMRES &)=delete
Broken assignment operator.
~ComplexGMRES()
Empty destructor.
Vector< CRDoubleMatrix * > Matrices_storage_pt
Vector of pointers to the real and imaginary part of the system matrix.
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,...
A vector in the mathematical sense, initially developed for linear algebra type applications....
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
virtual unsigned ndof_types() const
The number of types of degrees of freedom in this element are sub-divided into.
The FGMRES method, i.e. the flexible variant of the GMRES method which allows for nonconstant precond...
HelmholtzFGMRESMG(const HelmholtzFGMRESMG &)=delete
Broken copy constructor.
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...
virtual ~HelmholtzFGMRESMG()
Destructor (cleanup storage)
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
HelmholtzFGMRESMG()
Constructor (empty)
void update(const unsigned &k, const Vector< Vector< std::complex< double > > > &hessenberg, const Vector< std::complex< double > > &s, const Vector< Vector< DoubleVector > > &z_m, Vector< DoubleVector > &x)
Helper function to update the result vector.
void operator=(const HelmholtzFGMRESMG &)=delete
Broken assignment operator.
void set_preconditioner_LHS()
Overloaded function to let the user know that left preconditioning is not possible with FGMRES,...
The GMRES method for the Helmholtz solver.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
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 solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
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...
Vector< CRDoubleMatrix * > Matrices_storage_pt
Vector of pointers to the real and imaginary part of the system matrix.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Implementation of the pure virtual base class function. The function has been broken because this is ...
void apply_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Apply plane rotation. This is done using the update:
HelmholtzGMRESMG(const HelmholtzGMRESMG &)=delete
Broken copy constructor.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void set_preconditioner_LHS()
Set left preconditioning (the default)
virtual ~HelmholtzGMRESMG()
Destructor (cleanup storage)
void complex_matrix_multiplication(Vector< CRDoubleMatrix * > const matrices_pt, const Vector< DoubleVector > &x, Vector< DoubleVector > &soln)
Helper function to calculate a complex matrix-vector product. Assumes the matrix has been provided as...
void set_preconditioner_RHS()
Enable right preconditioning.
HelmholtzGMRESMG()
Constructor.
void update(const unsigned &k, const Vector< Vector< std::complex< double > > > &hessenberg, const Vector< std::complex< double > > &s, const Vector< Vector< DoubleVector > > &v, Vector< DoubleVector > &x)
Helper function to update the result vector.
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
unsigned iterations() const
Number of iterations taken.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
unsigned Iterations
Number of iterations taken.
void setup()
Implementation of the pure virtual base class function. This accompanies the preconditioner_solve fun...
void operator=(const HelmholtzGMRESMG &)=delete
Broken assignment operator.
void generate_plane_rotation(std::complex< double > &dx, std::complex< double > &dy, std::complex< double > &cs, std::complex< double > &sn)
Helper function: Generate a plane rotation. This is done by finding the value of (i....
Helmholtz smoother class: The smoother class is designed for the Helmholtz equation to be used in con...
void check_validity_of_solve_helper_inputs(CRDoubleMatrix *const &real_matrix_pt, CRDoubleMatrix *const &imag_matrix_pt, const Vector< DoubleVector > &rhs, Vector< 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 algorithm the residual norm does no...
virtual void complex_smoother_setup(Vector< CRDoubleMatrix * > matrix_pt)=0
Setup the smoother for the matrix specified by the pointer.
HelmholtzSmoother()
Empty constructor.
virtual void complex_smoother_solve(const Vector< DoubleVector > &rhs, Vector< DoubleVector > &result)=0
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
virtual ~HelmholtzSmoother()
Virtual empty destructor.
void complex_matrix_multiplication(Vector< CRDoubleMatrix * > matrices_pt, const Vector< DoubleVector > &x, Vector< DoubleVector > &soln)
Helper function to calculate a complex matrix-vector product. Assumes the matrix has been provided as...
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
double & tolerance()
Access to convergence tolerance.
double Jacobian_setup_time
Jacobian setup time.
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
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...
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
static OomphCommunicator * communicator_pt()
access to global communicator. This is the oomph-lib equivalent of MPI_COMM_WORLD
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....
virtual void set_comm_pt(const OomphCommunicator *const comm_pt)
Set the communicator pointer.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
virtual void set_matrix_pt(DoubleMatrixBase *matrix_pt)
Set the matrix pointer.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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...
Mesh *& mesh_pt()
Return a pointer to the global mesh.
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
TAdvectionDiffusionReactionElement()
Constructor: Call constructors for TElement and AdvectionDiffusionReaction equations.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
const double Pi
50 digits from maple
double timer()
returns the time in seconds after some point in past
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...