Reference documentation for deal.II version Git 92c1895bf8 2020-06-04 17:48:28 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mg_coarse.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2020 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/base/config.h>
21 
25 
27 
29 
32 
39 template <class VectorType = Vector<double>>
40 class MGCoarseGridApplySmoother : public MGCoarseGridBase<VectorType>
41 {
42 public:
47 
52 
56  void
57  clear();
58 
62  void
63  initialize(const MGSmootherBase<VectorType> &coarse_smooth);
64 
68  void
69  operator()(const unsigned int level,
70  VectorType & dst,
71  const VectorType & src) const override;
72 
73 private:
80 };
81 
82 
96 template <typename SolverType, class VectorType = Vector<double>>
98  : public MGCoarseGridBase<VectorType>
99 {
100 public:
105 
110  template <typename MatrixType, typename PreconditionerType>
112  const MatrixType &,
113  const PreconditionerType &);
114 
119 
123  template <typename MatrixType, typename PreconditionerType>
124  void
125  initialize(SolverType &, const MatrixType &, const PreconditionerType &);
126 
130  void
131  clear();
132 
137  void
138  operator()(const unsigned int level,
139  VectorType & dst,
140  const VectorType & src) const;
141 
146  template <typename MatrixType>
147  void
148  set_matrix(const MatrixType &);
149 
150 private:
156 
161 
166 };
167 
168 
169 
176 template <class VectorType,
177  class SolverType,
178  class MatrixType,
179  class PreconditionerType>
181 {
182 public:
187 
192  MGCoarseGridIterativeSolver(SolverType & solver,
193  const MatrixType & matrix,
194  const PreconditionerType &precondition);
195 
200  void
201  initialize(SolverType & solver,
202  const MatrixType & matrix,
203  const PreconditionerType &precondition);
204 
208  void
209  clear();
210 
215  virtual void
216  operator()(const unsigned int level,
217  VectorType & dst,
218  const VectorType & src) const override;
219 
220 private:
224  SmartPointer<SolverType,
226  SolverType,
227  MatrixType,
228  PreconditionerType>>
230 
234  SmartPointer<const MatrixType,
235  MGCoarseGridIterativeSolver<VectorType,
236  SolverType,
237  MatrixType,
238  PreconditionerType>>
240 
244  SmartPointer<const PreconditionerType,
245  MGCoarseGridIterativeSolver<VectorType,
246  SolverType,
247  MatrixType,
248  PreconditionerType>>
250 };
251 
252 
253 
264 template <typename number = double, class VectorType = Vector<number>>
265 class MGCoarseGridHouseholder : public MGCoarseGridBase<VectorType>
266 {
267 public:
271  MGCoarseGridHouseholder(const FullMatrix<number> *A = nullptr);
272 
276  void
278 
279  void
280  operator()(const unsigned int level,
281  VectorType & dst,
282  const VectorType & src) const override;
283 
284 private:
289 };
290 
299 template <typename number = double, class VectorType = Vector<number>>
300 class MGCoarseGridSVD : public MGCoarseGridBase<VectorType>
301 {
302 public:
306  MGCoarseGridSVD() = default;
307 
311  void
312  initialize(const FullMatrix<number> &A, const double threshold = 0);
313 
314  void
315  operator()(const unsigned int level,
316  VectorType & dst,
317  const VectorType & src) const;
318 
322  void
323  log() const;
324 
325 private:
330 };
331 
334 #ifndef DOXYGEN
335 /* ------------------ Functions for MGCoarseGridApplySmoother -----------*/
336 template <class VectorType>
338  : coarse_smooth(nullptr)
339 {}
340 
341 template <class VectorType>
343  const MGSmootherBase<VectorType> &coarse_smooth)
344  : coarse_smooth(nullptr)
345 {
346  initialize(coarse_smooth);
347 }
348 
349 
350 template <class VectorType>
351 void
353  const MGSmootherBase<VectorType> &coarse_smooth_)
354 {
355  coarse_smooth =
357  MGCoarseGridApplySmoother<VectorType>>(&coarse_smooth_,
358  typeid(*this).name());
359 }
360 
361 
362 template <class VectorType>
363 void
365 {
366  coarse_smooth = nullptr;
367 }
368 
369 
370 template <class VectorType>
371 void
373  VectorType & dst,
374  const VectorType & src) const
375 {
376  coarse_smooth->smooth(level, dst, src);
377 }
378 
379 /* ------------------ Functions for MGCoarseGridLACIteration ------------ */
380 
381 
382 template <typename SolverType, class VectorType>
384  : solver(0, typeid(*this).name())
385  , matrix(0)
386  , precondition(0)
387 {}
388 
389 
390 template <typename SolverType, class VectorType>
391 template <typename MatrixType, typename PreconditionerType>
393  SolverType & s,
394  const MatrixType & m,
395  const PreconditionerType &p)
396  : solver(&s, typeid(*this).name())
397 {
398  // Workaround: Unfortunately, not every "m" object has a rich enough
399  // interface to populate reinit_(domain|range)_vector. Thus, supply an
400  // empty LinearOperator exemplar.
401  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
402  precondition = linear_operator<VectorType>(matrix, p);
403 }
404 
405 
406 template <typename SolverType, class VectorType>
408 {
409  clear();
410 }
411 
412 
413 template <typename SolverType, class VectorType>
414 template <typename MatrixType, typename PreconditionerType>
415 void
417  SolverType & s,
418  const MatrixType & m,
419  const PreconditionerType &p)
420 {
421  solver = &s;
422  // Workaround: Unfortunately, not every "m" object has a rich enough
423  // interface to populate reinit_(domain|range)_vector. Thus, supply an
424  // empty LinearOperator exemplar.
425  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
426  precondition = linear_operator<VectorType>(matrix, p);
427 }
428 
429 
430 template <typename SolverType, class VectorType>
431 void
433 {
434  solver = nullptr;
436  precondition = LinearOperator<VectorType>();
437 }
438 
439 
440 template <typename SolverType, class VectorType>
441 void
443 operator()(const unsigned int /* level */,
444  VectorType & dst,
445  const VectorType &src) const
446 {
447  Assert(solver != nullptr, ExcNotInitialized());
448  solver->solve(matrix, dst, src, precondition);
449 }
450 
451 
452 template <typename SolverType, class VectorType>
453 template <typename MatrixType>
454 void
456  const MatrixType &m)
457 {
458  // Workaround: Unfortunately, not every "m" object has a rich enough
459  // interface to populate reinit_(domain|range)_vector. Thus, supply an
460  // empty LinearOperator exemplar.
461  matrix = linear_operator<VectorType>(LinearOperator<VectorType>(), m);
462 }
463 
464 
465 
466 /* ------------------ Functions for MGCoarseGridIterativeSolver ------------ */
467 
468 template <class VectorType,
469  class SolverType,
470  class MatrixType,
471  class PreconditionerType>
473  SolverType,
474  MatrixType,
475  PreconditionerType>::MGCoarseGridIterativeSolver()
476  : solver(0, typeid(*this).name())
477  , matrix(0, typeid(*this).name())
478  , preconditioner(0, typeid(*this).name())
479 {}
480 
481 
482 
483 template <class VectorType,
484  class SolverType,
485  class MatrixType,
486  class PreconditionerType>
488  SolverType,
489  MatrixType,
490  PreconditionerType>::
491  MGCoarseGridIterativeSolver(SolverType & solver,
492  const MatrixType & matrix,
493  const PreconditionerType &preconditioner)
494  : solver(&solver, typeid(*this).name())
495  , matrix(&matrix, typeid(*this).name())
496  , preconditioner(&preconditioner, typeid(*this).name())
497 {}
498 
499 
500 
501 template <class VectorType,
502  class SolverType,
503  class MatrixType,
504  class PreconditionerType>
505 void
507  VectorType,
508  SolverType,
509  MatrixType,
510  PreconditionerType>::initialize(SolverType & solver_,
511  const MatrixType & matrix_,
512  const PreconditionerType &preconditioner_)
513 {
514  solver = &solver_;
515  matrix = &matrix_;
516  preconditioner = &preconditioner_;
517 }
518 
519 
520 
521 template <class VectorType,
522  class SolverType,
523  class MatrixType,
524  class PreconditionerType>
525 void
527  SolverType,
528  MatrixType,
529  PreconditionerType>::clear()
530 {
531  solver = 0;
532  matrix = 0;
533  preconditioner = 0;
534 }
535 
536 
537 
538 template <class VectorType,
539  class SolverType,
540  class MatrixType,
541  class PreconditionerType>
542 void
544  SolverType,
545  MatrixType,
546  PreconditionerType>::
547 operator()(const unsigned int /*level*/,
548  VectorType & dst,
549  const VectorType &src) const
550 {
551  Assert(solver != nullptr, ExcNotInitialized());
552  Assert(matrix != nullptr, ExcNotInitialized());
553  Assert(preconditioner != nullptr, ExcNotInitialized());
554  solver->solve(*matrix, dst, src, *preconditioner);
555 }
556 
557 
558 
559 /* ------------------ Functions for MGCoarseGridHouseholder ------------ */
560 
561 template <typename number, class VectorType>
563  const FullMatrix<number> *A)
564 {
565  if (A != nullptr)
566  householder.initialize(*A);
567 }
568 
569 
570 
571 template <typename number, class VectorType>
572 void
574  const FullMatrix<number> &A)
575 {
576  householder.initialize(A);
577 }
578 
579 
580 
581 template <typename number, class VectorType>
582 void
584 operator()(const unsigned int /*level*/,
585  VectorType & dst,
586  const VectorType &src) const
587 {
588  householder.least_squares(dst, src);
589 }
590 
591 //---------------------------------------------------------------------------
592 
593 
594 
595 template <typename number, class VectorType>
596 void
598  double threshold)
599 {
600  matrix.reinit(A.n_rows(), A.n_cols());
601  matrix = A;
602  matrix.compute_inverse_svd(threshold);
603 }
604 
605 
606 template <typename number, class VectorType>
607 void
608 MGCoarseGridSVD<number, VectorType>::operator()(const unsigned int /*level*/,
609  VectorType & dst,
610  const VectorType &src) const
611 {
612  matrix.vmult(dst, src);
613 }
614 
615 
616 template <typename number, class VectorType>
617 void
619 {
620  const unsigned int n = std::min(matrix.n_rows(), matrix.n_cols());
621 
622  for (unsigned int i = 0; i < n; ++i)
623  deallog << ' ' << matrix.singular_value(i);
624  deallog << std::endl;
625 }
626 
627 
628 #endif // DOXYGEN
629 
631 
632 #endif
LogStream deallog
Definition: logstream.cc:37
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:165
static ::ExceptionBase & ExcNotInitialized()
SmartPointer< SolverType, MGCoarseGridIterativeSolver< VectorType, SolverType, MatrixType, PreconditionerType > > solver
Definition: mg_coarse.h:229
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:329
SmartPointer< const MatrixType, MGCoarseGridIterativeSolver< VectorType, SolverType, MatrixType, PreconditionerType > > matrix
Definition: mg_coarse.h:239
#define Assert(cond, exc)
Definition: exceptions.h:1408
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:361
unsigned int level
Definition: grid_out.cc:4355
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const =0
LinearOperator< VectorType > matrix
Definition: mg_coarse.h:160
SmartPointer< SolverType, MGCoarseGridLACIteration< SolverType, VectorType > > solver
Definition: mg_coarse.h:155
static const char A
void initialize(const MGSmootherBase< VectorType > &coarse_smooth)
MGCoarseGridHouseholder(const FullMatrix< number > *A=nullptr)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
T min(const T &t, const MPI_Comm &mpi_communicator)
Householder< number > householder
Definition: mg_coarse.h:288
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:249
SmartPointer< const MGSmootherBase< VectorType >, MGCoarseGridApplySmoother< VectorType > > coarse_smooth
Definition: mg_coarse.h:79
void operator()(const unsigned int level, VectorType &dst, const VectorType &src) const override
void initialize(const FullMatrix< number > &A)
#define DEAL_II_DEPRECATED
Definition: config.h:150
void set_matrix(const MatrixType &)
void log() const