Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__LIB_ALGEBRA__LAPACK_AGGLOMERATING_LU_OPERATOR__
34 : #define __H__LIB_ALGEBRA__LAPACK_AGGLOMERATING_LU_OPERATOR__
35 : #include <iostream>
36 : #include <sstream>
37 :
38 : #include "lib_algebra/operator/interface/operator_inverse.h"
39 : #include "lib_algebra/operator/interface/matrix_operator_inverse.h"
40 : #include "lib_algebra/operator/interface/preconditioner.h"
41 :
42 : #ifdef UG_PARALLEL
43 : #include "lib_algebra/parallelization/collect_matrix.h"
44 : #include "lib_algebra/parallelization/parallelization.h"
45 : #include "lib_algebra/parallelization/parallelization_util.h"
46 : #endif
47 :
48 : namespace ug{
49 :
50 :
51 : #ifdef UG_PARALLEL
52 : template<typename T>
53 : void GatherVectorOnOne(HorizontalAlgebraLayouts &agglomerationLayout,
54 : ParallelVector<T> &collectedVec, const ParallelVector<T> &vec,
55 : ParallelStorageType type)
56 : {
57 : bool bRoot = pcl::ProcRank() == vec.layouts()->proc_comm().get_proc_id(0);
58 : GatherVectorOnOne(agglomerationLayout.master(), agglomerationLayout.slave(), agglomerationLayout.comm(),
59 : collectedVec, vec, type, bRoot);
60 : }
61 :
62 :
63 :
64 : template<typename T>
65 : void BroadcastVectorFromOne(HorizontalAlgebraLayouts &agglomerationLayout,
66 : ParallelVector<T> &vec, const ParallelVector<T> &collectedVec,
67 : ParallelStorageType type)
68 : {
69 : bool bRoot = pcl::ProcRank() == vec.layouts()->proc_comm().get_proc_id(0);
70 : BroadcastVectorFromOne(agglomerationLayout.master(), agglomerationLayout.slave(),
71 : agglomerationLayout.comm(), vec, collectedVec, type, bRoot);
72 : }
73 : #endif
74 :
75 : template <typename TBase, typename TAlgebra>
76 0 : class AgglomeratingBase : public TBase
77 : {
78 : public:
79 : // Algebra type
80 : typedef TAlgebra algebra_type;
81 :
82 : // Vector type
83 : typedef typename TAlgebra::vector_type vector_type;
84 :
85 : // Matrix type
86 : typedef typename TAlgebra::matrix_type matrix_type;
87 :
88 : public:
89 : // Destructor
90 0 : virtual ~AgglomeratingBase() {};
91 :
92 : bool i_am_root()
93 : {
94 : return m_bRoot;
95 : //m_bRoot = pcl::ProcRank() == agglomerationLayout.proc_comm().get_proc_id(0);;
96 : }
97 :
98 : bool empty()
99 : {
100 : return m_bEmpty;
101 : }
102 :
103 : bool init_mat(const matrix_type &A)
104 : {
105 : try{
106 : #ifdef UG_PARALLEL
107 : m_bEmpty = A.layouts()->proc_comm().empty();
108 : if(m_bEmpty) return true;
109 : m_bRoot = pcl::ProcRank() == A.layouts()->proc_comm().get_proc_id(0);
110 : PROFILE_FUNC();
111 :
112 : m_spCollectedOp = make_sp(new MatrixOperator<matrix_type, vector_type>());
113 : matrix_type &collectedA = m_spCollectedOp->get_matrix();
114 :
115 : CollectMatrixOnOneProc(A, collectedA, agglomerationLayout.master(), agglomerationLayout.slave());
116 : agglomerationLayout.comm() = A.layouts()->comm();
117 : agglomerationLayout.proc_comm() = A.layouts()->proc_comm();
118 :
119 : m_spLocalAlgebraLayouts = CreateLocalAlgebraLayouts();
120 : collectedA.set_layouts(m_spLocalAlgebraLayouts);
121 : #else
122 : m_bEmpty = false;
123 : m_bRoot = true;
124 : #endif
125 : }UG_CATCH_THROW("AgglomeratingBase::" << __FUNCTION__ << " failed")
126 : return true;
127 : }
128 :
129 : #ifdef UG_PARALLEL
130 : void init_collected_vec(vector_type &collectedX)
131 : {
132 : if(i_am_root())
133 : {
134 : matrix_type &collectedA = m_spCollectedOp->get_matrix();
135 : collectedX.resize(collectedA.num_rows());
136 : collectedX.set_layouts(m_spLocalAlgebraLayouts);
137 : }
138 : }
139 :
140 : void gather_vector_on_one(vector_type &collectedB, const vector_type &b, ParallelStorageType type)
141 : {
142 : GatherVectorOnOne(agglomerationLayout, collectedB, b, PST_ADDITIVE);
143 : }
144 :
145 : void broadcast_vector_from_one(vector_type &x, const vector_type &collectedX, ParallelStorageType type)
146 : {
147 : BroadcastVectorFromOne(agglomerationLayout, x, collectedX, PST_CONSISTENT);
148 : }
149 : #endif
150 :
151 : virtual bool init_agglomerated(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op) = 0;
152 : virtual bool apply_agglomerated(vector_type& x, const vector_type& b) = 0;
153 :
154 0 : bool is_serial()
155 : {
156 0 : UG_COND_THROW(m_pMatrix == NULL, "ERROR: No Matrix given.");
157 : #ifdef UG_PARALLEL
158 : return m_pMatrix->layouts()->proc_comm().is_local() || m_pMatrix->layouts()->proc_comm().empty();
159 :
160 : #else
161 0 : return true;
162 : #endif
163 : }
164 :
165 : /// Preprocess routine
166 0 : bool base_init(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
167 : {
168 : try{
169 : PROFILE_FUNC();
170 : // get matrix of Operator
171 0 : m_pMatrix = &Op->get_matrix();
172 :
173 :
174 : #ifdef UG_PARALLEL
175 : if(is_serial())
176 : return init_agglomerated(Op);
177 :
178 : PROFILE_BEGIN(AGGLOMERATING_Solver_BARRIER);
179 : m_pMatrix->layouts()->proc_comm().barrier();
180 : PROFILE_END();
181 :
182 : // init LU operator
183 : if(!init_mat(*m_pMatrix))
184 : {UG_LOG("ERROR in LUOperator::init: Cannot init LU Decomposition.\n"); return false;}
185 :
186 : bool bSuccess=true;
187 : if(i_am_root())
188 : {
189 : init_collected_vec(collectedX);
190 : init_collected_vec(collectedB);
191 : bSuccess = init_agglomerated(m_spCollectedOp);
192 : //UG_DLOG(LIB_ALG_LINEAR_SOLVER, 1,
193 : if(collectedX.size() != m_pMatrix->num_rows()) {
194 : UG_LOG("Agglomerated on proc 0. Size is " << collectedX.size() << "(was on this proc: "
195 : << m_pMatrix->num_rows() << ")\n");
196 : }
197 : }
198 : if(pcl::AllProcsTrue(bSuccess, m_pMatrix->layouts()->proc_comm()) == false) return false;
199 : return true;
200 : #else
201 0 : return init_agglomerated(Op);
202 : #endif
203 0 : }UG_CATCH_THROW("AgglomeratingBase::" << __FUNCTION__ << " failed")
204 : }
205 :
206 :
207 0 : virtual bool base_init(SmartPtr<ILinearOperator<vector_type> > A)
208 : {
209 : // cast operator
210 0 : SmartPtr<MatrixOperator<matrix_type,vector_type> > op =
211 : A.template cast_dynamic<MatrixOperator<matrix_type,vector_type> >();
212 :
213 : // check if correct types are present
214 0 : if(op.invalid())
215 0 : UG_THROW("IMatrixOperatorInverse::init:"
216 : " Passed operator is not matrix-based.");
217 :
218 : // forward request
219 0 : return base_init(op);
220 : }
221 :
222 0 : virtual bool apply(vector_type& x, const vector_type& b)
223 : {
224 : try{
225 :
226 0 : if(is_serial())
227 0 : return apply_agglomerated(x, b);
228 : #if UG_PARALLEL
229 : if(empty()) return true;
230 :
231 : gather_vector_on_one(collectedB, b, PST_ADDITIVE);
232 :
233 : collectedX.set(0.0);
234 : if(i_am_root())
235 : apply_agglomerated(collectedX, collectedB);
236 :
237 : broadcast_vector_from_one(x, collectedX, PST_CONSISTENT);
238 : #endif
239 0 : }UG_CATCH_THROW("AgglomeratingBase::" << __FUNCTION__ << " failed")
240 : // we're done
241 : return true;
242 : }
243 :
244 : // Compute u = L^{-1} * f AND return defect f := f - L*u
245 :
246 0 : virtual bool apply_return_defect(vector_type& u, vector_type& f)
247 : {
248 : PROFILE_FUNC();
249 : // solve u
250 0 : if(!apply(u, f)) return false;
251 :
252 : // calculate defect
253 : f.set(0.0);
254 0 : if(!m_pMatrix->matmul_minus(f, u))
255 : {
256 : UG_LOG("ERROR in 'LUSolver::apply_return_defect': Cannot apply matmul_minus.\n");
257 : return false;
258 : }
259 :
260 : // we're done
261 0 : return true;
262 : }
263 :
264 0 : virtual bool apply_update_defect(vector_type& u, vector_type& f)
265 : {
266 : PROFILE_FUNC();
267 : // solve u
268 0 : if(!apply(u, f)) return false;
269 : // update defect
270 0 : if(!m_pMatrix->matmul_minus(f, u))
271 : {
272 : UG_LOG("ERROR in 'LUSolver::apply_return_defect': Cannot apply matmul_minus.\n");
273 : return false;
274 : }
275 :
276 : // we're done
277 0 : return true;
278 : }
279 :
280 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
281 : {
282 0 : return apply(c, d);
283 : }
284 :
285 :
286 :
287 : protected:
288 : // Operator to invert
289 :
290 : // matrix to invert
291 : matrix_type* m_pMatrix;
292 : #ifdef UG_PARALLEL
293 : vector_type collectedB, collectedX;
294 : HorizontalAlgebraLayouts agglomerationLayout;
295 : SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spCollectedOp;
296 : SmartPtr<AlgebraLayouts> m_spLocalAlgebraLayouts;
297 : #endif
298 :
299 : bool m_bRoot;
300 : bool m_bEmpty;
301 :
302 : };
303 :
304 :
305 :
306 : template <typename TAlgebra>
307 : class AgglomeratingSolver : public
308 : AgglomeratingBase<IMatrixOperatorInverse< typename TAlgebra::matrix_type, typename TAlgebra::vector_type>, TAlgebra >
309 : {
310 : public:
311 : // Algebra type
312 : typedef TAlgebra algebra_type;
313 :
314 : // Vector type
315 : typedef typename TAlgebra::vector_type vector_type;
316 :
317 : // Matrix type
318 : typedef typename TAlgebra::matrix_type matrix_type;
319 :
320 : /// Base type
321 : typedef AgglomeratingBase<IMatrixOperatorInverse<matrix_type,vector_type>, algebra_type > base_type;
322 :
323 :
324 : public:
325 0 : AgglomeratingSolver(SmartPtr<ILinearOperatorInverse<vector_type, vector_type> > linOpInverse)
326 0 : {
327 0 : UG_COND_THROW(linOpInverse.valid()==false, "linOpInverse has to be != NULL");
328 0 : m_pLinOpInverse = linOpInverse;
329 0 : m_name = std::string("AgglomeratingSolver(") + linOpInverse->name() + ")";
330 0 : };
331 :
332 0 : virtual bool init_agglomerated(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
333 : {
334 : try{
335 0 : return m_pLinOpInverse->init(Op);
336 0 : }UG_CATCH_THROW("AgglomeratingSolver::" << __FUNCTION__ << " failed")
337 : }
338 0 : virtual bool apply_agglomerated(vector_type& x, const vector_type& b)
339 : {
340 : try{
341 0 : return m_pLinOpInverse->apply(x, b);
342 0 : }UG_CATCH_THROW("AgglomeratingSolver::" << __FUNCTION__ << " failed")
343 : }
344 :
345 : // Destructor
346 0 : virtual ~AgglomeratingSolver() {};
347 :
348 :
349 0 : virtual const char* name() const
350 : {
351 0 : return m_name.c_str();
352 : }
353 :
354 0 : virtual bool supports_parallel() const { return true; }
355 :
356 0 : virtual std::string config_string() const
357 : {
358 0 : return std::string("AgglomeratingSolver: ") + m_pLinOpInverse->config_string();
359 : }
360 :
361 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > A)
362 : {
363 0 : return base_type::base_init(A);
364 : }
365 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > A, const vector_type& u)
366 : {
367 0 : return base_type::base_init(A);
368 : }
369 0 : virtual bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
370 : {
371 0 : return base_type::base_init(Op);
372 : }
373 :
374 :
375 : protected:
376 : SmartPtr<ILinearOperatorInverse<vector_type, vector_type> > m_pLinOpInverse;
377 : std::string m_name;
378 : };
379 :
380 :
381 :
382 : template <typename TAlgebra>
383 : class AgglomeratingIterator : public
384 : AgglomeratingBase<ILinearIterator<typename TAlgebra::vector_type>, TAlgebra >
385 : {
386 : public:
387 : // Algebra type
388 : typedef TAlgebra algebra_type;
389 :
390 : // Vector type
391 : typedef typename TAlgebra::vector_type vector_type;
392 :
393 : // Matrix type
394 : typedef typename TAlgebra::matrix_type matrix_type;
395 :
396 : /// Base type
397 : typedef AgglomeratingBase<ILinearIterator<vector_type>, algebra_type > base_type;
398 :
399 : public:
400 0 : AgglomeratingIterator(SmartPtr<ILinearIterator<vector_type> > splinIt)
401 0 : {
402 0 : UG_COND_THROW(splinIt.valid()==false, "linOpInverse has to be != NULL");
403 0 : m_splinIt = splinIt;
404 0 : m_name = std::string("AgglomeratingIterator(") + splinIt->name() + std::string(")");
405 0 : }
406 :
407 0 : virtual bool init_agglomerated(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
408 : {
409 : try{
410 0 : return m_splinIt->init(Op);
411 0 : }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
412 : }
413 0 : virtual bool apply_agglomerated(vector_type& x, const vector_type& b)
414 : {
415 : try{
416 0 : return m_splinIt->apply(x, b);
417 0 : }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
418 : }
419 :
420 : // Destructor
421 0 : virtual ~AgglomeratingIterator() {};
422 :
423 :
424 0 : virtual const char* name() const
425 : {
426 0 : return m_name.c_str();
427 : }
428 :
429 0 : virtual bool supports_parallel() const { return true; }
430 :
431 : /// Clone
432 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
433 : {
434 0 : SmartPtr<ILinearIterator<vector_type> > linIt = m_splinIt->clone();
435 0 : return make_sp(new AgglomeratingIterator<algebra_type>(linIt));
436 : }
437 :
438 0 : virtual std::string config_string() const
439 : {
440 0 : return std::string("AgglomeratingIterator: ") + m_splinIt->config_string();
441 : }
442 :
443 :
444 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > A)
445 : {
446 0 : return base_type::base_init(A);
447 : }
448 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > A, const vector_type& u)
449 : {
450 0 : return base_type::base_init(A);
451 : }
452 :
453 :
454 : protected:
455 : SmartPtr<ILinearIterator<vector_type> > m_splinIt;
456 : std::string m_name;
457 : };
458 :
459 : template <typename TAlgebra>
460 : class AgglomeratingPreconditioner: public
461 : AgglomeratingBase<IPreconditioner<TAlgebra>, TAlgebra >
462 : {
463 : public:
464 : // Algebra type
465 : typedef TAlgebra algebra_type;
466 :
467 : // Vector type
468 : typedef typename TAlgebra::vector_type vector_type;
469 :
470 : // Matrix type
471 : typedef typename TAlgebra::matrix_type matrix_type;
472 :
473 : /// Base type
474 : typedef AgglomeratingBase<IPreconditioner<TAlgebra>, TAlgebra > base_type;
475 :
476 : public:
477 0 : AgglomeratingPreconditioner(SmartPtr<ILinearIterator<vector_type> > splinIt)
478 0 : {
479 0 : UG_COND_THROW(splinIt.valid()==false, "linOpInverse has to be != NULL");
480 0 : m_splinIt = splinIt;
481 0 : m_name = std::string("AgglomeratingPreconditioner(") + splinIt->name() + std::string(")");
482 0 : }
483 :
484 0 : virtual bool init_agglomerated(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
485 : {
486 : try{
487 0 : return m_splinIt->init(Op);
488 0 : }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
489 : }
490 0 : virtual bool apply_agglomerated(vector_type& x, const vector_type& b)
491 : {
492 : try{
493 0 : return m_splinIt->apply(x, b);
494 0 : }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
495 : }
496 :
497 : // Destructor
498 0 : virtual ~AgglomeratingPreconditioner() {};
499 :
500 :
501 0 : virtual const char* name() const
502 : {
503 0 : return m_name.c_str();
504 : }
505 :
506 0 : virtual bool supports_parallel() const { return true; }
507 :
508 : /// Clone
509 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
510 : {
511 0 : SmartPtr<ILinearIterator<vector_type> > linIt = m_splinIt->clone();
512 0 : return make_sp(new AgglomeratingIterator<algebra_type>(linIt));
513 : }
514 :
515 0 : virtual std::string config_string() const
516 : {
517 0 : return std::string("AgglomeratingPreconditioner: ") + m_splinIt->config_string();
518 : }
519 :
520 0 : virtual bool postprocess() { return true; }
521 0 : virtual bool preprocess(SmartPtr<MatrixOperator<typename TAlgebra::matrix_type, typename TAlgebra::vector_type> > pOp)
522 : {
523 0 : return base_type::base_init(pOp);
524 : }
525 :
526 :
527 : protected:
528 : SmartPtr<ILinearIterator<vector_type> > m_splinIt;
529 : std::string m_name;
530 : };
531 :
532 :
533 :
534 : } // end namespace ug
535 :
536 : #endif /* __H__LIB_ALGEBRA__LAPACK_LU_OPERATOR__ */
|
