Line data Source code
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 :
40 : #include "lib_algebra/operator/interface/preconditioner.h"
41 :
42 : #include <vector>
43 : #include <algorithm>
44 :
45 : #ifdef UG_PARALLEL
46 : #include "pcl/pcl_util.h"
47 : #include "lib_algebra/parallelization/parallelization_util.h"
48 : #include "lib_algebra/parallelization/parallel_matrix_overlap_impl.h"
49 : #endif
50 :
51 :
52 : #define SCHUR_MOD
53 : namespace ug{
54 :
55 : template<typename TGroupObj, typename TDomain, typename TAlgebra>
56 0 : void ElementGaussSeidelStep(const typename TAlgebra::matrix_type& A,
57 : GridFunction<TDomain, TAlgebra>& c,
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
69 : DenseVector< VariableArray1<number> > s;
70 : DenseVector< VariableArray1<number> > x;
71 : DenseMatrix< VariableArray2<number> > mat;
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 0 : for(GroupObjIter iter = c.template begin<TGroupObj>();
85 0 : 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 0 : for(size_t i = 0; i < vElem.size(); ++i)
96 0 : c.algebra_indices(vElem[i], vInd, false);
97 :
98 : // check for doublicates
99 0 : if(vElem.size() > 1){
100 0 : std::sort(vInd.begin(), vInd.end());
101 0 : 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 0 : mat.resize(numIndex, numIndex);
110 0 : mat = 0.0;
111 0 : for (size_t j = 0; j<numIndex; j++){
112 0 : for (size_t k=0;k<numIndex;k++){
113 0 : const_row_iterator it = A.get_connection(vInd[j],vInd[k], bFound);
114 0 : if(bFound){
115 0 : 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 0 : s.resize(numIndex);
123 0 : for (size_t j = 0; j<numIndex; j++)
124 : {
125 0 : typename TAlgebra::vector_type::value_type sj = d[vInd[j]];
126 0 : for(const_row_iterator it = A.begin_row(vInd[j]); it != A.end_row(vInd[j]) ; ++it)
127 : {
128 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), c[it.index()]);
129 : };
130 :
131 0 : s.subassign(j*blockSize,sj);
132 : };
133 :
134 :
135 : #ifdef SCHUR_MOD
136 :
137 : // tokenize passed functions
138 : ConstSmartPtr<DoFDistributionInfo> ddinfo =
139 0 : c.approx_space()->dof_distribution_info();
140 :
141 : // get fct IDs of selected functions
142 : std::vector<size_t> vFct;
143 0 : for(size_t i = 0; i < cmp.size(); ++i)
144 0 : vFct.push_back(ddinfo->fct_id_by_name(cmp[i].c_str()));
145 :
146 :
147 : DenseVector< VariableArray1<number> > weights;
148 0 : weights.resize(numIndex);
149 :
150 0 : if (vFct.size()==1)
151 : {
152 : std::vector<DoFIndex> vSchurDoFIndex;
153 0 : 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 0 : for (size_t j = 0; j<numElim; j++)
163 : {
164 0 : 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 0 : for (size_t j = 0; j<numIndex; j++)
172 : {
173 0 : std::vector<size_t>::iterator it = find (vIndElim.begin(), vIndElim.end(), vInd[j]);
174 0 : if (it != vIndElim.end()) {
175 : // eliminate this row
176 0 : mapElim.push_back(j);
177 : } else {
178 : // keep this row
179 0 : 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)
190 : DenseMatrix< VariableArray2<number> > matKeep;
191 0 : matKeep.resize(numKeep, numKeep);
192 0 : matKeep = 0.0;
193 0 : for (size_t j = 0; j<numKeep; j++)
194 : {
195 0 : for (size_t k=0;k<numKeep;k++)
196 : {
197 0 : 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 0 : } // (vFct.size()==1)
287 :
288 : #endif
289 : // solve block
290 0 : x.resize(numIndex);
291 0 : InverseMatMult(x, 1.0, mat, s);
292 0 : for (size_t j=0;j<numIndex;j++)
293 : {
294 0 : c[vInd[j]] += relax*x[j];
295 : }
296 :
297 :
298 : }
299 0 : }
300 :
301 :
302 : /// ElementGaussSeidel Preconditioner
303 : template <typename TDomain, typename TAlgebra>
304 : class ElementGaussSeidel : public IPreconditioner<TAlgebra>
305 : {
306 : public:
307 : /// Algebra type
308 : typedef TAlgebra algebra_type;
309 :
310 : /// Vector type
311 : typedef typename TAlgebra::vector_type vector_type;
312 :
313 : /// Matrix type
314 : typedef typename TAlgebra::matrix_type matrix_type;
315 :
316 : /// Matrix Operator type
317 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
318 :
319 : /// Base type
320 : typedef IPreconditioner<TAlgebra> base_type;
321 :
322 : protected:
323 : using base_type::set_debug;
324 : using base_type::debug_writer;
325 : using base_type::write_debug;
326 :
327 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
328 :
329 : public:
330 : /// default constructor
331 0 : ElementGaussSeidel() : m_relax(1.0), m_type("element"), m_schur_alpha(1.0), m_elim_off_diag(false) {};
332 :
333 : /// constructor setting relaxation
334 0 : ElementGaussSeidel(number relax) : m_relax(relax), m_type("element"), m_schur_alpha(1.0), m_elim_off_diag(false) {};
335 :
336 : /// constructor setting type
337 0 : ElementGaussSeidel(const std::string& type) : m_relax(1.0), m_type(type), m_schur_alpha(1.0), m_elim_off_diag(false) {};
338 :
339 : /// constructor setting relaxation and type
340 0 : ElementGaussSeidel(number relax, const std::string& type) : m_relax(relax), m_type(type), m_schur_alpha(1.0), m_elim_off_diag(false) {};
341 :
342 : /// Clone
343 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
344 : {
345 : SmartPtr<ElementGaussSeidel<TDomain, TAlgebra> >
346 0 : newInst(new ElementGaussSeidel<TDomain, TAlgebra>());
347 0 : newInst->set_debug(debug_writer());
348 0 : newInst->set_damp(this->damping());
349 0 : newInst->set_relax(m_relax);
350 0 : newInst->set_type(m_type);
351 :
352 0 : newInst->m_schur_cmp =m_schur_cmp;
353 0 : newInst->m_schur_alpha = m_schur_alpha;
354 0 : return newInst;
355 : }
356 :
357 : /// Destructor
358 0 : virtual ~ElementGaussSeidel()
359 0 : {};
360 :
361 : /// returns if parallel solving is supported
362 0 : virtual bool supports_parallel() const {return true;}
363 :
364 : /// set relaxation parameter
365 0 : void set_relax(number omega){m_relax=omega;};
366 :
367 : /// set type
368 0 : void set_type(const std::string& type){m_type=type;};
369 :
370 0 : void select_schur_cmp(const std::vector<std::string>& cmp, number alpha)
371 : {
372 0 : m_schur_cmp =cmp;
373 0 : m_schur_alpha = alpha;
374 0 : };
375 :
376 0 : void set_elim_offdiag(bool elim) { m_elim_off_diag=elim;};
377 :
378 : protected:
379 : /// Name of preconditioner
380 0 : virtual const char* name() const {return "ElementGaussSeidel";}
381 :
382 : /// Preprocess routine
383 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
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 0 : return true;
397 : }
398 :
399 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
400 : {
401 : GridFunction<TDomain, TAlgebra>* pC
402 0 : = dynamic_cast<GridFunction<TDomain, TAlgebra>*>(&c);
403 0 : if(pC == NULL)
404 0 : 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 0 : matrix_type &A=*pOp;
431 0 : if (m_type == "element") ElementGaussSeidelStep<Element,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
432 0 : else if (m_type == "side") ElementGaussSeidelStep<Side,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
433 0 : else if (m_type == "face") ElementGaussSeidelStep<Face,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
434 0 : else if (m_type == "edge") ElementGaussSeidelStep<Edge,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
435 0 : else if (m_type == "vertex") ElementGaussSeidelStep<Vertex,TDomain,TAlgebra>(A, *pC, d, m_relax, m_schur_cmp, m_schur_alpha);
436 0 : 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 0 : return true;
443 : }
444 : }
445 :
446 : /// Postprocess routine
447 0 : virtual bool postprocess() {return true;}
448 :
449 : protected:
450 : #ifdef UG_PARALLEL
451 : matrix_type m_A;
452 : #endif
453 :
454 : number m_relax;
455 : std::string m_type;
456 :
457 : #ifdef SCHUR_MOD
458 : std::vector<std::string> m_schur_cmp;
459 : number m_schur_alpha;
460 : bool m_elim_off_diag;
461 : #endif
462 :
463 : };
464 :
465 : } // end namespace ug
466 :
467 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ELEMENT_GAUSS_SEIDEL__ */
|
