MFEM  v3.1
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 // Implementation of data types for sparse matrix smoothers
13 
14 #include <iostream>
15 #include "vector.hpp"
16 #include "matrix.hpp"
17 #include "sparsemat.hpp"
18 #include "sparsesmoothers.hpp"
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 
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 
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 
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:1653
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Definition: operator.hpp:41
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:106
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
Definition: operator.hpp:35
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const
Definition: sparsemat.cpp:1622
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with Jacobi smoother.
void mfem_error(const char *msg)
Definition: error.cpp:23
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:1514
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
Definition: sparsemat.cpp:1439
DSmoother(int t=0, double s=1., int it=1)
Create Jacobi smoother.
Vector data type.
Definition: vector.hpp:33
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:1685
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:1709
Abstract operator.
Definition: operator.hpp:21
const SparseMatrix * oper