ug4
vanka.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3  * Author: Christian Wehner
4  *
5  * This file is part of UG4.
6  *
7  * UG4 is free software: you can redistribute it and/or modify it under the
8  * terms of the GNU Lesser General Public License version 3 (as published by the
9  * Free Software Foundation) with the following additional attribution
10  * requirements (according to LGPL/GPL v3 §7):
11  *
12  * (1) The following notice must be displayed in the Appropriate Legal Notices
13  * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14  *
15  * (2) The following notice must be displayed at a prominent place in the
16  * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17  *
18  * (3) The following bibliography is recommended for citation and must be
19  * preserved in all covered files:
20  * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21  * parallel geometric multigrid solver on hierarchically distributed grids.
22  * Computing and visualization in science 16, 4 (2013), 151-164"
23  * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24  * flexible software system for simulating pde based models on high performance
25  * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26  *
27  * This program is distributed in the hope that it will be useful,
28  * but WITHOUT ANY WARRANTY; without even the implied warranty of
29  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  * GNU Lesser General Public License for more details.
31  */
32 
33 #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__VANKA__
34 #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__VANKA__
35 
38 
39 #ifdef UG_PARALLEL
40  #include "pcl/pcl_util.h"
43 #endif
44 
45 namespace ug{
46 
47 static const size_t MAXBLOCKSIZE = 53;
48 
49 template<typename Matrix_type, typename Vector_type>
50 bool Vanka_step(const Matrix_type &A, Vector_type &x, const Vector_type &b, number relax)
51 {
55 
56  size_t blockind[MAXBLOCKSIZE];
57 
58  for(size_t i=0; i < x.size(); i++)
59  {
60  x[i]=0;
61  };
62 
63  for(size_t i=0; i < x.size(); i++)
64  {
65  if (A(i,i)!=0){
66  /* // do usual gauss-seidel (would be needed in case of Dirichlet pressure condition)
67  typename Vector_type::value_type def = b[i];
68 
69  for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) && it.index() < i; ++it)
70  // s -= it.value() * x[it.index()];
71  MatMultAdd(def, 1.0, def, -1.0, it.value(), x[it.index()]);
72  // x[i] = s/A(i,i)
73  InverseMatMult(x[i], relax, A(i,i), s);*/
74  continue;
75  };
76 
77  size_t blocksize=0;
78 
79  for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
80  blockind[blocksize] = it.index();
81  blocksize++;
82  };
83  if (blocksize>MAXBLOCKSIZE) UG_THROW("MAXBLOCKSIZE too small\n");
84  mat.resize(blocksize,blocksize);
85  s.resize(blocksize);
86  localx.resize(blocksize);
87  for (size_t j=0;j<blocksize;j++){
88  // fill local block matrix
89  for (size_t k=0;k<blocksize;k++){
90  mat.subassign(j,k,A(blockind[j],blockind[k]));
91  };
92  // compute rhs
93  typename Vector_type::value_type sj = b[blockind[j]];
94  for(typename Matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
95  MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
96  };
97  s.subassign(j,sj);
98  };
99  // solve block
100  InverseMatMult(localx,1,mat,s);
101  for (size_t j=0;j<blocksize;j++){
102  x[blockind[j]] += relax*localx[j];
103  };
104  }
105 
106  return true;
107 }
108 
109 // Diagonal Vanka block smoother:
110 // When setting up the local block matrix the side-diagonal entries are left away, except for the pressure.
111 // The local block matrix therefore has the form
112 //
113 // x 0 ... 0 x
114 // 0 x 0 ... 0 x
115 // 0 0 x 0 ... x
116 // ...
117 // 0 ... 0 x x
118 // x x ... x x x
119 //
120 // The velocity off-diagonal entries are considered in the local defect vector.
121 // The local block system can be solved in O(n) time so that a step of this smoother is computationly cheaper than the
122 // full Vanka smoother.
123 template<typename Matrix_type, typename Vector_type>
124 bool Diag_Vanka_step(const Matrix_type &A, Vector_type &x, const Vector_type &b, number relax)
125 {
126  typedef typename Vector_type::value_type vector_block_type;
130  s.resize(MAXBLOCKSIZE);
131  typedef typename Matrix_type::value_type block_type;
132 
133  size_t blockind[MAXBLOCKSIZE];
134 
135  for(size_t i=0; i < x.size(); i++)
136  {
137  x[i]=0;
138  };
139 
140  for(size_t i=0; i < x.size(); i++)
141  {
142  if (A(i,i)!=0){
143  /* // do usual gauss-seidel (would be needed in case of Dirichlet pressure condition)
144  typename Vector_type::value_type def = b[i];
145 
146  for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) && it.index() < i; ++it)
147  // s -= it.value() * x[it.index()];
148  MatMultAdd(def, 1.0, def, -1.0, it.value(), x[it.index()]);
149  // x[i] = s/A(i,i)
150  InverseMatMult(x[i], relax, A(i,i), def);*/
151  continue;
152  };
153 
154  size_t blocksize=0;
155 
156  for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
157  if (it.index()==i) continue;
158  blockind[blocksize] = it.index();
159  s[blocksize] = b[blockind[blocksize]];
160  for(typename Matrix_type::const_row_iterator rowit = A.begin_row(blockind[blocksize]); rowit != A.end_row(blockind[blocksize]) ; ++rowit){
161  if ((rowit.index()==blockind[blocksize])||(rowit.index()==i)) continue;
162  // s[blocksize] -= a_ij*x_j
163  MatMultAdd(s[blocksize], 1.0, s[blocksize], -1.0, rowit.value(), x[rowit.index()]);
164  };
165  blocksize++;
166  };
167  if (blocksize>MAXBLOCKSIZE) UG_THROW("MAXBLOCKSIZE too small\n");
168  // remark: blocksize is without pressure variable, so actual blocksize is blocksize+1
169  block_type a_ii = A(i,i);
170  typename Vector_type::value_type s_i = b[i];
171  // Gauss elimination on local block matrix
172  for (size_t j=0;j<blocksize;j++){
173  block_type a_q = A(i,blockind[j]);
174  block_type a_jj = A(blockind[j],blockind[j]);
175  a_q /= a_jj;
176  // s_i -= a_ij/a_jj*s_j
177  MatMultAdd(s_i, 1.0, s_i, -1.0, a_q, s[j]);
178  // a_ii -= a_ij/a_jj*a_ji
179  a_ii-=a_q*A(blockind[j],i);
180  }
181  // solve diagonalized system
182  // x[i] = s_i/a_ii
183  InverseMatMult(x[i], 1.0, a_ii, s_i);
184  for (size_t j=0;j<blocksize;j++){
185  // s_j-=a_ji*x_i
186  MatMultAdd(s[j], 1.0, s[j], -1.0, A(blockind[j],i), x[i]);
187  // x_j=1/a_jj*s_j
188  InverseMatMult(x[blockind[j]], relax, A(blockind[j],blockind[j]),s[j]);
189  }
190  }
191  return true;
192 }
193 
194 
196 template <typename TAlgebra>
197 class Vanka : public IPreconditioner<TAlgebra>
198 {
199  public:
201  typedef TAlgebra algebra_type;
202 
204  typedef typename TAlgebra::vector_type vector_type;
205 
207  typedef typename TAlgebra::matrix_type matrix_type;
208 
211 
214 
215  protected:
216  using base_type::set_debug;
219 
220  public:
222  Vanka() {m_relax=1;};
223 
226  {
228  newInst->set_debug(debug_writer());
229  newInst->set_damp(this->damping());
230  return newInst;
231  }
232 
234  virtual ~Vanka()
235  {};
236 
237  virtual bool supports_parallel() const {return true;}
238 
240  public:
242 
243  protected:
245 
246  protected:
248  virtual const char* name() const {return "Vanka";}
249 
252  {
253 #ifdef UG_PARALLEL
254  if(pcl::NumProcs() > 1)
255  {
256  // copy original matrix
257  MakeConsistent(*pOp, m_A);
258  // set zero on slaves
259  std::vector<IndexLayout::Element> vIndex;
260  CollectUniqueElements(vIndex, m_A.layouts()->slave());
261  SetDirichletRow(m_A, vIndex);
262  }
263 #endif
264  return true;
265  }
266 
268  {
269 #ifdef UG_PARALLEL
270  if(pcl::NumProcs() > 1)
271  {
272  // make defect unique
273  // todo: change that copying
274  vector_type dhelp;
275  dhelp.resize(d.size()); dhelp = d;
276  dhelp.change_storage_type(PST_UNIQUE);
277 
278  if(!Vanka_step(m_A, c, dhelp, m_relax)) return false;
279 
280  c.set_storage_type(PST_UNIQUE);
281  return true;
282  }
283  else
284 #endif
285  {
286  if(!Vanka_step(*pOp, c, d, m_relax)) return false;
287 
288 #ifdef UG_PARALLEL
289  c.set_storage_type(PST_UNIQUE);
290 #endif
291  return true;
292  }
293  }
294 
296  virtual bool postprocess() {return true;}
297 
298  protected:
299 #ifdef UG_PARALLEL
301 #endif
302 
303 };
304 
306 template <typename TAlgebra>
307 class DiagVanka : public IPreconditioner<TAlgebra>
308 {
309  public:
311  typedef TAlgebra algebra_type;
312 
314  typedef typename TAlgebra::vector_type vector_type;
315 
317  typedef typename TAlgebra::matrix_type matrix_type;
318 
321 
324 
325  protected:
326  using base_type::set_debug;
329 
330  public:
333 
336  {
338  newInst->set_debug(debug_writer());
339  newInst->set_damp(this->damping());
340  return newInst;
341  }
342 
343  virtual bool supports_parallel() const {return true;}
344 
346  virtual ~DiagVanka()
347  {};
348 
350  public:
352 
353  protected:
355 
356  protected:
358  virtual const char* name() const {return "DiagVanka";}
359 
362  {
363 #ifdef UG_PARALLEL
364  if(pcl::NumProcs() > 1)
365  {
366  // copy original matrix
367  MakeConsistent(*pOp, m_A);
368  // set zero on slaves
369  std::vector<IndexLayout::Element> vIndex;
370  CollectUniqueElements(vIndex, m_A.layouts()->slave());
371  SetDirichletRow(m_A, vIndex);
372  }
373 #endif
374  return true;
375  }
376 
378  {
379 #ifdef UG_PARALLEL
380  if(pcl::NumProcs() > 1)
381  {
382  // make defect unique
383  // todo: change that copying
384  vector_type dhelp;
385  dhelp.resize(d.size()); dhelp = d;
386  dhelp.change_storage_type(PST_UNIQUE);
387 
388  if(!Diag_Vanka_step(m_A, c, dhelp, m_relax)) return false;
389 
390  c.set_storage_type(PST_UNIQUE);
391  return true;
392  }
393  else
394 #endif
395  {
396 
397  if(!Diag_Vanka_step(*pOp, c, d, m_relax)) return false;
398 
399 #ifdef UG_PARALLEL
400  c.set_storage_type(PST_UNIQUE);
401 #endif
402  return true;
403  }
404  }
405 
407  virtual bool postprocess() {return true;}
408 
409  protected:
410 #ifdef UG_PARALLEL
412 #endif
413 
414 };
415 
416 
417 } // end namespace ug
418 
419 #endif
parameterString s
Definition: smart_pointer.h:108
virtual void set_debug(SmartPtr< IDebugWriter< algebra_type > > spDebugWriter)
set debug writer
Definition: debug_writer.h:384
void write_debug(const matrix_type &mat, const char *filename)
write debug output for a matrix (if debug writer set)
Definition: debug_writer.h:399
TAlgebra::vector_type vector_type
type of vector
Definition: debug_writer.h:360
SmartPtr< IDebugWriter< algebra_type > > debug_writer()
returns the debug writer
Definition: debug_writer.h:391
TAlgebra::matrix_type matrix_type
type of matrix
Definition: debug_writer.h:363
Definition: densematrix.h:57
void subassign(size_t r, size_t c, const TStorage2 &t)
Definition: densematrix.h:215
Definition: densevector.h:101
Diagvanka Preconditioner, description see above diagvanka_step function.
Definition: vanka.h:308
number m_relax
Definition: vanka.h:351
virtual ~DiagVanka()
Destructor.
Definition: vanka.h:346
IPreconditioner< TAlgebra > base_type
Base type.
Definition: vanka.h:323
TAlgebra::matrix_type matrix_type
Matrix type.
Definition: vanka.h:317
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
Preprocess routine.
Definition: vanka.h:361
DiagVanka()
default constructor
Definition: vanka.h:332
virtual bool step(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp, vector_type &c, const vector_type &d)
computes a new correction c = B*d
Definition: vanka.h:377
TAlgebra algebra_type
Algebra type.
Definition: vanka.h:311
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition: vanka.h:320
virtual SmartPtr< ILinearIterator< vector_type > > clone()
Clone.
Definition: vanka.h:335
void set_relax(number omega)
set relaxation parameter
Definition: vanka.h:351
matrix_type m_A
Definition: vanka.h:411
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition: vanka.h:343
TAlgebra::vector_type vector_type
Vector type.
Definition: vanka.h:314
virtual const char * name() const
Name of preconditioner.
Definition: vanka.h:358
virtual bool postprocess()
Postprocess routine.
Definition: vanka.h:407
SmartPtr< IDamping< TAlgebra::vector_type, TAlgebra::vector_type > > damping()
returns the scaling
Definition: linear_iterator.h:173
describes a linear iterator that is based on a matrix operator
Definition: preconditioner.h:103
Definition: matrix_operator.h:49
Vanka Preconditioner.
Definition: vanka.h:198
Vanka()
default constructor
Definition: vanka.h:222
TAlgebra::matrix_type matrix_type
Matrix type.
Definition: vanka.h:207
TAlgebra algebra_type
Algebra type.
Definition: vanka.h:201
matrix_type m_A
Definition: vanka.h:300
virtual ~Vanka()
Destructor.
Definition: vanka.h:234
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition: vanka.h:237
virtual const char * name() const
Name of preconditioner.
Definition: vanka.h:248
TAlgebra::vector_type vector_type
Vector type.
Definition: vanka.h:204
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition: vanka.h:210
void set_relax(number omega)
set relaxation parameter
Definition: vanka.h:241
virtual SmartPtr< ILinearIterator< vector_type > > clone()
Clone.
Definition: vanka.h:225
number m_relax
Definition: vanka.h:241
IPreconditioner< TAlgebra > base_type
Base type.
Definition: vanka.h:213
virtual bool step(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp, vector_type &c, const vector_type &d)
computes a new correction c = B*d
Definition: vanka.h:267
virtual bool postprocess()
Postprocess routine.
Definition: vanka.h:296
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
Preprocess routine.
Definition: vanka.h:251
@ PST_UNIQUE
Definition: parallel_storage_type.h:70
void SetDirichletRow(TSparseMatrix &A, size_t i, size_t alpha)
Definition: sparsematrix_util.h:796
int NumProcs()
returns the number of processes
Definition: pcl_base.cpp:91
void CollectUniqueElements(std::vector< typename TLayout::Element > &elemsOut, const TLayout &layout)
writes all elements in the interfaces into the resulting vector. avoids doubles.
Definition: pcl_layout_util.h:142
number omega
#define UG_THROW(msg)
Definition: error.h:57
double number
Definition: types.h:124
the ug namespace
static const size_t MAXBLOCKSIZE
Definition: vanka.h:47
bool MakeConsistent(const ParallelMatrix< matrix_type > &_mat, ParallelMatrix< matrix_type > &newMat)
Definition: parallel_matrix_overlap_impl.h:438
bool MatMultAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition: operations_mat.h:68
bool Diag_Vanka_step(const Matrix_type &A, Vector_type &x, const Vector_type &b, number relax)
Definition: vanka.h:124
bool Vanka_step(const Matrix_type &A, Vector_type &x, const Vector_type &b, number relax)
Definition: vanka.h:50
bool InverseMatMult(number &dest, const double &beta, const TMat &mat, const TVec &vec)
you can implement this function with GetInverse and MatMult
T value_type
Definition: sparsematrix_interface.h:2