Reference documentation for deal.II version 8.1.0
solver_cg.h
1 // ---------------------------------------------------------------------
2 // @f$Id: solver_cg.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 1998 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__solver_cg_h
18 #define __deal2__solver_cg_h
19 
20 
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>
28 #include <cmath>
29 
31 
32 // forward declaration
34 
35 
38 
79 template <class VECTOR = Vector<double> >
80 class SolverCG : public Solver<VECTOR>
81 {
82 public:
87 
93  {
100 
109 
118 
126 
132  AdditionalData (const bool log_coefficients = false,
133  const bool compute_condition_number = false,
134  const bool compute_all_condition_numbers = false,
135  const bool compute_eigenvalues = false);
136  };
137 
141  SolverCG (SolverControl &cn,
143  const AdditionalData &data = AdditionalData());
144 
150  SolverCG (SolverControl &cn,
151  const AdditionalData &data=AdditionalData());
152 
156  virtual ~SolverCG ();
157 
162  template <class MATRIX, class PRECONDITIONER>
163  void
164  solve (const MATRIX &A,
165  VECTOR &x,
166  const VECTOR &b,
167  const PRECONDITIONER &precondition);
168 
169 protected:
176  virtual double criterion();
177 
187  virtual void print_vectors(const unsigned int step,
188  const VECTOR &x,
189  const VECTOR &r,
190  const VECTOR &d) const;
191 
198  VECTOR *Vr;
199  VECTOR *Vp;
200  VECTOR *Vz;
201 
212  double res2;
213 
218 
219 private:
220  void cleanup();
221 };
222 
225 /*------------------------- Implementation ----------------------------*/
226 
227 #ifndef DOXYGEN
228 
229 template <class VECTOR>
230 inline
232 AdditionalData (const bool log_coefficients,
233  const bool compute_condition_number,
235  const bool compute_eigenvalues)
236  :
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)
241 {}
242 
243 
244 
245 template <class VECTOR>
248  const AdditionalData &data)
249  :
250  Solver<VECTOR>(cn,mem),
251  additional_data(data)
252 {}
253 
254 
255 
256 template <class VECTOR>
258  const AdditionalData &data)
259  :
260  Solver<VECTOR>(cn),
261  additional_data(data)
262 {}
263 
264 
265 
266 template <class VECTOR>
268 {}
269 
270 
271 
272 template <class VECTOR>
273 double
275 {
276  return std::sqrt(res2);
277 }
278 
279 
280 
281 template <class VECTOR>
282 void
284 {
285  this->memory.free(Vr);
286  this->memory.free(Vp);
287  this->memory.free(Vz);
288  deallog.pop();
289 }
290 
291 
292 
293 template <class VECTOR>
294 void
295 SolverCG<VECTOR>::print_vectors(const unsigned int,
296  const VECTOR &,
297  const VECTOR &,
298  const VECTOR &) const
299 {}
300 
301 
302 
303 template <class VECTOR>
304 template <class MATRIX, class PRECONDITIONER>
305 void
307  VECTOR &x,
308  const VECTOR &b,
309  const PRECONDITIONER &precondition)
310 {
312 
313  deallog.push("cg");
314 
315  // Memory allocation
316  Vr = this->memory.alloc();
317  Vz = this->memory.alloc();
318  Vp = this->memory.alloc();
319  // Should we build the matrix for
320  // eigenvalue computations?
321  bool do_eigenvalues = additional_data.compute_condition_number
324  double eigen_beta_alpha = 0;
325 
326  // vectors used for eigenvalue
327  // computations
328  std::vector<double> diagonal;
329  std::vector<double> offdiagonal;
330 
331  try
332  {
333  // define some aliases for simpler access
334  VECTOR &g = *Vr;
335  VECTOR &d = *Vz;
336  VECTOR &h = *Vp;
337  // resize the vectors, but do not set
338  // the values since they'd be overwritten
339  // soon anyway.
340  g.reinit(x, true);
341  d.reinit(x, true);
342  h.reinit(x, true);
343  // Implementation taken from the DEAL
344  // library
345  int it=0;
346  double res,gh,alpha,beta;
347 
348  // compute residual. if vector is
349  // zero, then short-circuit the
350  // full computation
351  if (!x.all_zero())
352  {
353  A.vmult(g,x);
354  g.add(-1.,b);
355  }
356  else
357  g.equ(-1.,b);
358  res = g.l2_norm();
359 
360  conv = this->control().check(0,res);
361  if (conv)
362  {
363  cleanup();
364  return;
365  }
366 
368  {
369  precondition.vmult(h,g);
370 
371  d.equ(-1.,h);
372 
373  gh = g*h;
374  }
375  else
376  {
377  d.equ(-1.,g);
378  gh = res*res;
379  }
380 
381  while (conv == SolverControl::iterate)
382  {
383  it++;
384  A.vmult(h,d);
385 
386  alpha = d*h;
387  Assert(alpha != 0., ExcDivideByZero());
388  alpha = gh/alpha;
389 
390  g.add(alpha,h);
391  x.add(alpha,d);
392  res = g.l2_norm();
393 
394  print_vectors(it, x, g, d);
395 
396  conv = this->control().check(it,res);
397  if (conv != SolverControl::iterate)
398  break;
399 
401  == false)
402  {
403  precondition.vmult(h,g);
404 
405  beta = gh;
406  Assert(beta != 0., ExcDivideByZero());
407  gh = g*h;
408  beta = gh/beta;
409  d.sadd(beta,-1.,h);
410  }
411  else
412  {
413  beta = gh;
414  gh = res*res;
415  beta = gh/beta;
416  d.sadd(beta,-1.,g);
417  }
418 
420  deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl;
421  // set up the vectors
422  // containing the diagonal
423  // and the off diagonal of
424  // the projected matrix.
425  if (do_eigenvalues)
426  {
427  diagonal.push_back(1./alpha + eigen_beta_alpha);
428  eigen_beta_alpha = beta/alpha;
429  offdiagonal.push_back(std::sqrt(beta)/alpha);
430  }
431 
432  if (additional_data.compute_all_condition_numbers && (diagonal.size()>1))
433  {
434  TridiagonalMatrix<double> T(diagonal.size(), true);
435  for (size_type i=0; i<diagonal.size(); ++i)
436  {
437  T(i,i) = diagonal[i];
438  if (i< diagonal.size()-1)
439  T(i,i+1) = offdiagonal[i];
440  }
441  T.compute_eigenvalues();
442  deallog << "Condition number estimate: " <<
443  T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
444  }
445  }
446  }
447  catch (...)
448  {
449  cleanup();
450  throw;
451  }
452 
453  // Write eigenvalues or condition number
454  if (do_eigenvalues)
455  {
456  TridiagonalMatrix<double> T(diagonal.size(), true);
457  for (size_type i=0; i<diagonal.size(); ++i)
458  {
459  T(i,i) = diagonal[i];
460  if (i< diagonal.size()-1)
461  T(i,i+1) = offdiagonal[i];
462  }
463  T.compute_eigenvalues();
466  && (diagonal.size() > 1))
467  deallog << "Condition number estimate: " <<
468  T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
470  {
471  for (size_type i=0; i<T.n(); ++i)
472  deallog << ' ' << T.eigenvalue(i);
473  deallog << std::endl;
474  }
475  }
476 
477  // Deallocate Memory
478  cleanup();
479  // in case of failure: throw exception
480  if (this->control().last_check() != SolverControl::success)
481  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
482  this->control().last_value()));
483  // otherwise exit as normal
484 }
485 
486 #endif // DOXYGEN
487 
488 DEAL_II_NAMESPACE_CLOSE
489 
490 #endif
VectorMemory< VECTOR > & memory
Definition: solver.h:210
void pop()
virtual State check(const unsigned int step, const double check_value)
void vmult(VECTOR &u, const VECTOR &v) const
Continue iteration.
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
VECTOR * Vr
Definition: solver_cg.h:198
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
Definition: solver.h:239
unsigned int global_dof_index
Definition: types.h:100
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:299
virtual void print_vectors(const unsigned int step, const VECTOR &x, const VECTOR &r, const VECTOR &d) const
Definition: solver.h:147
double res2
Definition: solver_cg.h:212
void push(const std::string &text)
types::global_dof_index size_type
Definition: solver_cg.h:86
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
AdditionalData additional_data
Definition: solver_cg.h:217
virtual ~SolverCG()
SolverCG(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
::ExceptionBase & ExcDivideByZero()