17 #ifndef __deal2__solver_cg_h 18 #define __deal2__solver_cg_h 21 #include <deal.II/base/config.h> 22 #include <deal.II/lac/tridiagonal_matrix.h> 23 #include <deal.II/lac/solver.h> 24 #include <deal.II/lac/solver_control.h> 26 #include <deal.II/base/logstream.h> 27 #include <deal.II/base/subscriptor.h> 79 template <
class VECTOR = Vector<
double> >
133 const bool compute_condition_number =
false,
134 const bool compute_all_condition_numbers =
false,
135 const bool compute_eigenvalues =
false);
162 template <
class MATRIX,
class PRECONDITIONER>
167 const PRECONDITIONER &precondition);
190 const VECTOR &d)
const;
229 template <
class VECTOR>
237 log_coefficients (log_coefficients),
238 compute_condition_number(compute_condition_number),
239 compute_all_condition_numbers(compute_all_condition_numbers),
240 compute_eigenvalues(compute_eigenvalues)
245 template <
class VECTOR>
256 template <
class VECTOR>
266 template <
class VECTOR>
272 template <
class VECTOR>
276 return std::sqrt(
res2);
281 template <
class VECTOR>
293 template <
class VECTOR>
298 const VECTOR &)
const 303 template <
class VECTOR>
304 template <
class MATRIX,
class PRECONDITIONER>
309 const PRECONDITIONER &precondition)
317 Vz = this->
memory.alloc();
318 Vp = this->
memory.alloc();
324 double eigen_beta_alpha = 0;
328 std::vector<double> diagonal;
329 std::vector<double> offdiagonal;
346 double res,gh,alpha,beta;
369 precondition.vmult(h,g);
403 precondition.vmult(h,g);
420 deallog <<
"alpha-beta:" << alpha <<
'\t' << beta << std::endl;
427 diagonal.push_back(1./alpha + eigen_beta_alpha);
428 eigen_beta_alpha = beta/alpha;
429 offdiagonal.push_back(std::sqrt(beta)/alpha);
435 for (size_type i=0; i<diagonal.size(); ++i)
437 T(i,i) = diagonal[i];
438 if (i< diagonal.size()-1)
439 T(i,i+1) = offdiagonal[i];
441 T.compute_eigenvalues();
442 deallog <<
"Condition number estimate: " <<
443 T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
457 for (size_type i=0; i<diagonal.size(); ++i)
459 T(i,i) = diagonal[i];
460 if (i< diagonal.size()-1)
461 T(i,i+1) = offdiagonal[i];
463 T.compute_eigenvalues();
466 && (diagonal.size() > 1))
467 deallog <<
"Condition number estimate: " <<
468 T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
471 for (size_type i=0; i<T.n(); ++i)
472 deallog <<
' ' << T.eigenvalue(i);
473 deallog << std::endl;
482 this->
control().last_value()));
488 DEAL_II_NAMESPACE_CLOSE
VectorMemory< VECTOR > & memory
bool compute_condition_number
virtual State check(const unsigned int step, const double check_value)
void vmult(VECTOR &u, const VECTOR &v) const
#define AssertThrow(cond, exc)
AdditionalData(const bool log_coefficients=false, const bool compute_condition_number=false, const bool compute_all_condition_numbers=false, const bool compute_eigenvalues=false)
virtual double criterion()
SolverControl & control() const
unsigned int global_dof_index
Stop iteration, goal reached.
#define Assert(cond, exc)
bool compute_all_condition_numbers
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
void push(const std::string &text)
types::global_dof_index size_type
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
AdditionalData additional_data
SolverCG(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
::ExceptionBase & ExcDivideByZero()