Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mg_coarse.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_mg_coarse_h
17 #define dealii_mg_coarse_h
18 
19 
20 #include <deal.II/lac/full_matrix.h>
21 #include <deal.II/lac/householder.h>
22 #include <deal.II/lac/linear_operator.h>
23 #include <deal.II/lac/matrix_lib.h>
24 
25 #include <deal.II/multigrid/mg_base.h>
26 
27 DEAL_II_NAMESPACE_OPEN
28 
31 
38 template <class VectorType = Vector<double>>
39 class MGCoarseGridApplySmoother : public MGCoarseGridBase<VectorType>
40 {
41 public:
46 
51 
55  void
56  clear();
57 
61  void
63 
67  void
68  operator()(const unsigned int level,
69  VectorType & dst,
70  const VectorType & src) const override;
71 
72 private:
79 };
80 
81 
95 template <typename SolverType, class VectorType = Vector<double>>
96 class DEAL_II_DEPRECATED MGCoarseGridLACIteration
97  : public MGCoarseGridBase<VectorType>
98 {
99 public:
104 
109  template <typename MatrixType, typename PreconditionerType>
110  MGCoarseGridLACIteration(SolverType &,
111  const MatrixType &,
112  const PreconditionerType &);
113 
118 
122  template <typename MatrixType, typename PreconditionerType>
123  void
124  initialize(SolverType &, const MatrixType &, const PreconditionerType &);
125 
129  void
130  clear();
131 
136  void
137  operator()(const unsigned int level,
138  VectorType & dst,
139  const VectorType & src) const;
140 
145  template <typename MatrixType>
146  void
147  set_matrix(const MatrixType &);
148 
149 private:
155 
160 
165 };
166 
167 
168 
175 template <class VectorType,
176  class SolverType,
177  class MatrixType,
178  class PreconditionerType>
180 {
181 public:
186 
191  MGCoarseGridIterativeSolver(SolverType & solver,
192  const MatrixType & matrix,
193  const PreconditionerType &precondition);
194 
199  void
200  initialize(SolverType & solver,
201  const MatrixType & matrix,
202  const PreconditionerType &precondition);
203 
207  void
208  clear();
209 
214  virtual void
215  operator()(const unsigned int level,
216  VectorType & dst,
217  const VectorType & src) const override;
218 
219 private:
223  SmartPointer<SolverType,
224  MGCoarseGridIterativeSolver<VectorType,
225  SolverType,
226  MatrixType,
227  PreconditionerType>>
229 
233  SmartPointer<const MatrixType,
234  MGCoarseGridIterativeSolver<VectorType,
235  SolverType,
236  MatrixType,
237  PreconditionerType>>
239 
243  SmartPointer<const PreconditionerType,
244  MGCoarseGridIterativeSolver<VectorType,
245  SolverType,
246  MatrixType,
247  PreconditionerType>>
249 };
250 
251 
252 
263 template <typename number = double, class VectorType = Vector<number>>
264 class MGCoarseGridHouseholder : public MGCoarseGridBase<VectorType>
265 {
266 public:
270  MGCoarseGridHouseholder(const FullMatrix<number> *A = nullptr);
271 
275  void
276  initialize(const FullMatrix<number> &A);
277 
278  void
279  operator()(const unsigned int level,
280  VectorType & dst,
281  const VectorType & src) const override;
282 
283 private:
288 };
289 
298 template <typename number = double, class VectorType = Vector<number>>
299 class MGCoarseGridSVD : public MGCoarseGridBase<VectorType>
300 {
301 public:
305  MGCoarseGridSVD() = default;
306 
310  void
311  initialize(const FullMatrix<number> &A, const double threshold = 0);
312 
313  void
314  operator()(const unsigned int level,
315  VectorType & dst,
316  const VectorType & src) const;
317 
321  void
322  log() const;
323 
324 private:
329 };
330 
333 #ifndef DOXYGEN
334 /* ------------------ Functions for MGCoarseGridApplySmoother -----------*/
335 template <class VectorType>
337  : coarse_smooth(nullptr)
338 {}
339 
340 template <class VectorType>
342  const MGSmootherBase<VectorType> &coarse_smooth)
343  : coarse_smooth(nullptr)
344 {
345  initialize(coarse_smooth);
346 }
347 
348 
349 template <class VectorType>
350 void
352  const MGSmootherBase<VectorType> &coarse_smooth_)
353 {
354  coarse_smooth =
356  MGCoarseGridApplySmoother<VectorType>>(&coarse_smooth_,
357  typeid(*this).name());
358 }
359 
360 
361 template <class VectorType>
362 void
364 {
365  coarse_smooth = nullptr;
366 }
367 
368 
369 template <class VectorType>
370 void
372  VectorType & dst,
373  const VectorType & src) const
374 {
375  coarse_smooth->smooth(level, dst, src);
376 }
377 
378 /* ------------------ Functions for MGCoarseGridLACIteration ------------ */
379 
380 
381 template <typename SolverType, class VectorType>
383  : solver(0, typeid(*this).name())
384  , matrix(0)
385  , precondition(0)
386 {}
387 
388 
389 template <typename SolverType, class VectorType>
390 template <typename MatrixType, typename PreconditionerType>
392  SolverType & s,
393  const MatrixType & m,
394  const PreconditionerType &p)
395  : solver(&s, typeid(*this).name())
396 {
397  // Workaround: Unfortunately, not every "m" object has a rich enough
398  // interface to populate reinit_(domain|range)_vector. Thus, supply an
399  // empty LinearOperator exemplar.
400  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
401  precondition = linear_operator<VectorType>(matrix, p);
402 }
403 
404 
405 template <typename SolverType, class VectorType>
407 {
408  clear();
409 }
410 
411 
412 template <typename SolverType, class VectorType>
413 template <typename MatrixType, typename PreconditionerType>
414 void
416  SolverType & s,
417  const MatrixType & m,
418  const PreconditionerType &p)
419 {
420  solver = &s;
421  // Workaround: Unfortunately, not every "m" object has a rich enough
422  // interface to populate reinit_(domain|range)_vector. Thus, supply an
423  // empty LinearOperator exemplar.
424  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
425  precondition = linear_operator<VectorType>(matrix, p);
426 }
427 
428 
429 template <typename SolverType, class VectorType>
430 void
432 {
433  solver = nullptr;
435  precondition = LinearOperator<VectorType>();
436 }
437 
438 
439 template <typename SolverType, class VectorType>
440 void
442 operator()(const unsigned int /* level */,
443  VectorType & dst,
444  const VectorType &src) const
445 {
446  Assert(solver != nullptr, ExcNotInitialized());
447  solver->solve(matrix, dst, src, precondition);
448 }
449 
450 
451 template <typename SolverType, class VectorType>
452 template <typename MatrixType>
453 void
455  const MatrixType &m)
456 {
457  // Workaround: Unfortunately, not every "m" object has a rich enough
458  // interface to populate reinit_(domain|range)_vector. Thus, supply an
459  // empty LinearOperator exemplar.
460  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
461 }
462 
463 
464 
465 /* ------------------ Functions for MGCoarseGridIterativeSolver ------------ */
466 
467 template <class VectorType,
468  class SolverType,
469  class MatrixType,
470  class PreconditionerType>
471 MGCoarseGridIterativeSolver<VectorType,
472  SolverType,
473  MatrixType,
474  PreconditionerType>::MGCoarseGridIterativeSolver()
475  : solver(0, typeid(*this).name())
476  , matrix(0, typeid(*this).name())
477  , preconditioner(0, typeid(*this).name())
478 {}
479 
480 
481 
482 template <class VectorType,
483  class SolverType,
484  class MatrixType,
485  class PreconditionerType>
486 MGCoarseGridIterativeSolver<VectorType,
487  SolverType,
488  MatrixType,
489  PreconditionerType>::
490  MGCoarseGridIterativeSolver(SolverType & solver,
491  const MatrixType & matrix,
492  const PreconditionerType &preconditioner)
493  : solver(&solver, typeid(*this).name())
494  , matrix(&matrix, typeid(*this).name())
495  , preconditioner(&preconditioner, typeid(*this).name())
496 {}
497 
498 
499 
500 template <class VectorType,
501  class SolverType,
502  class MatrixType,
503  class PreconditionerType>
504 void
506  VectorType,
507  SolverType,
508  MatrixType,
509  PreconditionerType>::initialize(SolverType & solver_,
510  const MatrixType & matrix_,
511  const PreconditionerType &preconditioner_)
512 {
513  solver = &solver_;
514  matrix = &matrix_;
515  preconditioner = &preconditioner_;
516 }
517 
518 
519 
520 template <class VectorType,
521  class SolverType,
522  class MatrixType,
523  class PreconditionerType>
524 void
525 MGCoarseGridIterativeSolver<VectorType,
526  SolverType,
527  MatrixType,
528  PreconditionerType>::clear()
529 {
530  solver = 0;
531  matrix = 0;
532  preconditioner = 0;
533 }
534 
535 
536 
537 template <class VectorType,
538  class SolverType,
539  class MatrixType,
540  class PreconditionerType>
541 void
542 MGCoarseGridIterativeSolver<VectorType,
543  SolverType,
544  MatrixType,
545  PreconditionerType>::
546 operator()(const unsigned int /*level*/,
547  VectorType & dst,
548  const VectorType &src) const
549 {
550  Assert(solver != nullptr, ExcNotInitialized());
551  Assert(matrix != nullptr, ExcNotInitialized());
552  Assert(preconditioner != nullptr, ExcNotInitialized());
553  solver->solve(*matrix, dst, src, *preconditioner);
554 }
555 
556 
557 
558 /* ------------------ Functions for MGCoarseGridHouseholder ------------ */
559 
560 template <typename number, class VectorType>
562  const FullMatrix<number> *A)
563 {
564  if (A != nullptr)
565  householder.initialize(*A);
566 }
567 
568 
569 
570 template <typename number, class VectorType>
571 void
573  const FullMatrix<number> &A)
574 {
575  householder.initialize(A);
576 }
577 
578 
579 
580 template <typename number, class VectorType>
581 void
583 operator()(const unsigned int /*level*/,
584  VectorType & dst,
585  const VectorType &src) const
586 {
587  householder.least_squares(dst, src);
588 }
589 
590 //---------------------------------------------------------------------------
591 
592 
593 
594 template <typename number, class VectorType>
595 void
597  double threshold)
598 {
599  matrix.reinit(A.n_rows(), A.n_cols());
600  matrix = A;
601  matrix.compute_inverse_svd(threshold);
602 }
603 
604 
605 template <typename number, class VectorType>
606 void
607 MGCoarseGridSVD<number, VectorType>::operator()(const unsigned int /*level*/,
608  VectorType & dst,
609  const VectorType &src) const
610 {
611  matrix.vmult(dst, src);
612 }
613 
614 
615 template <typename number, class VectorType>
616 void
618 {
619  const unsigned int n = std::min(matrix.n_rows(), matrix.n_cols());
620 
621  for (unsigned int i = 0; i < n; ++i)
622  deallog << ' ' << matrix.singular_value(i);
623  deallog << std::endl;
624 }
625 
626 
627 #endif // DOXYGEN
628 
629 DEAL_II_NAMESPACE_CLOSE
630 
631 #endif
virtual void operator()(const unsigned int level, VectorType &dst, const VectorType &src) const =0
void initialize(const FullMatrix< number > &A, const double threshold=0)
Contents is actually a matrix.
void operator()(const unsigned int level, VectorType &dst, const VectorType &src) const override
LinearOperator< VectorType > precondition
Definition: mg_coarse.h:164
static ::ExceptionBase & ExcNotInitialized()
virtual void operator()(const unsigned int level, VectorType &dst, const VectorType &src) const override
SmartPointer< SolverType, MGCoarseGridIterativeSolver< VectorType, SolverType, MatrixType, PreconditionerType > > solver
Definition: mg_coarse.h:228
void operator()(const unsigned int level, VectorType &dst, const VectorType &src) const
void initialize(SolverType &, const MatrixType &, const PreconditionerType &)
LAPACKFullMatrix< number > matrix
Definition: mg_coarse.h:328
SmartPointer< const MatrixType, MGCoarseGridIterativeSolver< VectorType, SolverType, MatrixType, PreconditionerType > > matrix
Definition: mg_coarse.h:238
#define Assert(cond, exc)
Definition: exceptions.h:1407
MGCoarseGridSVD()=default
void initialize(SolverType &solver, const MatrixType &matrix, const PreconditionerType &precondition)
LinearOperator< VectorType > matrix
Definition: mg_coarse.h:159
SmartPointer< SolverType, MGCoarseGridLACIteration< SolverType, VectorType > > solver
Definition: mg_coarse.h:154
void initialize(const MGSmootherBase< VectorType > &coarse_smooth)
MGCoarseGridHouseholder(const FullMatrix< number > *A=nullptr)
Householder< number > householder
Definition: mg_coarse.h:287
void operator()(const unsigned int level, VectorType &dst, const VectorType &src) const
SmartPointer< const PreconditionerType, MGCoarseGridIterativeSolver< VectorType, SolverType, MatrixType, PreconditionerType > > preconditioner
Definition: mg_coarse.h:248
SmartPointer< const MGSmootherBase< VectorType >, MGCoarseGridApplySmoother< VectorType > > coarse_smooth
Definition: mg_coarse.h:78
void operator()(const unsigned int level, VectorType &dst, const VectorType &src) const override
void initialize(const FullMatrix< number > &A)
void set_matrix(const MatrixType &)
void log() const