MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
sparsesmoothers.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
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
20namespace 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.
35void 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 {
46 }
47 if (type != 1)
48 {
50 }
51 }
52}
53
54/// Create the Jacobi smoother.
57{
58 type = t;
59 scale = s;
60 iterations = it;
61}
62
63/// Matrix vector multiplication with Jacobi smoother.
64void DSmoother::Mult(const Vector &x, Vector &y) const
65{
66 if (!iterative_mode && type == 0 && iterations == 1)
67 {
69 return;
70 }
71
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, use_abs_diag);
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}
DSmoother(int t=0, real_t s=1., int it=1)
Create Jacobi smoother.
bool use_abs_diag
Uses abs values of the diagonal entries. Relevant only when type = 0.
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.
Abstract operator.
Definition operator.hpp:25
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
Data type sparse matrix.
Definition sparsemat.hpp:51
void Gauss_Seidel_back(const Vector &x, Vector &y) const
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, real_t sc=1.0) const
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, real_t sc, bool use_abs_diag=false) const
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
void DiagScale(const Vector &b, Vector &x, real_t sc=1.0, bool use_abs_diag=false) const
x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, real_t sc=1.0) const
virtual void SetOperator(const Operator &a)
Set/update the solver for the given operator.
const SparseMatrix * oper
Vector data type.
Definition vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t a
Definition lissajous.cpp:41
void Swap(Array< T > &, Array< T > &)
Definition array.hpp:648
void mfem_error(const char *msg)
Definition error.cpp:154
float real_t
Definition config.hpp:43
real_t p(const Vector &x, real_t t)
RefCoord t[3]
RefCoord s[3]