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>
233 const bool compute_condition_number,
234 const bool compute_all_condition_numbers,
235 const bool compute_eigenvalues)
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>
248 const AdditionalData &data)
251 additional_data(data)
256 template <
class VECTOR>
258 const AdditionalData &data)
261 additional_data(data)
266 template <
class VECTOR>
272 template <
class VECTOR>
276 return std::sqrt(res2);
281 template <
class VECTOR>
285 this->memory.free(Vr);
286 this->memory.free(Vp);
287 this->memory.free(Vz);
293 template <
class VECTOR>
298 const VECTOR &)
const
303 template <
class VECTOR>
304 template <
class MATRIX,
class PRECONDITIONER>
309 const PRECONDITIONER &precondition)
316 Vr = this->memory.alloc();
317 Vz = this->memory.alloc();
318 Vp = this->memory.alloc();
321 bool do_eigenvalues = additional_data.compute_condition_number
322 | additional_data.compute_all_condition_numbers
323 | additional_data.compute_eigenvalues;
324 double eigen_beta_alpha = 0;
328 std::vector<double> diagonal;
329 std::vector<double> offdiagonal;
346 double res,gh,alpha,beta;
360 conv = this->control().check(0,res);
369 precondition.vmult(h,g);
394 print_vectors(it, x, g, d);
396 conv = this->control().check(it,res);
403 precondition.vmult(h,g);
419 if (additional_data.log_coefficients)
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);
432 if (additional_data.compute_all_condition_numbers && (diagonal.size()>1))
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();
464 if (additional_data.compute_condition_number
465 && ! additional_data.compute_all_condition_numbers
466 && (diagonal.size() > 1))
467 deallog <<
"Condition number estimate: " <<
468 T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
469 if (additional_data.compute_eigenvalues)
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
bool compute_condition_number
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()
unsigned int global_dof_index
#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()
Stop iteration, goal reached.