1#ifndef DUNE_FEM_GMRES_HH
2#define DUNE_FEM_GMRES_HH
24 Matrix(
int n,
int m) : n_(n), m_(m), data_( n*m, T(0) ) {}
25 Matrix(
int n) : Matrix(n,n) {}
28 T& operator()(
int i,
int j)
30 assert(i>=0 && i<n_ && j>=0 && j<m_);
31 return data_[i*m_ + j];
35 T operator()(
int i,
int j)
const
37 assert(i>=0 && i<n_ && j>=0 && j<m_);
38 return data_[i*m_ + j];
42 operator T*() {
return data_.data(); }
43 operator const T *()
const {
return data_.data(); }
47 std::vector< T > data_;
52 template <
class FieldType>
53 FieldType
scalarProduct(
const int dim,
const FieldType *x,
const FieldType* y )
56 for(
int i=0; i<dim; ++i )
58 scp += x[ i ] * y[ i ];
64 template <
class Communication,
class FieldType,
class DiscreteFunction>
65 void gemv(
const Communication& comm,
67 std::vector< DiscreteFunction >& v,
68 const DiscreteFunction& vjp,
72 for(
int l=0; l<m; ++l)
77 const auto& auxiliaryDofs = vjp.space().auxiliaryDofs();
78 const auto& vj = vjp.dofVector();
80 const size_t numAuxiliarys = auxiliaryDofs.size();
81 for(
size_t auxiliary = 0, i = 0 ; auxiliary < numAuxiliarys; ++auxiliary )
83 const size_t nextAuxiliary = auxiliaryDofs[ auxiliary ];
84 for(; i < nextAuxiliary; ++i )
86 for(
int l=0; l<m; ++l)
88 y[ l ] += (vj[ i ] * v[ l ].dofVector()[ i ]);
98 template<
class FieldType>
100 FieldType* x, FieldType* y,
101 const FieldType c,
const FieldType s)
106 const FieldType _x = *x;
107 const FieldType _y = *y;
119 template <
class Operator,
class Preconditioner,
class DiscreteFunction>
121 std::vector< DiscreteFunction >& v,
123 const DiscreteFunction& b,
125 const double tolerance,
126 const int maxIterations,
127 const int toleranceCriteria,
128 std::ostream* os =
nullptr )
130 typedef typename DiscreteFunction :: RangeFieldType FieldType;
132 const auto& comm = u.space().gridPart().comm();
134 detail::Matrix< FieldType > H( m+1, m );
135 std::vector< FieldType > g_( 6*m, 0.0 );
137 FieldType* g = g_.data();
138 FieldType* s = g + (m+1);
139 FieldType* c = s + m;
140 FieldType* y = c + m;
142 DiscreteFunction& v0 = v[ 0 ];
144 std::vector< FieldType > global_dot( m+1, FieldType(0) );
147 double _tolerance = tolerance;
150 global_dot[ 0 ] = b.scalarProductDofs( b );
167 global_dot[ 0 ] = v0.scalarProductDofs( v0 );
170 FieldType res =
std::sqrt(global_dot[0]);
179 (*os) <<
"Fem::GMRES outer iteration : " << res << std::endl;
182 if (res < _tolerance)
break;
185 for(
int i=1; i<=m; i++) g[i] = 0.0;
193 for(
int j=0; j<m; j++)
195 DiscreteFunction& vj = v[ j ];
196 DiscreteFunction& vjp = v[ j + 1 ];
202 DiscreteFunction& z = v[ m+1 ];
203 (*preconditioner)(vj, z );
214 gemv(comm, j+1, v, vjp, global_dot.data());
216 for(
int i=0; i<=j; i++) H(i,j) = global_dot[i];
223 for(
int l=0; l<j+1; ++l)
225 vjp.axpy( -global_dot[l], v[l] );
228 global_dot[ 0 ] = vjp.scalarProductDofs( vjp );
237 for(
int i=0; i<j; i++)
239 rotate(1, &H(i+1,j), &H(i,j), c[i], s[i]);
242 const FieldType h_j_j = H(j,j);
243 const FieldType h_jp_j = H(j+1,j);
244 const FieldType norm =
std::sqrt(h_j_j*h_j_j + h_jp_j*h_jp_j);
246 s[j] = -h_jp_j / norm;
247 rotate(1, &H(j+1,j), &H(j,j), c[j], s[j]);
248 rotate(1, &g[j+1], &g[j], c[j], s[j]);
252 (*os) <<
"Fem::GMRES it: " << iterations <<
" : " <<
std::abs(g[j+1]) << std::endl;
257 || iterations >= maxIterations )
break;
264 int last = iterations%m;
265 if (last == 0) last = m;
268 for(
int i=last-1; i>=0; --i)
270 const FieldType dot =
scalarProduct( last-(i+1), &H(i,i)+1, &y[i+1] );
271 y[i] = (g[i] - dot)/ H(i,i);
278 DiscreteFunction& u_tmp = v[ m ];
279 DiscreteFunction& z = v[ m+1 ];
283 for(
int i=0; i<last; ++i)
285 u_tmp.axpy( y[ i ], v[ i ] );
288 (*preconditioner)(u_tmp, z);
293 for(
int i=0; i<last; ++i)
295 u.axpy( y[ i ], v[ i ] );
299 if (
std::abs(g[last]) < _tolerance)
break;
304 (*os) <<
"Fem::GMRES: number of iterations: "
309 return (iterations < maxIterations) ? iterations : -iterations;
double sqrt(const Dune::Fem::Double &v)
Definition: double.hh:977
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
void gemv(const Communication &comm, const int m, std::vector< DiscreteFunction > &v, const DiscreteFunction &vjp, FieldType *y)
Definition: gmres.hh:65
int gmres(Operator &op, Preconditioner *preconditioner, std::vector< DiscreteFunction > &v, DiscreteFunction &u, const DiscreteFunction &b, const int m, const double tolerance, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: gmres.hh:120
void rotate(const int dim, FieldType *x, FieldType *y, const FieldType c, const FieldType s)
dblas_rotate with inc = 1
Definition: gmres.hh:99
FieldType scalarProduct(const int dim, const FieldType *x, const FieldType *y)
return x * y
Definition: gmres.hh:53
abstract operator
Definition: operator.hh:34
static const int relative
Definition: cg.hh:18
static const int residualReduction
Definition: cg.hh:19