Line data Source code
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 : /*
34 : * Gauss-Seidel type smoothers using alternating directions
35 : * for handling of anisotropic problems
36 : * Steps forward and backward in all coordinate directions
37 : */
38 :
39 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINE_SMOOTHERS__
40 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINE_SMOOTHERS__
41 :
42 : #include "common/util/smart_pointer.h"
43 : #include "lib_algebra/operator/interface/preconditioner.h"
44 :
45 : #include "lib_algebra/algebra_common/core_smoothers.h"
46 : #include "lib_disc/function_spaces/dof_position_util.h"
47 : #include "lib_disc/function_spaces/approximation_space.h"
48 : //#include "lib_disc/dof_manager/ordering/lexorder.h"
49 : #include "lib_disc/ordering_strategies/algorithms/lexorder.h"
50 : #include "lib_grid/algorithms/attachment_util.h"
51 : #include "lib_disc/common/groups_util.h"
52 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
53 : #ifdef UG_PARALLEL
54 : #include "lib_algebra/parallelization/parallelization.h"
55 : #include "lib_algebra/parallelization/parallel_matrix_overlap_impl.h"
56 : #endif
57 :
58 : namespace ug{
59 :
60 : // extern definition (in cpp file, avoiding conflicts)
61 : template<int dim>
62 : extern bool ComparePosDimYDir(const std::pair<MathVector<dim>, size_t> &p1,
63 : const std::pair<MathVector<dim>, size_t> &p2);
64 :
65 : template<int dim>
66 : extern bool ComparePosDimZDir(const std::pair<MathVector<dim>, size_t> &p1,
67 : const std::pair<MathVector<dim>, size_t> &p2);
68 :
69 :
70 : // ORDER IN Y DIRECTION
71 :
72 : template<int dim>
73 0 : void ComputeDirectionYOrder(std::vector<std::pair<MathVector<dim>, size_t> >& vPos,std::vector<size_t>& indY)
74 : {
75 0 : indY.resize(indY.size()+vPos.size());
76 : // sort indices based on their position
77 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimYDir<dim>);
78 : // write mapping
79 0 : for (size_t i=0; i < vPos.size(); ++i){
80 0 : indY[i] = vPos[i].second;
81 : };
82 0 : }
83 :
84 : template <typename TDomain>
85 0 : void OrderDirectionYForDofDist(SmartPtr<DoFDistribution> dd,
86 : ConstSmartPtr<TDomain> domain,std::vector<size_t>& indY)
87 : {
88 : // Lex Ordering is only possible in this cases:
89 : // b) Same number of DoFs on each geometric object (or no DoFs on object)
90 : // --> in this case we can order all dofs
91 : // a) different trial spaces, but DoFs for each trial spaces only on separate
92 : // geometric objects. (e.g. one space only vertices, one space only on edges)
93 : // --> in this case we can order all geometric objects separately
94 :
95 : // a) check for same number of DoFs on every geometric object
96 : bool bEqualNumDoFOnEachGeomObj = true;
97 : int numDoFOnGeomObj = -1;
98 0 : for(int si = 0; si < dd->num_subsets(); ++si){
99 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
100 0 : const int numDoF = dd->num_dofs((ReferenceObjectID)roid, si);
101 :
102 0 : if(numDoF == 0) continue;
103 :
104 0 : if(numDoFOnGeomObj == -1){
105 : numDoFOnGeomObj = numDoF;
106 : }
107 : else{
108 0 : if(numDoFOnGeomObj != numDoF)
109 : bEqualNumDoFOnEachGeomObj = false;
110 : }
111 : }
112 : }
113 :
114 : // b) check for non-mixed spaces
115 0 : std::vector<bool> bSingleSpaceUsage(NUM_REFERENCE_OBJECTS, true);
116 0 : std::vector<bool> vHasDoFs(NUM_REFERENCE_OBJECTS, false);
117 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
118 :
119 0 : LFEID lfeID = dd->local_finite_element_id(fct);
120 0 : const CommonLocalDoFSet& locDoF = LocalFiniteElementProvider::get_dofs(lfeID);
121 :
122 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
123 : const int numDoF = locDoF.num_dof((ReferenceObjectID)roid);
124 :
125 0 : if(numDoF <= 0) continue;
126 :
127 0 : if(vHasDoFs[roid] == false){
128 : vHasDoFs[roid] = true;
129 : }
130 : else{
131 : bSingleSpaceUsage[roid] = false;
132 : }
133 : }
134 : }
135 0 : std::vector<bool> bSortableComp(dd->num_fct(), true);
136 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
137 :
138 0 : LFEID lfeID = dd->local_finite_element_id(fct);
139 0 : const CommonLocalDoFSet& locDoF = LocalFiniteElementProvider::get_dofs(lfeID);
140 :
141 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
142 0 : if(locDoF.num_dof((ReferenceObjectID)roid) != 0){
143 0 : if(bSingleSpaceUsage[roid] == false)
144 : bSortableComp[fct] = false;
145 : }
146 : }
147 : }
148 :
149 : // get position attachment
150 : typedef typename std::pair<MathVector<TDomain::dim>, size_t> pos_type;
151 :
152 : // get positions of indices
153 : std::vector<pos_type> vPositions;
154 :
155 : // a) we can order globally
156 0 : if(bEqualNumDoFOnEachGeomObj)
157 : {
158 0 : ExtractPositions(domain, dd, vPositions);
159 :
160 0 : ComputeDirectionYOrder<TDomain::dim>(vPositions, indY);
161 : }
162 : // b) we can only order some spaces
163 : else
164 : {
165 : // UG_LOG("OrderLex: Cannot order globally, trying to order some components:\n");
166 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
167 0 : if(bSortableComp[fct] == false){
168 : // UG_LOG("OrderLex: "<<dd->name(fct)<<" NOT SORTED.\n");
169 0 : continue;
170 : }
171 :
172 0 : ExtractPositions(domain, dd, fct, vPositions);
173 :
174 0 : ComputeDirectionYOrder<TDomain::dim>(vPositions, indY);
175 : }
176 : }
177 0 : };
178 :
179 :
180 : // ORDER IN Z DIRECTION
181 :
182 :
183 : template<int dim>
184 0 : void ComputeDirectionZOrder(std::vector<std::pair<MathVector<dim>, size_t> >& vPos,std::vector<size_t>& indZ)
185 : {
186 0 : indZ.resize(indZ.size()+vPos.size());
187 0 : std::sort(vPos.begin(), vPos.end(), ComparePosDimZDir<dim>);
188 0 : for (size_t i=0; i < vPos.size(); ++i){
189 0 : indZ[i] = vPos[i].second;
190 : };
191 0 : }
192 :
193 : template <typename TDomain>
194 0 : void OrderDirectionZForDofDist(SmartPtr<DoFDistribution> dd,
195 : ConstSmartPtr<TDomain> domain,std::vector<size_t>& indZ)
196 : {
197 : // Lex Ordering is only possible in this cases:
198 : // b) Same number of DoFs on each geometric object (or no DoFs on object)
199 : // --> in this case we can order all dofs
200 : // a) different trial spaces, but DoFs for each trial spaces only on separate
201 : // geometric objects. (e.g. one space only vertices, one space only on edges)
202 : // --> in this case we can order all geometric objects separately
203 :
204 : // a) check for same number of DoFs on every geometric object
205 : bool bEqualNumDoFOnEachGeomObj = true;
206 : int numDoFOnGeomObj = -1;
207 0 : for(int si = 0; si < dd->num_subsets(); ++si){
208 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
209 0 : const int numDoF = dd->num_dofs((ReferenceObjectID)roid, si);
210 :
211 0 : if(numDoF == 0) continue;
212 :
213 0 : if(numDoFOnGeomObj == -1){
214 : numDoFOnGeomObj = numDoF;
215 : }
216 : else{
217 0 : if(numDoFOnGeomObj != numDoF)
218 : bEqualNumDoFOnEachGeomObj = false;
219 : }
220 : }
221 : }
222 :
223 : // b) check for non-mixed spaces
224 0 : std::vector<bool> bSingleSpaceUsage(NUM_REFERENCE_OBJECTS, true);
225 0 : std::vector<bool> vHasDoFs(NUM_REFERENCE_OBJECTS, false);
226 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
227 :
228 0 : LFEID lfeID = dd->local_finite_element_id(fct);
229 0 : const CommonLocalDoFSet& locDoF = LocalFiniteElementProvider::get_dofs(lfeID);
230 :
231 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
232 : const int numDoF = locDoF.num_dof((ReferenceObjectID)roid);
233 :
234 0 : if(numDoF <= 0) continue;
235 :
236 0 : if(vHasDoFs[roid] == false){
237 : vHasDoFs[roid] = true;
238 : }
239 : else{
240 : bSingleSpaceUsage[roid] = false;
241 : }
242 : }
243 : }
244 0 : std::vector<bool> bSortableComp(dd->num_fct(), true);
245 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
246 :
247 0 : LFEID lfeID = dd->local_finite_element_id(fct);
248 0 : const CommonLocalDoFSet& locDoF = LocalFiniteElementProvider::get_dofs(lfeID);
249 :
250 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid){
251 0 : if(locDoF.num_dof((ReferenceObjectID)roid) != 0){
252 0 : if(bSingleSpaceUsage[roid] == false)
253 : bSortableComp[fct] = false;
254 : }
255 : }
256 : }
257 :
258 : // get position attachment
259 : typedef typename std::pair<MathVector<TDomain::dim>, size_t> pos_type;
260 :
261 : // get positions of indices
262 : std::vector<pos_type> vPositions;
263 :
264 : // a) we can order globally
265 0 : if(bEqualNumDoFOnEachGeomObj)
266 : {
267 0 : ExtractPositions(domain, dd, vPositions);
268 :
269 : // get mapping: old -> new index
270 0 : ComputeDirectionZOrder<TDomain::dim>(vPositions, indZ);
271 : }
272 : // b) we can only order some spaces
273 : else
274 : {
275 : // UG_LOG("OrderLex: Cannot order globally, trying to order some components:\n");
276 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct){
277 0 : if(bSortableComp[fct] == false){
278 : // UG_LOG("OrderLex: "<<dd->name(fct)<<" NOT SORTED.\n");
279 0 : continue;
280 : }
281 :
282 0 : ExtractPositions(domain, dd, fct, vPositions);
283 :
284 : // get mapping: old -> new index
285 0 : ComputeDirectionZOrder<TDomain::dim>(vPositions, indZ);
286 : }
287 : }
288 0 : };
289 :
290 :
291 : template <typename TDomain, typename TBaseElem>
292 : void collectStretchedElementIndices(ConstSmartPtr<TDomain> domain,
293 : ConstSmartPtr<DoFDistribution> dd,
294 : std::vector<size_t>& indarray,number alpha){
295 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
296 : // vector for positions
297 : std::vector<MathVector<TDomain::dim> > vElemPos;
298 :
299 : // algebra indices vector
300 : std::vector<DoFIndex> ind;
301 :
302 : const typename TDomain::position_accessor_type& aaPos = domain->position_accessor();
303 :
304 : // loop all subsets
305 : for(int si = 0; si < dd->num_subsets(); ++si)
306 : {
307 : // get iterators
308 : iter = dd->template begin<TBaseElem>(si);
309 : iterEnd = dd->template end<TBaseElem>(si);
310 :
311 : // loop all elements
312 : for(;iter != iterEnd; ++iter)
313 : {
314 : // get vertex
315 : TBaseElem* elem = *iter;
316 :
317 : // loop all functions
318 : for(size_t fct = 0; fct < dd->num_fct(); ++fct)
319 : {
320 : // skip non-used function
321 : if(!dd->is_def_in_subset(fct,si)) continue;
322 : std::vector<Edge*> vEdge;
323 : CollectEdgesSorted(vEdge, domain->grid, elem);
324 : std::vector<number> edgeLength(vEdge.size());
325 : std::vector<DoFIndex> ind;
326 : MathVector<TDomain::dim> vCoord[2];
327 : for (size_t i=0;i<vEdge.size();i++){
328 : for (size_t j=0;j<2;j++) vCoord[i] = aaPos[vEdge[i]->vertex(i)];
329 : edgeLength[i] = VecDistance(vCoord[0],vCoord[1]);
330 : };
331 : number minedgelength=edgeLength[0];
332 : number maxedgelength=edgeLength[1];
333 : for (size_t i=1;i<vEdge.size();i++){
334 : if (edgeLength[i]<minedgelength) minedgelength=edgeLength[i];
335 : if (edgeLength[i]>maxedgelength) maxedgelength=edgeLength[i];
336 : };
337 : if (maxedgelength/minedgelength>alpha){
338 : // load indices associated with element function
339 : dd->inner_dof_indices(elem, fct, ind);
340 :
341 : // load positions associated with element and function
342 : InnerDoFPosition(vElemPos, elem, *(const_cast<TDomain*>(domain.get())),
343 : dd->local_finite_element_id(fct));
344 :
345 : // check correct size
346 : UG_ASSERT(ind.size() == vElemPos.size(), "Num MultiIndex ("<<ind.size()
347 : <<") and Num Position ("<<vElemPos.size()<<") must match."
348 : "GeomObject dim="<<geometry_traits<TBaseElem>::BASE_OBJECT_ID);
349 :
350 : // write position
351 : for(size_t sh = 0; sh < ind.size(); ++sh)
352 : {
353 : const size_t index = ind[sh][0];
354 : indarray.push_back(index);
355 : }
356 : };
357 : };
358 : };
359 : };
360 : };
361 :
362 :
363 :
364 : template <typename TDomain,typename TAlgebra>
365 : class LineGaussSeidel : public IPreconditioner<TAlgebra>
366 : {
367 : public:
368 : /// Domain
369 : typedef TDomain domain_type;
370 :
371 : // Algebra type
372 : typedef TAlgebra algebra_type;
373 :
374 : // Vector type
375 : typedef typename TAlgebra::vector_type vector_type;
376 :
377 : // Matrix type
378 : typedef typename TAlgebra::matrix_type matrix_type;
379 :
380 : /// Matrix Operator type
381 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
382 :
383 : /// Base type
384 : typedef IPreconditioner<TAlgebra> base_type;
385 :
386 : protected:
387 : using base_type::set_debug;
388 : using base_type::debug_writer;
389 : using base_type::write_debug;
390 :
391 : protected:
392 : // approximation space for level and surface grid
393 : SmartPtr<ApproximationSpace<TDomain> > m_spApproxSpace;
394 :
395 : // index sets
396 : std::vector<size_t> indY;
397 : std::vector<size_t> indZ;
398 : size_t m_ind_end;
399 :
400 : size_t m_nr_forwardx;
401 : size_t m_nr_backwardx;
402 : size_t m_nr_forwardy;
403 : size_t m_nr_backwardy;
404 : size_t m_nr_forwardz;
405 : size_t m_nr_backwardz;
406 :
407 : /// world dimension
408 : static const int dim = domain_type::dim;
409 :
410 : bool m_init;
411 :
412 : public:
413 : // Constructor
414 0 : LineGaussSeidel(SmartPtr<ApproximationSpace<TDomain> > approxSpace){
415 0 : m_spApproxSpace = approxSpace;
416 0 : m_init = false;
417 0 : m_nr_forwardx=1;
418 0 : m_nr_backwardx=1;
419 0 : m_nr_forwardy=1;
420 0 : m_nr_backwardy=1;
421 0 : m_nr_forwardz=1;
422 0 : m_nr_backwardz=1;
423 0 : OrderLex<TDomain>(*m_spApproxSpace,"lr");
424 0 : };
425 :
426 : // Update ordering indices
427 0 : bool update(size_t xsize){
428 0 : indY.resize(0);
429 0 : indZ.resize(0);
430 : // lexicographic ordering of unknowns
431 : if (m_nr_forwardx+m_nr_backwardx>0){
432 : // OrderLex<TDomain>(*m_spApproxSpace,"lr");
433 : }
434 0 : if (m_nr_forwardy+m_nr_backwardy+m_nr_forwardz+m_nr_backwardz>0){
435 : size_t level=3289578756;
436 0 : for (size_t i=0;i<m_spApproxSpace->num_levels();i++){
437 0 : if (m_spApproxSpace->dof_distribution(GridLevel(i, GridLevel::LEVEL, true))->num_indices()==xsize){
438 : level = i;
439 : break;
440 : };
441 : };
442 0 : if (level==3289578756){
443 : return false;
444 : }
445 0 : if ((dim>1)&&(m_nr_forwardy+m_nr_backwardy>0)){
446 0 : OrderDirectionYForDofDist<TDomain>(m_spApproxSpace->dof_distribution(GridLevel(level, GridLevel::LEVEL, true)), m_spApproxSpace->domain(),indY);
447 : }
448 0 : if ((dim>2)&&(m_nr_forwardz+m_nr_backwardz>0)){
449 0 : OrderDirectionZForDofDist<TDomain>(m_spApproxSpace->dof_distribution(GridLevel(level, GridLevel::LEVEL, true)), m_spApproxSpace->domain(),indZ);
450 : }
451 : };
452 0 : m_ind_end = indY.size();
453 0 : m_init = true;
454 0 : return true;
455 : }
456 :
457 : // Destructor
458 0 : ~LineGaussSeidel(){
459 : indY.clear();
460 : indZ.clear();
461 0 : };
462 :
463 0 : void set_num_steps(size_t forwardx,size_t backwardx,size_t forwardy,size_t backwardy,size_t forwardz,size_t backwardz){
464 0 : m_nr_forwardx=forwardx;
465 0 : m_nr_backwardx=backwardx;
466 0 : m_nr_forwardy=forwardy;
467 0 : m_nr_backwardy=backwardy;
468 0 : m_nr_forwardz=forwardz;
469 0 : m_nr_backwardz=backwardz;
470 0 : };
471 :
472 0 : void set_num_steps(size_t forwardx,size_t backwardx,size_t forwardy,size_t backwardy){
473 0 : m_nr_forwardx=forwardx;
474 0 : m_nr_backwardx=backwardx;
475 0 : m_nr_forwardy=forwardy;
476 0 : m_nr_backwardy=backwardy;
477 0 : m_nr_forwardz=0;
478 0 : m_nr_backwardz=0;
479 0 : };
480 :
481 0 : void set_num_steps(size_t forwardx,size_t backwardx){
482 0 : m_nr_forwardx=forwardx;
483 0 : m_nr_backwardx=backwardx;
484 0 : m_nr_forwardy=0;
485 0 : m_nr_backwardy=0;
486 0 : m_nr_forwardz=0;
487 0 : m_nr_backwardz=0;
488 0 : };
489 :
490 : // Clone
491 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
492 : {
493 0 : SmartPtr<LineGaussSeidel<domain_type,algebra_type> > newInst(new LineGaussSeidel<domain_type,algebra_type>(m_spApproxSpace));
494 0 : newInst->set_debug(debug_writer());
495 0 : newInst->set_damp(this->damping());
496 : newInst->set_num_steps(this->get_num_forwardx(),this->get_num_backwardx(),this->get_num_forwardy(),this->get_num_backwardy(),
497 : this->get_num_forwardz(),this->get_num_backwardz());
498 0 : return newInst;
499 : }
500 :
501 : /// returns if parallel solving is supported
502 0 : virtual bool supports_parallel() const {return true;}
503 :
504 : public:
505 0 : size_t get_num_forwardx(){return m_nr_forwardx;}
506 0 : size_t get_num_backwardx(){return m_nr_backwardx;}
507 0 : size_t get_num_forwardy(){return m_nr_forwardy;}
508 0 : size_t get_num_backwardy(){return m_nr_backwardy;}
509 0 : size_t get_num_forwardz(){return m_nr_forwardz;}
510 0 : size_t get_num_backwardz(){return m_nr_backwardz;}
511 :
512 : protected:
513 : // Name of preconditioner
514 0 : virtual const char* name() const {return "Line Gauss-Seidel";}
515 :
516 : // Preprocess routine
517 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
518 : {
519 : #ifdef UG_PARALLEL
520 : if(pcl::NumProcs() > 1)
521 : {
522 : // copy original matrix
523 : MakeConsistent(*pOp, m_A);
524 : // set zero on slaves
525 : std::vector<IndexLayout::Element> vIndex;
526 : CollectUniqueElements(vIndex, m_A.layouts()->slave());
527 : SetDirichletRow(m_A, vIndex);
528 : }
529 : #endif
530 0 : return true;
531 : }
532 :
533 : // Postprocess routine
534 0 : virtual bool postprocess() {return true;}
535 :
536 : // Stepping routine
537 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
538 : {
539 : #ifdef UG_PARALLEL
540 : if(pcl::NumProcs() > 1)
541 : {
542 : // make defect unique
543 : // todo: change that copying
544 : vector_type dhelp;
545 : dhelp.resize(d.size()); dhelp = d;
546 : dhelp.change_storage_type(PST_UNIQUE);
547 : if (!linegs_step(m_A, c, dhelp)) return false;
548 : // if(!gs_step_LL(m_A, c, dhelp)) return false;
549 : c.set_storage_type(PST_UNIQUE);
550 : return true;
551 : }
552 : else
553 : #endif
554 : {
555 0 : if(!linegs_step(*pOp, c, d)) return false;
556 : // if(!gs_step_LL(mat, c, d)) return false;
557 : #ifdef UG_PARALLEL
558 : c.set_storage_type(PST_UNIQUE);
559 : #endif
560 : return true;
561 : }
562 : }
563 :
564 : protected:
565 : #ifdef UG_PARALLEL
566 : matrix_type m_A;
567 : #endif
568 :
569 0 : bool linegs_step(const matrix_type &A, vector_type &x, const vector_type &b)
570 : {
571 : // gs LL has preconditioning matrix
572 : // (D-L)^{-1}
573 0 : if (m_init==false) update(x.size());
574 :
575 : size_t i;
576 :
577 : typename vector_type::value_type s;
578 :
579 : // forward in x direction
580 0 : if (m_nr_forwardx>0){
581 :
582 0 : for(i=0; i < x.size(); i++)
583 : {
584 0 : s = b[i];
585 :
586 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) && it.index() < i; ++it)
587 : // s -= it.value() * x[it.index()];
588 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
589 :
590 : // x[i] = s/A(i,i)
591 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
592 : }
593 0 : for (size_t count=1;count<m_nr_forwardx;count++){
594 0 : for(i=0; i < x.size(); i++)
595 : {
596 0 : s = b[i];
597 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it){
598 0 : if (it.index()==i) continue;
599 : // s -= it.value() * x[it.index()];
600 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
601 : }
602 : // x[i] = s/A(i,i)
603 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
604 : }
605 : };
606 : };
607 :
608 : // backward in x direction
609 0 : for (size_t count=0;count<m_nr_backwardx;count++){
610 0 : for (i=x.size()-1; (int)i>= 0; i--)
611 : {
612 0 : s = b[i];
613 :
614 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it){
615 0 : if (it.index()==i) continue;
616 : // s -= it.value() * x[it.index()];
617 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
618 : }
619 : // x[i] = s/A(i,i)
620 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
621 0 : if (i==0) break;
622 : };
623 : };
624 :
625 : if (dim==1) return true;
626 :
627 0 : if (m_nr_forwardy+m_nr_backwardy+m_nr_forwardz+m_nr_backwardz==0) return true;
628 :
629 : // forward in y direction
630 0 : for (size_t count=0;count<m_nr_forwardy;count++){
631 0 : for (size_t j=0;j < m_ind_end; j++){
632 0 : i = indY[j];
633 :
634 0 : s = b[i];
635 :
636 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it){
637 0 : if (it.index()==i) continue;
638 : // s -= it.value() * x[it.index()];
639 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
640 : }
641 : // x[i] = s/A(i,i)
642 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
643 : }
644 : };
645 :
646 : // backward in y direction
647 0 : for (size_t count=0;count<m_nr_backwardy;count++){
648 0 : for (size_t j=m_ind_end-1;(int)j >= 0; j--){
649 0 : i = indY[j];
650 :
651 0 : s = b[i];
652 :
653 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it){
654 0 : if (it.index()==i) continue;
655 : // s -= it.value() * x[it.index()];
656 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
657 : }
658 : // x[i] = s/A(i,i)
659 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
660 0 : if (j==0) break;
661 : }
662 : };
663 : if (dim==2) return true;
664 :
665 : // forward in z direction
666 0 : for (size_t count=0;count<m_nr_forwardz;count++){
667 0 : for (size_t j=0;j < m_ind_end; j++){
668 0 : i = indZ[j];
669 :
670 0 : s = b[i];
671 :
672 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
673 0 : if (it.index()==i) continue;
674 : // s -= it.value() * x[it.index()];
675 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
676 : };
677 : // x[i] = s/A(i,i)
678 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
679 : }
680 : }
681 :
682 : // backward in z direction
683 0 : for (size_t count=0;count<m_nr_backwardz;count++){
684 0 : for (size_t j=m_ind_end-1;(int)j >= 0; j--){
685 0 : i = indZ[j];
686 :
687 0 : s = b[i];
688 :
689 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
690 0 : if (it.index()==i) continue;
691 : // s -= it.value() * x[it.index()];
692 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
693 : };
694 : // x[i] = s/A(i,i)
695 : InverseMatMult(x[i], 1.0, A(i,i), s);
696 0 : if (j==0) break;
697 : }
698 : }
699 : return true;
700 : }
701 :
702 : };
703 :
704 : template <typename TDomain,typename TAlgebra>
705 : class LineVanka : public IPreconditioner<TAlgebra>
706 : {
707 : public:
708 : /// Domain
709 : typedef TDomain domain_type;
710 :
711 : // Algebra type
712 : typedef TAlgebra algebra_type;
713 :
714 : // Vector type
715 : typedef typename TAlgebra::vector_type vector_type;
716 :
717 : // Matrix type
718 : typedef typename TAlgebra::matrix_type matrix_type;
719 :
720 : /// Matrix Operator type
721 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
722 :
723 : /// Base type
724 : typedef IPreconditioner<TAlgebra> base_type;
725 :
726 : protected:
727 : using base_type::set_debug;
728 : using base_type::debug_writer;
729 : using base_type::write_debug;
730 :
731 : protected:
732 : // approximation space for level and surface grid
733 : SmartPtr<ApproximationSpace<TDomain> > m_spApproxSpace;
734 :
735 : // index sets
736 : std::vector<size_t> indY;
737 : std::vector<size_t> indZ;
738 : size_t m_ind_end;
739 :
740 : size_t m_nr_forwardx;
741 : size_t m_nr_backwardx;
742 : size_t m_nr_forwardy;
743 : size_t m_nr_backwardy;
744 : size_t m_nr_forwardz;
745 : size_t m_nr_backwardz;
746 :
747 : /// world dimension
748 : static const int dim = domain_type::dim;
749 :
750 : bool m_init;
751 :
752 : public:
753 : /// set m_relaxation parameter
754 0 : void set_relax(number omega){m_relax=omega;};
755 0 : number relax(){return m_relax;};
756 :
757 : protected:
758 : number m_relax;
759 :
760 : public:
761 : // Constructor
762 0 : LineVanka(SmartPtr<ApproximationSpace<TDomain> > approxSpace){
763 0 : m_spApproxSpace = approxSpace;
764 0 : m_init = false;
765 0 : m_nr_forwardx=1;
766 0 : m_nr_backwardx=1;
767 0 : m_nr_forwardy=1;
768 0 : m_nr_backwardy=1;
769 0 : m_nr_forwardz=1;
770 0 : m_nr_backwardz=1;
771 0 : m_relax=1;
772 : // OrderLex<TDomain>(*m_spApproxSpace,"lr");
773 0 : };
774 :
775 : // Update ordering indices
776 0 : bool update(size_t xsize){
777 0 : indY.resize(0);
778 0 : indZ.resize(0);
779 : // lexicographic ordering of unknowns
780 : if (m_nr_forwardx+m_nr_backwardx>0){
781 : // OrderLex<TDomain>(*m_spApproxSpace,"lr");
782 : }
783 0 : if (m_nr_forwardy+m_nr_backwardy+m_nr_forwardz+m_nr_backwardz>0){
784 : size_t level=3289578756;
785 0 : for (size_t i=0;i<m_spApproxSpace->num_levels();i++){
786 0 : if (m_spApproxSpace->dof_distribution(GridLevel(i, GridLevel::LEVEL, true))->num_indices()==xsize){
787 : level = i;
788 : break;
789 : };
790 : };
791 0 : if (level==3289578756){
792 : return false;
793 : }
794 0 : if ((dim>1)&&(m_nr_forwardy+m_nr_backwardy>0)){
795 0 : OrderDirectionYForDofDist<TDomain>(m_spApproxSpace->dof_distribution(GridLevel(level, GridLevel::LEVEL, true)), m_spApproxSpace->domain(),indY);
796 : }
797 0 : if ((dim>2)&&(m_nr_forwardz+m_nr_backwardz>0)){
798 0 : OrderDirectionZForDofDist<TDomain>(m_spApproxSpace->dof_distribution(GridLevel(level, GridLevel::LEVEL, true)), m_spApproxSpace->domain(),indZ);
799 : }
800 : };
801 0 : m_ind_end = indY.size();
802 0 : m_init = true;
803 0 : return true;
804 : }
805 :
806 : // Destructor
807 0 : ~LineVanka(){
808 : indY.clear();
809 : indZ.clear();
810 0 : };
811 :
812 0 : void set_num_steps(size_t forwardx,size_t backwardx,size_t forwardy,size_t backwardy,size_t forwardz,size_t backwardz){
813 0 : m_nr_forwardx=forwardx;
814 0 : m_nr_backwardx=backwardx;
815 0 : m_nr_forwardy=forwardy;
816 0 : m_nr_backwardy=backwardy;
817 0 : m_nr_forwardz=forwardz;
818 0 : m_nr_backwardz=backwardz;
819 0 : };
820 :
821 0 : void set_num_steps(size_t forwardx,size_t backwardx,size_t forwardy,size_t backwardy){
822 0 : m_nr_forwardx=forwardx;
823 0 : m_nr_backwardx=backwardx;
824 0 : m_nr_forwardy=forwardy;
825 0 : m_nr_backwardy=backwardy;
826 0 : m_nr_forwardz=0;
827 0 : m_nr_backwardz=0;
828 0 : };
829 :
830 0 : void set_num_steps(size_t forwardx,size_t backwardx){
831 0 : m_nr_forwardx=forwardx;
832 0 : m_nr_backwardx=backwardx;
833 0 : m_nr_forwardy=0;
834 0 : m_nr_backwardy=0;
835 0 : m_nr_forwardz=0;
836 0 : m_nr_backwardz=0;
837 0 : };
838 :
839 : // Clone
840 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
841 : {
842 0 : SmartPtr<LineVanka<domain_type,algebra_type> > newInst(new LineVanka<domain_type,algebra_type>(m_spApproxSpace));
843 0 : newInst->set_debug(debug_writer());
844 0 : newInst->set_damp(this->damping());
845 : newInst->set_num_steps(this->get_num_forwardx(),this->get_num_backwardx(),this->get_num_forwardy(),this->get_num_backwardy(),
846 : this->get_num_forwardz(),this->get_num_backwardz());
847 : newInst->set_relax(this->relax());
848 0 : return newInst;
849 : }
850 :
851 : /// returns if parallel solving is supported
852 0 : virtual bool supports_parallel() const {return true;}
853 :
854 : public:
855 0 : size_t get_num_forwardx(){return m_nr_forwardx;}
856 0 : size_t get_num_backwardx(){return m_nr_backwardx;}
857 0 : size_t get_num_forwardy(){return m_nr_forwardy;}
858 0 : size_t get_num_backwardy(){return m_nr_backwardy;}
859 0 : size_t get_num_forwardz(){return m_nr_forwardz;}
860 0 : size_t get_num_backwardz(){return m_nr_backwardz;}
861 :
862 : protected:
863 : // Name of preconditioner
864 0 : virtual const char* name() const {return "Line Vanka";}
865 :
866 : // Preprocess routine
867 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
868 : {
869 : #ifdef UG_PARALLEL
870 : if(pcl::NumProcs() > 1)
871 : {
872 : // copy original matrix
873 : MakeConsistent(*pOp, m_A);
874 : // set zero on slaves
875 : std::vector<IndexLayout::Element> vIndex;
876 : CollectUniqueElements(vIndex, m_A.layouts()->slave());
877 : SetDirichletRow(m_A, vIndex);
878 : }
879 : #endif
880 0 : return true;
881 : }
882 :
883 : // Postprocess routine
884 0 : virtual bool postprocess() {return true;}
885 :
886 : // Stepping routine
887 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
888 : {
889 : #ifdef UG_PARALLEL
890 : if(pcl::NumProcs() > 1)
891 : {
892 : // make defect unique
893 : // todo: change that copying
894 : vector_type dhelp;
895 : dhelp.resize(d.size()); dhelp = d;
896 : dhelp.change_storage_type(PST_UNIQUE);
897 : if (!linevanka_step(m_A, c, dhelp)) return false;
898 : // if(!gs_step_LL(m_A, c, dhelp)) return false;
899 : c.set_storage_type(PST_UNIQUE);
900 : return true;
901 : }
902 : else
903 : #endif
904 : {
905 0 : if(!linevanka_step(*pOp, c, d)) return false;
906 : // if(!gs_step_LL(mat, c, d)) return false;
907 : #ifdef UG_PARALLEL
908 : c.set_storage_type(PST_UNIQUE);
909 : #endif
910 : return true;
911 : }
912 : }
913 :
914 : protected:
915 : #ifdef UG_PARALLEL
916 : matrix_type m_A;
917 : #endif
918 :
919 0 : bool linevanka_step(const matrix_type &A, vector_type &x, const vector_type &b)
920 : {
921 : DenseVector< VariableArray1<number> > s;
922 : DenseVector< VariableArray1<number> > localx;
923 : DenseMatrix< VariableArray2<number> > mat;
924 :
925 : static const int MAXBLOCKSIZE = 19;
926 : size_t blockind[MAXBLOCKSIZE];
927 : size_t blocksize;
928 : // gs LL has preconditioning matrix
929 : // (D-L)^{-1}
930 0 : if (m_init==false) update(x.size());
931 :
932 : size_t i;
933 : size_t nrofelements=0;
934 :
935 0 : for(i=0; i < x.size(); i++)
936 : {
937 0 : x[i]=0;
938 0 : if (A(i,i)==0) nrofelements++;
939 : };
940 :
941 : // forward in x direction
942 0 : for (size_t count=0;count<m_nr_forwardx;count++){
943 0 : for(i=0; i < x.size(); i++){
944 0 : if (A(i,i)!=0) continue;
945 : blocksize=0;
946 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
947 0 : blockind[blocksize] = it.index();
948 0 : x[it.index()] = 0;
949 0 : blocksize++;
950 : };
951 0 : mat.resize(blocksize,blocksize);
952 0 : s.resize(blocksize);
953 0 : localx.resize(blocksize);
954 0 : for (size_t j=0;j<blocksize;j++){
955 : // fill local block matrix
956 0 : for (size_t k=0;k<blocksize;k++){
957 0 : mat.subassign(j,k,A(blockind[j],blockind[k]));
958 : };
959 : // compute rhs
960 0 : typename vector_type::value_type sj = b[blockind[j]];
961 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
962 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
963 : };
964 : s.subassign(j,sj);
965 : };
966 : // solve block
967 0 : InverseMatMult(localx,1,mat,s);
968 0 : for (size_t j=0;j<blocksize;j++){
969 0 : x[blockind[j]] = m_relax*localx[j];
970 : };
971 : };
972 : };
973 : // backward in x direction
974 0 : for (size_t count=0;count<m_nr_backwardx;count++){
975 0 : for (i=x.size()-1;(int)i>= 0; i--)
976 : {
977 0 : if (A(i,i)==0){
978 : blocksize=0;
979 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
980 0 : blockind[blocksize] = it.index();
981 0 : x[it.index()] = 0;
982 0 : blocksize++;
983 : };
984 0 : mat.resize(blocksize,blocksize);
985 0 : s.resize(blocksize);
986 0 : localx.resize(blocksize);
987 0 : for (size_t j=0;j<blocksize;j++){
988 : // fill local block matrix
989 0 : for (size_t k=0;k<blocksize;k++){
990 0 : mat.subassign(j,k,A(blockind[j],blockind[k]));
991 : };
992 : // compute rhs
993 0 : typename vector_type::value_type sj = b[blockind[j]];
994 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
995 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
996 : };
997 : s.subassign(j,sj);
998 : };
999 : // solve block
1000 0 : InverseMatMult(localx,1,mat,s);
1001 0 : for (size_t j=0;j<blocksize;j++){
1002 0 : x[blockind[j]] = m_relax*localx[j];
1003 : };
1004 : };
1005 0 : if (i==0) break;
1006 : };
1007 : };
1008 :
1009 : if (dim==1) return true;
1010 :
1011 0 : if (m_nr_forwardy+m_nr_backwardy+m_nr_forwardz+m_nr_backwardz==0) return true;
1012 :
1013 : // forward in y direction
1014 0 : for (size_t count=0;count<m_nr_forwardy;count++){
1015 0 : for (size_t sortedi=0;sortedi < m_ind_end; sortedi++){
1016 0 : i = indY[sortedi];
1017 : blocksize=0;
1018 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
1019 0 : blockind[blocksize] = it.index();
1020 0 : x[it.index()] = 0;
1021 0 : blocksize++;
1022 : };
1023 0 : mat.resize(blocksize,blocksize);
1024 0 : s.resize(blocksize);
1025 0 : localx.resize(blocksize);
1026 0 : for (size_t j=0;j<blocksize;j++){
1027 : // fill local block matrix
1028 0 : for (size_t k=0;k<blocksize;k++){
1029 0 : mat.subassign(j,k,A(blockind[j],blockind[k]));
1030 : };
1031 : // compute rhs
1032 0 : typename vector_type::value_type sj = b[blockind[j]];
1033 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
1034 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
1035 : };
1036 : s.subassign(j,sj);
1037 : };
1038 : // solve block
1039 0 : InverseMatMult(localx,1,mat,s);
1040 0 : for (size_t j=0;j<blocksize;j++){
1041 0 : x[blockind[j]] = m_relax*localx[j];
1042 : };
1043 : }
1044 : };
1045 :
1046 : // backward in y direction
1047 0 : for (size_t count=0;count<m_nr_backwardy;count++){
1048 0 : for (size_t sortedi=m_ind_end-1;(int)sortedi >= 0; sortedi--){
1049 0 : i = indY[sortedi];
1050 : blocksize=0;
1051 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
1052 0 : blockind[blocksize] = it.index();
1053 0 : x[it.index()] = 0;
1054 0 : blocksize++;
1055 : };
1056 0 : mat.resize(blocksize,blocksize);
1057 0 : s.resize(blocksize);
1058 0 : localx.resize(blocksize);
1059 0 : for (size_t j=0;j<blocksize;j++){
1060 : // fill local block matrix
1061 0 : for (size_t k=0;k<blocksize;k++){
1062 0 : mat.subassign(j,k,A(blockind[j],blockind[k]));
1063 : };
1064 : // compute rhs
1065 0 : typename vector_type::value_type sj = b[blockind[j]];
1066 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
1067 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
1068 : };
1069 : s.subassign(j,sj);
1070 : };
1071 : // solve block
1072 0 : InverseMatMult(localx,1,mat,s);
1073 0 : for (size_t j=0;j<blocksize;j++){
1074 0 : x[blockind[j]] = m_relax*localx[j];
1075 : };
1076 0 : if (sortedi==0) break;
1077 : }
1078 : };
1079 : if (dim==2) return true;
1080 :
1081 : // forward in z direction
1082 0 : for (size_t count=0;count<m_nr_forwardz;count++){
1083 0 : for (size_t sortedi=0;sortedi < m_ind_end; sortedi++){
1084 0 : i = indZ[sortedi];
1085 : blocksize=0;
1086 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
1087 0 : blockind[blocksize] = it.index();
1088 0 : x[it.index()] = 0;
1089 0 : blocksize++;
1090 : };
1091 0 : mat.resize(blocksize,blocksize);
1092 0 : s.resize(blocksize);
1093 0 : localx.resize(blocksize);
1094 0 : for (size_t j=0;j<blocksize;j++){
1095 : // fill local block matrix
1096 0 : for (size_t k=0;k<blocksize;k++){
1097 0 : mat.subassign(j,k,A(blockind[j],blockind[k]));
1098 : };
1099 : // compute rhs
1100 0 : typename vector_type::value_type sj = b[blockind[j]];
1101 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
1102 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
1103 : };
1104 : s.subassign(j,sj);
1105 : };
1106 : // solve block
1107 0 : InverseMatMult(localx,1,mat,s);
1108 0 : for (size_t j=0;j<blocksize;j++){
1109 0 : x[blockind[j]] = m_relax*localx[j];
1110 : };
1111 : }
1112 : }
1113 :
1114 : // backward in z direction
1115 0 : for (size_t count=0;count<m_nr_backwardz;count++){
1116 0 : for (size_t sortedi=m_ind_end-1;(int)sortedi >= 0; sortedi--){
1117 0 : i = indZ[sortedi];
1118 : blocksize=0;
1119 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) ; ++it){
1120 0 : blockind[blocksize] = it.index();
1121 0 : x[it.index()] = 0;
1122 0 : blocksize++;
1123 : };
1124 0 : mat.resize(blocksize,blocksize);
1125 0 : s.resize(blocksize);
1126 0 : localx.resize(blocksize);
1127 0 : for (size_t j=0;j<blocksize;j++){
1128 : // fill local block matrix
1129 0 : for (size_t k=0;k<blocksize;k++){
1130 0 : mat.subassign(j,k,A(blockind[j],blockind[k]));
1131 : };
1132 : // compute rhs
1133 0 : typename vector_type::value_type sj = b[blockind[j]];
1134 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(blockind[j]); it != A.end_row(blockind[j]) ; ++it){
1135 0 : MatMultAdd(sj, 1.0, sj, -1.0, it.value(), x[it.index()]);
1136 : };
1137 : s.subassign(j,sj);
1138 : };
1139 : // solve block
1140 0 : InverseMatMult(localx,1,mat,s);
1141 0 : for (size_t j=0;j<blocksize;j++){
1142 0 : x[blockind[j]] = m_relax*localx[j];
1143 : };
1144 0 : if (sortedi==0) break;
1145 : }
1146 : }
1147 : return true;
1148 : }
1149 :
1150 : };
1151 :
1152 : } // end namespace ug
1153 :
1154 : #endif // __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINE_SMOOTHERS__
|
