MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
sparsesmoothers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 // Implementation of data types for sparse matrix smoothers
13 
14 #include "vector.hpp"
15 #include "matrix.hpp"
16 #include "sparsemat.hpp"
17 #include "sparsesmoothers.hpp"
18 #include <iostream>
19 
20 namespace mfem
21 {
22 
24 {
25  oper = dynamic_cast<const SparseMatrix*>(&a);
26  if (oper == NULL)
27  {
28  mfem_error("SparseSmoother::SetOperator : not a SparseMatrix!");
29  }
30  height = oper->Height();
31  width = oper->Width();
32 }
33 
34 /// Matrix vector multiplication with GS Smoother.
35 void GSSmoother::Mult(const Vector &x, Vector &y) const
36 {
37  if (!iterative_mode)
38  {
39  y = 0.0;
40  }
41  for (int i = 0; i < iterations; i++)
42  {
43  if (type != 2)
44  {
45  oper->Gauss_Seidel_forw(x, y);
46  }
47  if (type != 1)
48  {
49  oper->Gauss_Seidel_back(x, y);
50  }
51  }
52 }
53 
54 /// Create the Jacobi smoother.
55 DSmoother::DSmoother(const SparseMatrix &a, int t, double s, int it)
56  : SparseSmoother(a)
57 {
58  type = t;
59  scale = s;
60  iterations = it;
61 }
62 
63 /// Matrix vector multiplication with Jacobi smoother.
64 void DSmoother::Mult(const Vector &x, Vector &y) const
65 {
66  if (!iterative_mode && type == 0 && iterations == 1)
67  {
68  oper->DiagScale(x, y, scale);
69  return;
70  }
71 
72  z.SetSize(width);
73 
74  Vector *r = &y, *p = &z;
75 
76  if (iterations % 2 == 0)
77  {
78  Swap<Vector*>(r, p);
79  }
80 
81  if (!iterative_mode)
82  {
83  *p = 0.0;
84  }
85  else if (iterations % 2)
86  {
87  *p = y;
88  }
89  for (int i = 0; i < iterations; i++)
90  {
91  if (type == 0)
92  {
93  oper->Jacobi(x, *p, *r, scale);
94  }
95  else if (type == 1)
96  {
97  oper->Jacobi2(x, *p, *r, scale);
98  }
99  else if (type == 2)
100  {
101  oper->Jacobi3(x, *p, *r, scale);
102  }
103  else
104  {
105  mfem_error("DSmoother::Mult wrong type");
106  }
107  Swap<Vector*>(r, p);
108  }
109 }
110 
111 }
void DiagScale(const Vector &b, Vector &x, double sc=1.0) const
Definition: sparsemat.cpp:2144
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
virtual void SetOperator(const Operator &a)
Set/update the solver for the given operator.
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:504
Data type sparse matrix.
Definition: sparsemat.hpp:40
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const
Definition: sparsemat.cpp:2113
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with Jacobi smoother.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with GS Smoother.
void Gauss_Seidel_back(const Vector &x, Vector &y) const
Definition: sparsemat.cpp:1997
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
Definition: sparsemat.cpp:1913
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
double a
Definition: lissajous.cpp:41
DSmoother(int t=0, double s=1., int it=1)
Create Jacobi smoother.
Vector data type.
Definition: vector.hpp:48
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:2176
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:2200
Abstract operator.
Definition: operator.hpp:24
const SparseMatrix * oper
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28