Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
element_gauss_seidel.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 * Authors: Christian Wehner, Andreas Vogel
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/*
34 * Remark: Base on Vanka idea and implementation
35 */
36
37#ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ELEMENT_GAUSS_SEIDEL__
38#define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ELEMENT_GAUSS_SEIDEL__
39
41
42#include <vector>
43#include <algorithm>
44
45#ifdef UG_PARALLEL
46 #include "pcl/pcl_util.h"
49#endif
50
51
52#define SCHUR_MOD
53namespace ug{
54
55template<typename TGroupObj, typename TDomain, typename TAlgebra>
56void ElementGaussSeidelStep(const typename TAlgebra::matrix_type& A,
58 const typename TAlgebra::vector_type& d,
59 number relax
60#ifdef SCHUR_MOD
61 , const std::vector<std::string>& cmp, number alpha, bool elim_off_diag=false
62#endif
63 )
64{
65 typedef typename TAlgebra::matrix_type::const_row_iterator const_row_iterator;
66 const static int blockSize = TAlgebra::blockSize;
67
68 // memory for local algebra
72 std::vector<size_t> vInd;
73
74 // set all vector entries to zero
75 c.set(0.0);
76#ifdef UG_PARALLEL
77 c.set_storage_type(PST_ADDITIVE);
78#endif
79 typedef typename GridFunction<TDomain, TAlgebra>::element_type Element;
80 std::vector<Element*> vElem;
81
82 // loop all grouping objects
83 typedef typename GridFunction<TDomain, TAlgebra>::template traits<TGroupObj>::const_iterator GroupObjIter;
84 for(GroupObjIter iter = c.template begin<TGroupObj>();
85 iter != c.template end<TGroupObj>(); ++iter){
86
87 // get grouping obj
88 TGroupObj* groupObj = *iter;
89
90 // collect elems associated to grouping object
91 c.collect_associated(vElem, groupObj);
92
93 // get all algebraic indices on element
94 vInd.clear();
95 for(size_t i = 0; i < vElem.size(); ++i)
96 c.algebra_indices(vElem[i], vInd, false);
97
98 // check for doublicates
99 if(vElem.size() > 1){
100 std::sort(vInd.begin(), vInd.end());
101 vInd.erase(std::unique(vInd.begin(), vInd.end()), vInd.end());
102 }
103
104 // get number of indices on patch
105 const size_t numIndex = vInd.size();
106
107 // fill local block matrix
108 bool bFound;
109 mat.resize(numIndex, numIndex);
110 mat = 0.0;
111 for (size_t j = 0; j<numIndex; j++){
112 for (size_t k=0;k<numIndex;k++){
113 const_row_iterator it = A.get_connection(vInd[j],vInd[k], bFound);
114 if(bFound){
115 mat.subassign(j*blockSize,k*blockSize,it.value());
116 }
117 };
118 }
119
120 // compute s[j] := d[j] - sum_k A(j,k)*c[k]
121 // note: the loop over k is the whole matrix row (not only selected indices)
122 s.resize(numIndex);
123 for (size_t j = 0; j<numIndex; j++)
124 {
125 typename TAlgebra::vector_type::value_type sj = d[vInd[j]];
126 for(const_row_iterator it = A.begin_row(vInd[j]); it != A.end_row(vInd[j]) ; ++it)
127 {
128 MatMultAdd(sj, 1.0, sj, -1.0, it.value(), c[it.index()]);
129 };
130
131 s.subassign(j*blockSize,sj);
132 };
133
134
135#ifdef SCHUR_MOD
136
137 // tokenize passed functions
139 c.approx_space()->dof_distribution_info();
140
141 // get fct IDs of selected functions
142 std::vector<size_t> vFct;
143 for(size_t i = 0; i < cmp.size(); ++i)
144 vFct.push_back(ddinfo->fct_id_by_name(cmp[i].c_str()));
145
146
148 weights.resize(numIndex);
149
150 if (vFct.size()==1)
151 {
152 std::vector<DoFIndex> vSchurDoFIndex;
153 c.dof_indices(groupObj, vFct[0], vSchurDoFIndex, true);
154
155 // identify schur rows
156 UG_ASSERT(blockSize==1, "Element GS: Elimination does only work for blockSize==1!")
157
158 std::vector<size_t> vIndElim; // global indices for elimination
159 // std::vector<size_t> vIndKeep; // global indices kept
160
161 const size_t numElim = vSchurDoFIndex.size();
162 for (size_t j = 0; j<numElim; j++)
163 {
164 vIndElim.push_back(vSchurDoFIndex[j][0]);
165 }
166
167 std::vector<size_t> mapElim; // local indices for elimination (j_elim -> j_local)
168 std::vector<size_t> mapKeep; // local indices for elimination (j_nonelim -> j_local)
169
170 // compute weights & fill map
171 for (size_t j = 0; j<numIndex; j++)
172 {
173 std::vector<size_t>::iterator it = find (vIndElim.begin(), vIndElim.end(), vInd[j]);
174 if (it != vIndElim.end()) {
175 // eliminate this row
176 mapElim.push_back(j);
177 } else {
178 // keep this row
179 mapKeep.push_back(j); // vIndKeep.push_back(vInd[j]);
180 }
181 }
182
183 const size_t numKeep = mapKeep.size();// vIndKeep.size();
184 UG_ASSERT((numElim+numKeep == numIndex), "Map elim does not match!");
185 UG_ASSERT(mapElim.size()==numElim, "Map elim does not match!");
186 UG_ASSERT(mapKeep.size()==numKeep, "Map mon elim does not match!");
187
188
189 // extract mat_keep (needed for inversion)
191 matKeep.resize(numKeep, numKeep);
192 matKeep = 0.0;
193 for (size_t j = 0; j<numKeep; j++)
194 {
195 for (size_t k=0;k<numKeep;k++)
196 {
197 matKeep(j,k) = mat(mapKeep[j], mapKeep[k]);
198 }
199 }
200
201 if (false) {
202
203 // compute mat elim
204 // compute contribution Bi^T A^-1 Cj to schur complement
205 // std::cout << "S (" << numElim << "):" << std::endl;
206
207 // compute
208 DenseVector< VariableArray1<number> > schur_cj; schur_cj.resize(numKeep);
209 DenseVector< VariableArray1<number> > schur_yj; schur_yj.resize(numKeep);
210
211 DenseMatrix< VariableArray2<number> > matElim; matElim.resize(numElim, numElim);
212 matElim = 0.0;
213
214
215 for (size_t j = 0; j<numElim; j++)
216 {
217 for (size_t k=0; k<numKeep; k++)
218 {
219 schur_cj[k] = mat(mapKeep[k], mapElim[j]);
220 }
221
222 InverseMatMult(schur_yj, 1.0, matKeep, schur_cj);
223
224 for (size_t i = 0; i<numElim; i++)
225 {
226 // compute elim_ij
227 matElim(i,j) = 0.0;
228 for (size_t k=0; k<numKeep; k++)
229 {
230 number schur_bik = mat(mapElim[i], mapKeep[k]);
231 matElim(i,j) += schur_bik*schur_yj[k];
232
233 }
234
235 // replace/update matrix
236 // std::cout << mat(mapElim[i], mapElim[j]) << "+" << matElim(i,j) << "=";
237 mat(mapElim[i], mapElim[j]) -= alpha*matElim(i,j); //* alpha
238 // std::cout << mat(mapElim[i], mapElim[j]) << std::endl;
239
240 // std::cout << matElim;
241 }
242
243
244
245 }
246 }
247
248
249 if (false) // eliminate off-diag rows B, C?
250 {
251 for (size_t j = 0; j<numElim; j++)
252 {
253 for (size_t k=0; k<numKeep; k++)
254 {
255 mat(mapElim[j], mapKeep[k]) = 0.0;
256 mat(mapKeep[k], mapElim[j]) = 0.0;
257 }
258
259 for (size_t k=0; k<numElim; k++)
260 {
261
262 if (j==k) continue;
263 mat(mapElim[j], mapElim[k]) = 0.0;
264 mat(mapElim[k], mapElim[j]) = 0.0;
265 }
266 }
267 }
268
269
270
271 if (false) // rescale elim part of matrix?
272 {
273 for (size_t i = 0; i<numElim; i++)
274 {
275 for (size_t j=0; j<numElim; j++)
276 {
277 mat(mapElim[i], mapElim[j]) *= alpha;
278
279 }
280 }
281 }
282
283 // std::cout << mat;
284
285 // std::cout << std::endl;
286 } // (vFct.size()==1)
287
288#endif
289 // solve block
290 x.resize(numIndex);
291 InverseMatMult(x, 1.0, mat, s);
292 for (size_t j=0;j<numIndex;j++)
293 {
294 c[vInd[j]] += relax*x[j];
295 }
296
297
298 }
299}
300
301
303template <typename TDomain, typename TAlgebra>
304class ElementGaussSeidel : public IPreconditioner<TAlgebra>
305{
306 public:
308 typedef TAlgebra algebra_type;
309
311 typedef typename TAlgebra::vector_type vector_type;
312
314 typedef typename TAlgebra::matrix_type matrix_type;
315
318
321
322 protected:
326
328
329 public:
331 ElementGaussSeidel() : m_relax(1.0), m_type("element"), m_schur_alpha(1.0), m_elim_off_diag(false) {};
332
334 ElementGaussSeidel(number relax) : m_relax(relax), m_type("element"), m_schur_alpha(1.0), m_elim_off_diag(false) {};
335
337 ElementGaussSeidel(const std::string& type) : m_relax(1.0), m_type(type), m_schur_alpha(1.0), m_elim_off_diag(false) {};
338
340 ElementGaussSeidel(number relax, const std::string& type) : m_relax(relax), m_type(type), m_schur_alpha(1.0), m_elim_off_diag(false) {};
341
344 {
347 newInst->set_debug(debug_writer());
348 newInst->set_damp(this->damping());
349 newInst->set_relax(m_relax);
350 newInst->set_type(m_type);
351
352 newInst->m_schur_cmp =m_schur_cmp;
353 newInst->m_schur_alpha = m_schur_alpha;
354 return newInst;
355 }
356
359 {};
360
362 virtual bool supports_parallel() const {return true;}
363
365 void set_relax(number omega){m_relax=omega;};
366
368 void set_type(const std::string& type){m_type=type;};
369
370 void select_schur_cmp(const std::vector<std::string>& cmp, number alpha)
371 {
372 m_schur_cmp =cmp;
373 m_schur_alpha = alpha;
374 };
375
376 void set_elim_offdiag(bool elim) { m_elim_off_diag=elim;};
377
378 protected:
380 virtual const char* name() const {return "ElementGaussSeidel";}
381
384 {
385#ifdef UG_PARALLEL
386 if(pcl::NumProcs() > 1)
387 {
388 // copy original matrix
389 MakeConsistent(*pOp, m_A);
390 // set zero on slaves
391 // std::vector<IndexLayout::Element> vIndex;
392 //CollectUniqueElements(vIndex, m_A.layouts()->slave());
393 //SetDirichletRow(m_A, vIndex);
394 }
395#endif
396 return true;
397 }
398
400 {
402 = dynamic_cast<GridFunction<TDomain, TAlgebra>*>(&c);
403 if(pC == NULL)
404 UG_THROW("ElementGaussSeidel expects correction to be a GridFunction.");
405
406
407 typedef typename GridFunction<TDomain, TAlgebra>::element_type Element;
408 typedef typename GridFunction<TDomain, TAlgebra>::side_type Side;
409
410
411#ifdef UG_PARALLEL
412 if(pcl::NumProcs() > 1){
413 // make defect unique
414 SmartPtr<vector_type> spDtmp = d.clone();
415 spDtmp->change_storage_type(PST_UNIQUE);
416
417 // execute step
418 if (m_type == "element") ElementGaussSeidelStep<Element,TDomain,TAlgebra>(m_A, *pC, *spDtmp, m_relax, m_schur_cmp, m_schur_alpha);
419 else UG_THROW("ElementGaussSeidel: wrong patch type '"<<m_type<<"'."
420 " Options: element, side, face, edge, vertex.");
421
422 // make correction consistent
423 // pC->set_storage_type(PST_CONSISTENT);
424 pC->change_storage_type(PST_CONSISTENT);
425 return true;
426 }
427 else
428#endif
429 {
430 matrix_type &A=*pOp;
431 if (m_type == "element") ElementGaussSeidelStep<Element,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
432 else if (m_type == "side") ElementGaussSeidelStep<Side,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
433 else if (m_type == "face") ElementGaussSeidelStep<Face,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
434 else if (m_type == "edge") ElementGaussSeidelStep<Edge,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
435 else if (m_type == "vertex") ElementGaussSeidelStep<Vertex,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
436 else UG_THROW("ElementGaussSeidel: wrong patch type '"<<m_type<<"'."
437 " Options: element, side, face, edge, vertex.")
438#ifdef UG_PARALLEL
439 pC->set_storage_type(PST_CONSISTENT);
440#endif
441
442 return true;
443 }
444 }
445
447 virtual bool postprocess() {return true;}
448
449 protected:
450#ifdef UG_PARALLEL
452#endif
453
455 std::string m_type;
456
457#ifdef SCHUR_MOD
458 std::vector<std::string> m_schur_cmp;
461#endif
462
463};
464
465} // end namespace ug
466
467#endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ELEMENT_GAUSS_SEIDEL__ */
parameterString s
Definition smart_pointer.h:296
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
SmartPtr< IDebugWriter< algebra_type > > debug_writer()
returns the debug writer
Definition debug_writer.h:391
Definition densematrix.h:57
void subassign(size_t r, size_t c, const TStorage2 &t)
Definition densematrix.h:215
Definition densevector.h:101
ElementGaussSeidel Preconditioner.
Definition element_gauss_seidel.h:305
number m_schur_alpha
Definition element_gauss_seidel.h:459
bool m_elim_off_diag
Definition element_gauss_seidel.h:460
void select_schur_cmp(const std::vector< std::string > &cmp, number alpha)
Definition element_gauss_seidel.h:370
ElementGaussSeidel(number relax, const std::string &type)
constructor setting relaxation and type
Definition element_gauss_seidel.h:340
virtual ~ElementGaussSeidel()
Destructor.
Definition element_gauss_seidel.h:358
void set_type(const std::string &type)
set type
Definition element_gauss_seidel.h:368
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
Preprocess routine.
Definition element_gauss_seidel.h:383
virtual const char * name() const
Name of preconditioner.
Definition element_gauss_seidel.h:380
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition element_gauss_seidel.h:362
TAlgebra::vector_type vector_type
Vector type.
Definition element_gauss_seidel.h:311
number m_relax
Definition element_gauss_seidel.h:454
ElementGaussSeidel()
default constructor
Definition element_gauss_seidel.h:331
GridFunction< TDomain, TAlgebra > grid_function_type
Definition element_gauss_seidel.h:327
TAlgebra algebra_type
Algebra type.
Definition element_gauss_seidel.h:308
ElementGaussSeidel(const std::string &type)
constructor setting type
Definition element_gauss_seidel.h:337
matrix_type m_A
Definition element_gauss_seidel.h:451
virtual SmartPtr< ILinearIterator< vector_type > > clone()
Clone.
Definition element_gauss_seidel.h:343
std::vector< std::string > m_schur_cmp
Definition element_gauss_seidel.h:458
void set_relax(number omega)
set relaxation parameter
Definition element_gauss_seidel.h:365
virtual bool postprocess()
Postprocess routine.
Definition element_gauss_seidel.h:447
std::string m_type
Definition element_gauss_seidel.h:455
void set_elim_offdiag(bool elim)
Definition element_gauss_seidel.h:376
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition element_gauss_seidel.h:317
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 element_gauss_seidel.h:399
ElementGaussSeidel(number relax)
constructor setting relaxation
Definition element_gauss_seidel.h:334
TAlgebra::matrix_type matrix_type
Matrix type.
Definition element_gauss_seidel.h:314
IPreconditioner< TAlgebra > base_type
Base type.
Definition element_gauss_seidel.h:320
represents numerical solutions on a grid using an algebraic vector
Definition grid_function.h:121
size_t dof_indices(TElem *elem, size_t fct, std::vector< DoFIndex > &ind, bool bHang=false, bool bClear=true) const
get multi indices on an finite element in canonical order
Definition grid_function.h:322
size_t algebra_indices(TElem *elem, std::vector< size_t > &ind, bool bClear=true) const
get algebra indices on an geometric object in canonical order
Definition grid_function.h:332
void collect_associated(std::vector< TBaseElem * > &vAssElem, TElem *elem, bool clearContainer=true) const
returns the adjacend elements
Definition grid_function.h:290
SmartPtr< ApproximationSpace< TDomain > > approx_space()
returns approx space
Definition grid_function.h:348
dim_traits< dim >::grid_base_object element_type
Definition grid_function.h:159
dim_traits< dim-1 >::grid_base_object side_type
Definition grid_function.h:163
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
#define SCHUR_MOD
Definition element_gauss_seidel.h:52
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
@ PST_UNIQUE
Definition parallel_storage_type.h:70
@ PST_ADDITIVE
Definition parallel_storage_type.h:69
int NumProcs()
returns the number of processes
Definition pcl_base.cpp:91
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_THROW(msg)
Definition error.h:57
double number
Definition types.h:124
the ug namespace
IndexLayout::Interface::iterator find(IndexLayout::Interface &interface, size_t i)
Definition parallel_index_layout.h:77
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
void ElementGaussSeidelStep(const typename TAlgebra::matrix_type &A, GridFunction< TDomain, TAlgebra > &c, const typename TAlgebra::vector_type &d, number relax, const std::vector< std::string > &cmp, number alpha, bool elim_off_diag=false)
Definition element_gauss_seidel.h:56
bool InverseMatMult(number &dest, const double &beta, const TMat &mat, const TVec &vec)
you can implement this function with GetInverse and MatMult