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__COMPONENT_GAUSS_SEIDEL__
38 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__COMPONENT_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 : namespace ug{
52 :
53 : /// ComponentGaussSeidel Preconditioner
54 : template <typename TDomain, typename TAlgebra>
55 : class ComponentGaussSeidel : public IPreconditioner<TAlgebra>
56 : {
57 : public:
58 : /// World dimension
59 : typedef GridFunction<TDomain, TAlgebra> GF;
60 : typedef typename GF::element_type Element;
61 : typedef typename GF::side_type Side;
62 : static const int dim = TDomain::dim;
63 :
64 : /// Algebra types
65 : /// \{
66 : typedef TAlgebra algebra_type;
67 : typedef typename TAlgebra::vector_type vector_type;
68 : typedef typename TAlgebra::matrix_type matrix_type;
69 : /// \}
70 :
71 : protected:
72 : typedef IPreconditioner<TAlgebra> base_type;
73 : using base_type::set_debug;
74 : using base_type::debug_writer;
75 : using base_type::write_debug;
76 :
77 : public:
78 : /// default constructor
79 0 : ComponentGaussSeidel(const std::vector<std::string>& vFullRowCmp)
80 0 : : m_bInit(false), m_relax(1), m_alpha(1.0), m_beta(1.0), m_bWeighted(false), m_vFullRowCmp(vFullRowCmp) {}
81 :
82 : /// constructor setting relaxation and type
83 0 : ComponentGaussSeidel(number relax, const std::vector<std::string>& vFullRowCmp)
84 0 : : m_bInit(false), m_relax(relax), m_alpha(1.0), m_beta(1.0), m_bWeighted(false), m_vFullRowCmp(vFullRowCmp){}
85 :
86 : /// constructor setting relaxation and type
87 0 : ComponentGaussSeidel(number relax, const std::vector<std::string>& vFullRowCmp,
88 : const std::vector<int>& vSmooth, const std::vector<number>& vDamp)
89 0 : : m_bInit(false), m_relax(relax), m_alpha(1.0), m_beta(1.0), m_bWeighted(false), m_vFullRowCmp(vFullRowCmp),
90 0 : m_vGroupObj(vSmooth), m_vDamp(vDamp) {};
91 :
92 : /// Clone
93 : virtual SmartPtr<ILinearIterator<vector_type> > clone();
94 :
95 : /// returns if parallel solving is supported
96 0 : virtual bool supports_parallel() const {return true;}
97 :
98 0 : void set_alpha(number alpha)
99 0 : {m_alpha = alpha;}
100 :
101 0 : void set_beta(number beta)
102 0 : {m_beta = beta;}
103 :
104 0 : void set_weights(bool b)
105 0 : {m_bWeighted = b;}
106 :
107 : protected:
108 : /// Name of preconditioner
109 0 : virtual const char* name() const {return "ComponentGaussSeidel";}
110 :
111 : /// Preprocess routine
112 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp);
113 :
114 : struct DimCache{
115 : std::vector<std::vector<DoFIndex> > vvDoFIndex;
116 : std::vector<DenseMatrix< VariableArray2<number> > > vBlockInverse;
117 : };
118 :
119 : void apply_blocks(const matrix_type& A, GF& c,
120 : const vector_type& d, number relax,
121 : const DimCache& dimCache,
122 : bool bReverse);
123 :
124 : void apply_blocks_weighted(const matrix_type& A, GF& c,
125 : const vector_type& d, number relax,
126 : const DimCache& dimCache,
127 : bool bReverse);
128 :
129 : void extract_blocks(const matrix_type& A, const GF& c);
130 :
131 : /// step method
132 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d);
133 :
134 : /// Postprocess routine
135 0 : virtual bool postprocess() {return true;}
136 :
137 : protected:
138 : #ifdef UG_PARALLEL
139 : matrix_type m_A;
140 : #endif
141 :
142 : /// init flag
143 : bool m_bInit;
144 :
145 : /// relaxing parameter
146 : number m_relax;
147 : number m_alpha;
148 : number m_beta;
149 :
150 : bool m_bWeighted;
151 : SmartPtr<GF> m_weight;
152 :
153 : /// components, whose matrix row must be fulfilled completely on gs-block
154 : std::vector<std::string> m_vFullRowCmp;
155 :
156 : /// smooth order
157 : std::vector<int> m_vGroupObj;
158 : std::vector<number> m_vDamp;
159 :
160 : /// caching storage
161 : DimCache m_vDimCache[VOLUME+1];
162 : };
163 :
164 : template <typename TDomain, typename TAlgebra>
165 0 : bool ComponentGaussSeidel<TDomain, TAlgebra>::
166 : preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
167 : {
168 : #ifdef UG_PARALLEL
169 : if(pcl::NumProcs() > 1)
170 : {
171 : // copy original matrix
172 : MakeConsistent(*pOp, m_A);
173 : // set zero on slaves
174 : std::vector<IndexLayout::Element> vIndex;
175 : CollectUniqueElements(vIndex, m_A.layouts()->slave());
176 : SetDirichletRow(m_A, vIndex);
177 : }
178 : #endif
179 :
180 0 : m_bInit = false;
181 :
182 0 : return true;
183 : }
184 :
185 : template <typename TDomain, typename TAlgebra>
186 0 : void ComponentGaussSeidel<TDomain, TAlgebra>::
187 : apply_blocks(const matrix_type& A, GF& c,
188 : const vector_type& d, number relax,
189 : const DimCache& dimCache,
190 : bool bReverse)
191 : {
192 : // memory for local algebra
193 : DenseVector< VariableArray1<number> > s;
194 : DenseVector< VariableArray1<number> > x;
195 :
196 : // get caches
197 : const std::vector<std::vector<DoFIndex> >& vvDoFIndex = dimCache.vvDoFIndex;
198 : const std::vector<DenseMatrix< VariableArray2<number> > >& vBlockInverse = dimCache.vBlockInverse;
199 :
200 : // loop blocks
201 0 : for(size_t b = 0; b < vvDoFIndex.size(); ++b)
202 : {
203 0 : size_t block = (!bReverse) ? b : (vvDoFIndex.size()-1 - b);
204 :
205 : // get storage
206 : const std::vector<DoFIndex>& vDoFIndex = vvDoFIndex[block];
207 : const DenseMatrix<VariableArray2<number> >& BlockInv = vBlockInverse[block];
208 : const size_t numIndex = vDoFIndex.size();
209 :
210 : // compute s[j] := d[j] - sum_k A(j,k)*c[k]
211 : // note: the loop over k is the whole matrix row (not only selected indices)
212 : typedef typename matrix_type::const_row_iterator const_row_iterator;
213 : const static int blockSize = TAlgebra::blockSize;
214 0 : s.resize(numIndex);
215 0 : for (size_t j = 0; j<numIndex; j++){
216 0 : typename vector_type::value_type sj = d[vDoFIndex[j][0]];
217 : for(const_row_iterator it = A.begin_row(vDoFIndex[j][0]);
218 0 : it != A.end_row(vDoFIndex[j][0]) ; ++it){
219 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), c[it.index()]);
220 : };
221 0 : s.subassign(j*blockSize,sj);
222 : };
223 :
224 : // solve block
225 0 : x.resize(numIndex);
226 0 : MatMult(x, relax, BlockInv, s);
227 :
228 : // add to global correction
229 0 : for (size_t j=0;j<numIndex;j++)
230 0 : DoFRef(c, vDoFIndex[j]) += x[j];
231 : }
232 0 : }
233 :
234 :
235 : template <typename TDomain, typename TAlgebra>
236 0 : void ComponentGaussSeidel<TDomain, TAlgebra>::
237 : apply_blocks_weighted(const matrix_type& A, GF& c,
238 : const vector_type& d, number relax,
239 : const DimCache& dimCache,
240 : bool bReverse)
241 : {
242 : // memory for local algebra
243 : DenseVector< VariableArray1<number> > s;
244 : DenseVector< VariableArray1<number> > x;
245 :
246 : // get caches
247 : const std::vector<std::vector<DoFIndex> >& vvDoFIndex = dimCache.vvDoFIndex;
248 : const std::vector<DenseMatrix< VariableArray2<number> > >& vBlockInverse = dimCache.vBlockInverse;
249 :
250 : // loop blocks
251 0 : for(size_t b = 0; b < vvDoFIndex.size(); ++b)
252 : {
253 0 : size_t block = (!bReverse) ? b : (vvDoFIndex.size()-1 - b);
254 :
255 : // get storage
256 : const std::vector<DoFIndex>& vDoFIndex = vvDoFIndex[block];
257 : const DenseMatrix<VariableArray2<number> >& BlockInv = vBlockInverse[block];
258 : const size_t numIndex = vDoFIndex.size();
259 :
260 : // compute s[i] := d[i] - sum_k A(i,k)*c[k]
261 : // note: the loop over k is the whole matrix row (not only selected indices)
262 : typedef typename matrix_type::const_row_iterator const_row_iterator;
263 : const static int blockSize = TAlgebra::blockSize;
264 0 : s.resize(numIndex);
265 0 : for (size_t i = 0; i<numIndex; i++)
266 : {
267 0 : typename vector_type::value_type si = d[vDoFIndex[i][0]];
268 : for(const_row_iterator it = A.begin_row(vDoFIndex[i][0]);
269 0 : it != A.end_row(vDoFIndex[i][0]) ; ++it)
270 : {
271 0 : MatMultAdd(si, 1.0, si, -1.0, it.value(), c[it.index()]);
272 : }
273 0 : s.subassign(i*blockSize, si);
274 : }
275 :
276 : // restrict
277 0 : for (size_t i = 0; i<numIndex; ++i)
278 : {
279 0 : number wi = sqrt(DoFRef(*m_weight, vDoFIndex[i]));
280 : UG_ASSERT(wi!=0.0, "Huhh: Zero weights found!" << i << "dofIndex" << vDoFIndex[i])
281 0 : s[i] /= wi;
282 : }
283 :
284 : // solve block
285 0 : x.resize(numIndex);
286 0 : MatMult(x, relax, BlockInv, s);
287 :
288 : // add to global correction
289 0 : for (size_t i=0; i<numIndex;i++)
290 : {
291 0 : number wi = sqrt(DoFRef(*m_weight, vDoFIndex[i]));
292 0 : DoFRef(c, vDoFIndex[i]) += x[i]/wi;
293 : }
294 : }
295 0 : }
296 :
297 : // extract_by_grouping(std::vector<std::vector<DoFIndex> >& vvDoFIndex,
298 : template<typename TGroupObj, typename GF>
299 0 : void extract_by_grouping(std::vector<std::vector<DoFIndex> >& vvDoFIndex,
300 : const GF& c,
301 : const std::vector<size_t>& vFullRowCmp,
302 : const std::vector<size_t>& vRemainCmp)
303 : {
304 :
305 : typedef typename GF::element_type Element;
306 :
307 : // memory for local algebra
308 : std::vector<DoFIndex> vFullRowDoFIndex;
309 : std::vector<Element*> vElem;
310 :
311 : // clear indices
312 : // \todo: improve this by precomputing size
313 : vvDoFIndex.clear();
314 :
315 : // loop all grouping objects
316 : typedef typename GF::template traits<TGroupObj>::const_iterator GroupObjIter;
317 0 : for(GroupObjIter iter = c.template begin<TGroupObj>();
318 0 : iter != c.template end<TGroupObj>(); ++iter)
319 : {
320 : // get grouping obj
321 : TGroupObj* groupObj = *iter;
322 :
323 : // get all dof indices on obj associated to cmps
324 : vFullRowDoFIndex.clear();
325 0 : for(size_t f = 0; f < vFullRowCmp.size(); ++f)
326 0 : c.inner_dof_indices(groupObj, vFullRowCmp[f], vFullRowDoFIndex, false);
327 :
328 : // check if equation present
329 0 : if(vFullRowDoFIndex.empty())
330 0 : UG_THROW("extract_by_grouping: Should not happen.");
331 :
332 : // get all algebraic indices on element
333 : if(TGroupObj::dim <= VERTEX){
334 : std::vector<Edge*> vSub;
335 : c.collect_associated(vSub, groupObj);
336 0 : for(size_t i = 0; i < vSub.size(); ++i)
337 0 : for(size_t f = 0; f < vFullRowCmp.size(); ++f)
338 0 : c.inner_dof_indices(vSub[i], vFullRowCmp[f], vFullRowDoFIndex, false);
339 0 : }
340 : if(TGroupObj::dim <= EDGE){
341 : std::vector<Face*> vSub;
342 : c.collect_associated(vSub, groupObj);
343 0 : for(size_t i = 0; i < vSub.size(); ++i)
344 0 : for(size_t f = 0; f < vFullRowCmp.size(); ++f)
345 0 : c.inner_dof_indices(vSub[i], vFullRowCmp[f], vFullRowDoFIndex, false);
346 0 : }
347 :
348 : // collect elems associated to grouping object
349 : c.collect_associated(vElem, groupObj);
350 :
351 : // get all algebraic indices on element
352 : std::vector<DoFIndex> vDoFIndex;
353 0 : for(size_t i = 0; i < vElem.size(); ++i)
354 0 : for(size_t f = 0; f < vRemainCmp.size(); ++f)
355 0 : c.dof_indices(vElem[i], vRemainCmp[f], vDoFIndex, false, false);
356 :
357 : // check for dublicates
358 0 : if(vElem.size() > 1){
359 0 : std::sort(vDoFIndex.begin(), vDoFIndex.end());
360 0 : vDoFIndex.erase(std::unique(vDoFIndex.begin(), vDoFIndex.end()), vDoFIndex.end());
361 : }
362 :
363 : // concat all indices
364 0 : vDoFIndex.insert(vDoFIndex.end(), vFullRowDoFIndex.begin(), vFullRowDoFIndex.end());
365 :
366 : // add
367 0 : vvDoFIndex.push_back(vDoFIndex);
368 : }
369 0 : }
370 :
371 : template <typename TDomain, typename TAlgebra>
372 0 : void ComponentGaussSeidel<TDomain, TAlgebra>::
373 : extract_blocks(const matrix_type& A, const GF& c)
374 : {
375 :
376 : ConstSmartPtr<DoFDistributionInfo> ddinfo =
377 0 : c.approx_space()->dof_distribution_info();
378 :
379 : // tokenize passed functions
380 0 : if(m_vFullRowCmp.empty())
381 0 : UG_THROW("ComponentGaussSeidel: No components set.")
382 :
383 : // get ids of selected functions
384 : std::vector<size_t> vFullRowCmp;
385 0 : for(size_t i = 0; i < m_vFullRowCmp.size(); ++i)
386 0 : vFullRowCmp.push_back(ddinfo->fct_id_by_name(m_vFullRowCmp[i].c_str()));
387 :
388 : // compute remaining functions
389 : std::vector<size_t> vRemainCmp;
390 0 : for(size_t f = 0; f < ddinfo->num_fct(); ++f)
391 0 : if(std::find(vFullRowCmp.begin(), vFullRowCmp.end(), f) == vFullRowCmp.end())
392 0 : vRemainCmp.push_back(f);
393 :
394 : /// Create 'm_vGroupObj' contain all dimensions containg objects from 'vFullRowCmp'
395 0 : if(m_vGroupObj.empty()){
396 0 : for(int d = VERTEX; d <= VOLUME; ++d){
397 : bool bCarryDoFs = false;
398 0 : for(size_t f = 0; f < vFullRowCmp.size(); ++f)
399 0 : if(ddinfo->max_fct_dofs(vFullRowCmp[f], d) > 0)
400 : bCarryDoFs = true;
401 0 : if(bCarryDoFs)
402 0 : m_vGroupObj.push_back(d);
403 : }
404 : }
405 :
406 : // extract for each dim-grouping objects
407 0 : for(int d = VERTEX; d <= VOLUME; ++d)
408 : {
409 : // only extract if needed
410 0 : if(std::find(m_vGroupObj.begin(), m_vGroupObj.end(), d) == m_vGroupObj.end())
411 0 : continue;
412 :
413 : // only extract if carrying dofs
414 : bool bCarryDoFs = false;
415 0 : for(size_t f = 0; f < vFullRowCmp.size(); ++f)
416 0 : if(ddinfo->max_fct_dofs(vFullRowCmp[f], d) > 0)
417 : bCarryDoFs = true;
418 0 : if(!bCarryDoFs)
419 0 : continue;
420 :
421 : // extract
422 0 : std::vector<std::vector<DoFIndex> >& vvDoFIndex = m_vDimCache[d].vvDoFIndex;
423 0 : switch(d){
424 0 : case VERTEX: extract_by_grouping<Vertex, GF>(vvDoFIndex, c, vFullRowCmp, vRemainCmp); break;
425 0 : case EDGE: extract_by_grouping<Edge, GF>(vvDoFIndex, c, vFullRowCmp, vRemainCmp); break;
426 0 : case FACE: extract_by_grouping<Face, GF>(vvDoFIndex, c, vFullRowCmp, vRemainCmp); break;
427 0 : case VOLUME: extract_by_grouping<Volume, GF>(vvDoFIndex, c, vFullRowCmp, vRemainCmp); break;
428 0 : default: UG_THROW("wrong dim");
429 : }
430 : }
431 :
432 : // compute weights v(i) (= number of groups j that i belongs to)
433 0 : if (m_bWeighted)
434 : {
435 0 : m_weight = c.clone_without_values();
436 : m_weight->set(0.0);
437 0 : for(int d = VERTEX; d <= VOLUME; ++d)
438 : {
439 : std::vector<std::vector<DoFIndex> >& vvDoFIndex = m_vDimCache[d].vvDoFIndex;
440 :
441 0 : for(size_t j = 0; j < vvDoFIndex.size(); ++j)
442 : {
443 :
444 : std::vector<DoFIndex>& vDoFIndex = vvDoFIndex[j];
445 : const size_t numIndex = vDoFIndex.size();
446 :
447 : // count number of connections for this row
448 0 : DoFRef(*m_weight, vDoFIndex[numIndex-1]) = 1.0;
449 0 : for (size_t i =0; i<numIndex-1; ++i)
450 0 : DoFRef(*m_weight, vDoFIndex[i]) += 1.0;
451 : }
452 :
453 : }
454 : }
455 :
456 : // extract local matrices
457 0 : for(int d = VERTEX; d <= VOLUME; ++d)
458 : {
459 : std::vector<std::vector<DoFIndex> >& vvDoFIndex = m_vDimCache[d].vvDoFIndex;
460 0 : std::vector<DenseMatrix<VariableArray2<number> > >& vBlockInverse = m_vDimCache[d].vBlockInverse;
461 :
462 0 : vBlockInverse.resize(vvDoFIndex.size());
463 0 : for(size_t b = 0; b < vvDoFIndex.size(); ++b)
464 : {
465 : // get storage
466 : std::vector<DoFIndex>& vDoFIndex = vvDoFIndex[b];
467 : DenseMatrix<VariableArray2<number> >& BlockInv = vBlockInverse[b];
468 :
469 : // get number of indices on patch
470 : const size_t numIndex = vDoFIndex.size();
471 : // std::cerr << "locSize" << numIndex << std::endl;
472 :
473 : // fill local block matrix
474 0 : BlockInv.resize(numIndex, numIndex);
475 0 : BlockInv = 0.0;
476 :
477 : // copy matrix rows (only including cols of selected indices)
478 0 : for (size_t i = 0; i < numIndex; i++)
479 : {
480 : // Diag (A)
481 0 : BlockInv(i,i) = m_alpha*DoFRef(A, vDoFIndex[i], vDoFIndex[i]);
482 0 : BlockInv(numIndex-1,i) = DoFRef(A, vDoFIndex[numIndex-1], vDoFIndex[i]);
483 0 : BlockInv(i,numIndex-1) = DoFRef(A, vDoFIndex[i], vDoFIndex[numIndex-1]);
484 :
485 : //for (size_t k = 0; k < numIndex; k++)
486 : // BlockInv(i,k) = DoFRef(A, vDoFIndex[i], vDoFIndex[k]);
487 : }
488 :
489 :
490 : // schur complement correction
491 :
492 :
493 : number schur = 0.0;
494 0 : if (m_weight.invalid())
495 : {
496 0 : for (size_t i = 0; i < numIndex-1; i++)
497 : {
498 0 : schur += BlockInv(numIndex-1,i)/BlockInv(i,i)*BlockInv(i,numIndex-1);
499 : }
500 : // std::cerr << "unweighted:" << schur << " " << BlockInv(numIndex-1, numIndex-1)<< m_beta<< "=>";
501 :
502 0 : BlockInv(numIndex-1, numIndex-1) = schur + (BlockInv(numIndex-1, numIndex-1)-schur)/m_beta;
503 :
504 : //std::cerr << BlockInv(numIndex-1, numIndex-1) << std::endl;
505 :
506 : } else {
507 0 : for (size_t i = 0; i < numIndex-1; i++)
508 : {
509 0 : number wi2 = DoFRef(*m_weight, vDoFIndex[i]);
510 0 : schur += wi2*BlockInv(numIndex-1,i)/BlockInv(i,i)*BlockInv(i,numIndex-1);
511 : }
512 : //std::cerr << "weighted:"<< schur << " " << BlockInv(numIndex-1, numIndex-1)<< m_beta<< "=>";
513 :
514 : // BlockInv(numIndex-1, numIndex-1) = - schur - (BlockInv(numIndex-1, numIndex-1)-schur)/m_beta; // worked with alpha=1, beta=-2.0
515 0 : BlockInv(numIndex-1, numIndex-1) = schur + (BlockInv(numIndex-1, numIndex-1) - schur)/m_beta;
516 :
517 : //std::cerr << BlockInv(numIndex-1, numIndex-1) << std::endl;
518 : }
519 :
520 :
521 : //BlockInv(numIndex-1, numIndex-1) = - 2.0*schur;
522 :
523 :
524 :
525 :
526 : // get inverse
527 0 : if(!Invert(BlockInv)){
528 0 : std::stringstream ss;
529 0 : ss << d <<"dim - Elem-Mat "<<b<<" singular. size: "<<numIndex<<"\n";
530 0 : for (size_t j = 0; j < numIndex; j++)
531 0 : ss << vDoFIndex[j][0] << ", ";
532 0 : ss << "\n";
533 :
534 0 : BlockInv = 0.0;
535 0 : for (size_t j = 0; j < numIndex; j++)
536 0 : for (size_t k = 0; k < numIndex; k++)
537 0 : BlockInv(j,k) = DoFRef(A, vDoFIndex[j], vDoFIndex[k]);
538 :
539 0 : ss << BlockInv;
540 :
541 0 : UG_THROW(ss.str());
542 0 : }
543 : }
544 : }
545 :
546 0 : m_bInit = true;
547 0 : }
548 :
549 : template <typename TDomain, typename TAlgebra>
550 0 : bool ComponentGaussSeidel<TDomain, TAlgebra>::
551 : step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
552 : {
553 : // check that grid funtion passed
554 : GridFunction<TDomain, TAlgebra>* pC
555 0 : = dynamic_cast<GridFunction<TDomain, TAlgebra>*>(&c);
556 0 : if(pC == NULL)
557 0 : UG_THROW("ComponentGaussSeidel: expects correction to be a GridFunction.");
558 :
559 : const vector_type* pD = &d;
560 : const matrix_type* pMat = pOp.get();
561 : #ifdef UG_PARALLEL
562 : SmartPtr<vector_type> spDtmp;
563 : if(pcl::NumProcs() > 1){
564 : // make defect unique
565 : spDtmp = d.clone();
566 : spDtmp->change_storage_type(PST_UNIQUE);
567 : pD = spDtmp.get();
568 : pMat = &m_A;
569 : }
570 : #endif
571 :
572 : // check if initialized
573 0 : if(!m_bInit)
574 0 : extract_blocks(*pMat, *pC);
575 :
576 : // set correction to zero
577 : pC->set(0.0);
578 :
579 : // loop Grouping objects
580 0 : for(size_t i = 0; i < m_vGroupObj.size(); ++i)
581 : {
582 : // get its dimension
583 0 : const int d = m_vGroupObj[i];
584 :
585 : // get relax param
586 0 : number damp = m_relax;
587 0 : if(d < (int)m_vDamp.size()) damp *= m_vDamp[d];
588 :
589 : // apply
590 0 : if (m_weight.invalid()) {
591 0 : apply_blocks(*pMat, *pC, *pD, damp, m_vDimCache[d], false);
592 0 : apply_blocks(*pMat, *pC, *pD, damp, m_vDimCache[d], true);
593 : } else {
594 0 : apply_blocks_weighted(*pMat, *pC, *pD, damp, m_vDimCache[d], false);
595 0 : apply_blocks_weighted(*pMat, *pC, *pD, damp, m_vDimCache[d], true);
596 : }
597 :
598 : }
599 :
600 : #ifdef UG_PARALLEL
601 : pC->set_storage_type(PST_UNIQUE);
602 : #endif
603 :
604 0 : return true;
605 : }
606 :
607 :
608 : template <typename TDomain, typename TAlgebra>
609 : SmartPtr<ILinearIterator<typename TAlgebra::vector_type> >
610 0 : ComponentGaussSeidel<TDomain, TAlgebra>::clone()
611 : {
612 : SmartPtr<ComponentGaussSeidel<TDomain, TAlgebra> >
613 0 : newInst(new ComponentGaussSeidel<TDomain, TAlgebra>(m_relax, m_vFullRowCmp, m_vGroupObj, m_vDamp));
614 0 : newInst->set_debug(debug_writer());
615 0 : newInst->set_damp(this->damping());
616 0 : newInst->set_alpha(this->m_alpha);
617 0 : newInst->set_beta(this->m_beta);
618 0 : newInst->set_weights(this->m_bWeighted);
619 0 : return newInst;
620 : }
621 :
622 :
623 : } // end namespace ug
624 :
625 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__COMPONENT_GAUSS_SEIDEL__ */
|
