Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
agglomerating_solver.h
Go to the documentation of this file.
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
41
42#ifdef UG_PARALLEL
46#endif
47
48namespace ug{
49
50
51#ifdef UG_PARALLEL
52template<typename T>
54 ParallelVector<T> &collectedVec, const ParallelVector<T> &vec,
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
64template<typename T>
66 ParallelVector<T> &vec, const ParallelVector<T> &collectedVec,
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
75template <typename TBase, typename TAlgebra>
76class 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 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
113 matrix_type &collectedA = m_spCollectedOp->get_matrix();
114
116 agglomerationLayout.comm() = A.layouts()->comm();
117 agglomerationLayout.proc_comm() = A.layouts()->proc_comm();
118
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
131 {
132 if(i_am_root())
133 {
134 matrix_type &collectedA = m_spCollectedOp->get_matrix();
135 collectedX.resize(collectedA.num_rows());
137 }
138 }
139
144
149#endif
150
152 virtual bool apply_agglomerated(vector_type& x, const vector_type& b) = 0;
153
155 {
156 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 return true;
162#endif
163 }
164
167 {
168 try{
169 PROFILE_FUNC();
170 // get matrix of Operator
171 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 {
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 return init_agglomerated(Op);
202#endif
203 }UG_CATCH_THROW("AgglomeratingBase::" << __FUNCTION__ << " failed")
204 }
205
206
208 {
209 // cast operator
211 A.template cast_dynamic<MatrixOperator<matrix_type,vector_type> >();
212
213 // check if correct types are present
214 if(op.invalid())
215 UG_THROW("IMatrixOperatorInverse::init:"
216 " Passed operator is not matrix-based.");
217
218 // forward request
219 return base_init(op);
220 }
221
222 virtual bool apply(vector_type& x, const vector_type& b)
223 {
224 try{
225
226 if(is_serial())
227 return apply_agglomerated(x, b);
228#if UG_PARALLEL
229 if(empty()) return true;
230
232
233 collectedX.set(0.0);
234 if(i_am_root())
236
238#endif
239 }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
247 {
248 PROFILE_FUNC();
249 // solve u
250 if(!apply(u, f)) return false;
251
252 // calculate defect
253 f.set(0.0);
254 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 return true;
262 }
263
265 {
266 PROFILE_FUNC();
267 // solve u
268 if(!apply(u, f)) return false;
269 // update defect
270 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 return true;
278 }
279
281 {
282 return apply(c, d);
283 }
284
285
286
287 protected:
288 // Operator to invert
289
290 // matrix to invert
292#ifdef UG_PARALLEL
297#endif
298
301
302};
303
304
305
306template <typename TAlgebra>
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
322
323
324 public:
326 {
327 UG_COND_THROW(linOpInverse.valid()==false, "linOpInverse has to be != NULL");
328 m_pLinOpInverse = linOpInverse;
329 m_name = std::string("AgglomeratingSolver(") + linOpInverse->name() + ")";
330 };
331
333 {
334 try{
335 return m_pLinOpInverse->init(Op);
336 }UG_CATCH_THROW("AgglomeratingSolver::" << __FUNCTION__ << " failed")
337 }
338 virtual bool apply_agglomerated(vector_type& x, const vector_type& b)
339 {
340 try{
341 return m_pLinOpInverse->apply(x, b);
342 }UG_CATCH_THROW("AgglomeratingSolver::" << __FUNCTION__ << " failed")
343 }
344
345 // Destructor
347
348
349 virtual const char* name() const
350 {
351 return m_name.c_str();
352 }
353
354 virtual bool supports_parallel() const { return true; }
355
356 virtual std::string config_string() const
357 {
358 return std::string("AgglomeratingSolver: ") + m_pLinOpInverse->config_string();
359 }
360
362 {
363 return base_type::base_init(A);
364 }
366 {
367 return base_type::base_init(A);
368 }
370 {
371 return base_type::base_init(Op);
372 }
373
374
375 protected:
377 std::string m_name;
378};
379
380
381
382template <typename TAlgebra>
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
398
399 public:
401 {
402 UG_COND_THROW(splinIt.valid()==false, "linOpInverse has to be != NULL");
403 m_splinIt = splinIt;
404 m_name = std::string("AgglomeratingIterator(") + splinIt->name() + std::string(")");
405 }
406
408 {
409 try{
410 return m_splinIt->init(Op);
411 }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
412 }
413 virtual bool apply_agglomerated(vector_type& x, const vector_type& b)
414 {
415 try{
416 return m_splinIt->apply(x, b);
417 }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
418 }
419
420 // Destructor
422
423
424 virtual const char* name() const
425 {
426 return m_name.c_str();
427 }
428
429 virtual bool supports_parallel() const { return true; }
430
437
438 virtual std::string config_string() const
439 {
440 return std::string("AgglomeratingIterator: ") + m_splinIt->config_string();
441 }
442
443
445 {
446 return base_type::base_init(A);
447 }
449 {
450 return base_type::base_init(A);
451 }
452
453
454 protected:
456 std::string m_name;
457};
458
459template <typename TAlgebra>
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
475
476 public:
478 {
479 UG_COND_THROW(splinIt.valid()==false, "linOpInverse has to be != NULL");
480 m_splinIt = splinIt;
481 m_name = std::string("AgglomeratingPreconditioner(") + splinIt->name() + std::string(")");
482 }
483
485 {
486 try{
487 return m_splinIt->init(Op);
488 }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
489 }
490 virtual bool apply_agglomerated(vector_type& x, const vector_type& b)
491 {
492 try{
493 return m_splinIt->apply(x, b);
494 }UG_CATCH_THROW("AgglomeratingIterator::" << __FUNCTION__ << " failed")
495 }
496
497 // Destructor
499
500
501 virtual const char* name() const
502 {
503 return m_name.c_str();
504 }
505
506 virtual bool supports_parallel() const { return true; }
507
514
515 virtual std::string config_string() const
516 {
517 return std::string("AgglomeratingPreconditioner: ") + m_splinIt->config_string();
518 }
519
520 virtual bool postprocess() { return true; }
525
526
527 protected:
529 std::string m_name;
530};
531
532
533
534} // end namespace ug
535
536#endif /* __H__LIB_ALGEBRA__LAPACK_LU_OPERATOR__ */
Definition smart_pointer.h:108
bool invalid() const
returns true if the pointer is invalid, false if not.
Definition smart_pointer.h:212
Definition agglomerating_solver.h:77
virtual bool apply(vector_type &x, const vector_type &b)
Definition agglomerating_solver.h:222
virtual bool apply_update_defect(vector_type &u, vector_type &f)
Definition agglomerating_solver.h:264
virtual ~AgglomeratingBase()
Definition agglomerating_solver.h:90
bool is_serial()
Definition agglomerating_solver.h:154
HorizontalAlgebraLayouts agglomerationLayout
Definition agglomerating_solver.h:294
bool m_bEmpty
Definition agglomerating_solver.h:300
void gather_vector_on_one(vector_type &collectedB, const vector_type &b, ParallelStorageType type)
Definition agglomerating_solver.h:140
TAlgebra algebra_type
Definition agglomerating_solver.h:80
TAlgebra::vector_type vector_type
Definition agglomerating_solver.h:83
SmartPtr< MatrixOperator< matrix_type, vector_type > > m_spCollectedOp
Definition agglomerating_solver.h:295
TAlgebra::matrix_type matrix_type
Definition agglomerating_solver.h:86
virtual bool init_agglomerated(SmartPtr< MatrixOperator< matrix_type, vector_type > > Op)=0
bool base_init(SmartPtr< MatrixOperator< matrix_type, vector_type > > Op)
Preprocess routine.
Definition agglomerating_solver.h:166
vector_type collectedB
Definition agglomerating_solver.h:293
SmartPtr< AlgebraLayouts > m_spLocalAlgebraLayouts
Definition agglomerating_solver.h:296
bool init_mat(const matrix_type &A)
Definition agglomerating_solver.h:103
virtual bool step(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp, vector_type &c, const vector_type &d)
Definition agglomerating_solver.h:280
matrix_type * m_pMatrix
Definition agglomerating_solver.h:291
virtual bool apply_agglomerated(vector_type &x, const vector_type &b)=0
vector_type collectedX
Definition agglomerating_solver.h:293
virtual bool apply_return_defect(vector_type &u, vector_type &f)
Definition agglomerating_solver.h:246
bool empty()
Definition agglomerating_solver.h:98
bool m_bRoot
Definition agglomerating_solver.h:299
void init_collected_vec(vector_type &collectedX)
Definition agglomerating_solver.h:130
void broadcast_vector_from_one(vector_type &x, const vector_type &collectedX, ParallelStorageType type)
Definition agglomerating_solver.h:145
bool i_am_root()
Definition agglomerating_solver.h:92
virtual bool base_init(SmartPtr< ILinearOperator< vector_type > > A)
Definition agglomerating_solver.h:207
Definition agglomerating_solver.h:385
virtual bool init(SmartPtr< ILinearOperator< vector_type > > A)
Definition agglomerating_solver.h:444
AgglomeratingIterator(SmartPtr< ILinearIterator< vector_type > > splinIt)
Definition agglomerating_solver.h:400
virtual const char * name() const
returns the name of iterator
Definition agglomerating_solver.h:424
virtual bool init(SmartPtr< ILinearOperator< vector_type > > A, const vector_type &u)
Definition agglomerating_solver.h:448
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition agglomerating_solver.h:429
AgglomeratingBase< ILinearIterator< vector_type >, algebra_type > base_type
Base type.
Definition agglomerating_solver.h:397
virtual std::string config_string() const
Definition agglomerating_solver.h:438
SmartPtr< ILinearIterator< vector_type > > m_splinIt
Definition agglomerating_solver.h:455
virtual ~AgglomeratingIterator()
Definition agglomerating_solver.h:421
virtual bool apply_agglomerated(vector_type &x, const vector_type &b)
Definition agglomerating_solver.h:413
virtual SmartPtr< ILinearIterator< vector_type > > clone()
Clone.
Definition agglomerating_solver.h:432
TAlgebra::vector_type vector_type
Definition agglomerating_solver.h:391
virtual bool init_agglomerated(SmartPtr< MatrixOperator< matrix_type, vector_type > > Op)
Definition agglomerating_solver.h:407
TAlgebra::matrix_type matrix_type
Definition agglomerating_solver.h:394
std::string m_name
Definition agglomerating_solver.h:456
TAlgebra algebra_type
Definition agglomerating_solver.h:388
Definition agglomerating_solver.h:462
TAlgebra::matrix_type matrix_type
Definition agglomerating_solver.h:471
virtual const char * name() const
returns the name of iterator
Definition agglomerating_solver.h:501
AgglomeratingPreconditioner(SmartPtr< ILinearIterator< vector_type > > splinIt)
Definition agglomerating_solver.h:477
virtual bool postprocess()
cleans the operator
Definition agglomerating_solver.h:520
SmartPtr< ILinearIterator< vector_type > > m_splinIt
Definition agglomerating_solver.h:528
virtual ~AgglomeratingPreconditioner()
Definition agglomerating_solver.h:498
std::string m_name
Definition agglomerating_solver.h:529
virtual bool preprocess(SmartPtr< MatrixOperator< typename TAlgebra::matrix_type, typename TAlgebra::vector_type > > pOp)
initializes the preconditioner
Definition agglomerating_solver.h:521
virtual bool apply_agglomerated(vector_type &x, const vector_type &b)
Definition agglomerating_solver.h:490
virtual SmartPtr< ILinearIterator< vector_type > > clone()
Clone.
Definition agglomerating_solver.h:509
TAlgebra algebra_type
Definition agglomerating_solver.h:465
TAlgebra::vector_type vector_type
Definition agglomerating_solver.h:468
virtual bool init_agglomerated(SmartPtr< MatrixOperator< matrix_type, vector_type > > Op)
Definition agglomerating_solver.h:484
AgglomeratingBase< IPreconditioner< TAlgebra >, TAlgebra > base_type
Base type.
Definition agglomerating_solver.h:474
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition agglomerating_solver.h:506
virtual std::string config_string() const
Definition agglomerating_solver.h:515
Definition agglomerating_solver.h:309
virtual bool init_agglomerated(SmartPtr< MatrixOperator< matrix_type, vector_type > > Op)
Definition agglomerating_solver.h:332
virtual bool init(SmartPtr< ILinearOperator< vector_type > > A, const vector_type &u)
Definition agglomerating_solver.h:365
virtual ~AgglomeratingSolver()
Definition agglomerating_solver.h:346
virtual bool apply_agglomerated(vector_type &x, const vector_type &b)
Definition agglomerating_solver.h:338
virtual bool init(SmartPtr< ILinearOperator< vector_type > > A)
Definition agglomerating_solver.h:361
virtual bool init(SmartPtr< MatrixOperator< matrix_type, vector_type > > Op)
Definition agglomerating_solver.h:369
SmartPtr< ILinearOperatorInverse< vector_type, vector_type > > m_pLinOpInverse
Definition agglomerating_solver.h:376
virtual std::string config_string() const
returns information about configuration parameters
Definition agglomerating_solver.h:356
TAlgebra::vector_type vector_type
Definition agglomerating_solver.h:315
std::string m_name
Definition agglomerating_solver.h:377
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition agglomerating_solver.h:354
AgglomeratingBase< IMatrixOperatorInverse< matrix_type, vector_type >, algebra_type > base_type
Base type.
Definition agglomerating_solver.h:321
AgglomeratingSolver(SmartPtr< ILinearOperatorInverse< vector_type, vector_type > > linOpInverse)
Definition agglomerating_solver.h:325
TAlgebra algebra_type
Definition agglomerating_solver.h:312
TAlgebra::matrix_type matrix_type
Definition agglomerating_solver.h:318
virtual const char * name() const
returns the name of the operator inverse
Definition agglomerating_solver.h:349
Holds Interfaces and communicators for horizontal communication.
Definition algebra_layouts.h:48
const pcl::ProcessCommunicator & proc_comm() const
returns process communicator
Definition algebra_layouts.h:68
const IndexLayout & master() const
Definition algebra_layouts.h:61
pcl::InterfaceCommunicator< IndexLayout > & comm() const
returns (non-const !!!) communicator
Definition algebra_layouts.h:78
const IndexLayout & slave() const
Definition algebra_layouts.h:63
describes a linear iterator
Definition linear_iterator.h:81
describes a linear mapping X->Y
Definition linear_operator.h:80
describes an inverse linear mapping X->Y
Definition linear_operator_inverse.h:80
Definition matrix_operator.h:49
virtual M & get_matrix()
Definition matrix_operator.h:74
Definition parallel_vector.h:60
ConstSmartPtr< AlgebraLayouts > layouts() const
returns the algebra layouts
Definition parallel_vector.h:97
ParallelStorageType
Definition parallel_storage_type.h:66
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
@ PST_ADDITIVE
Definition parallel_storage_type.h:69
int ProcRank()
returns the rank of the process
Definition pcl_base.cpp:83
bool AllProcsTrue(bool bFlag, ProcessCommunicator comm)
Definition pcl_util.cpp:54
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
the ug namespace
void CollectMatrixOnOneProc(const matrix_type &A, matrix_type &collectedA, IndexLayout &masterLayout, IndexLayout &slaveLayout)
Definition collect_matrix.h:246
void GatherVectorOnOne(HorizontalAlgebraLayouts &agglomerationLayout, ParallelVector< T > &collectedVec, const ParallelVector< T > &vec, ParallelStorageType type)
Definition agglomerating_solver.h:53
void BroadcastVectorFromOne(HorizontalAlgebraLayouts &agglomerationLayout, ParallelVector< T > &vec, const ParallelVector< T > &collectedVec, ParallelStorageType type)
Definition agglomerating_solver.h:65
SmartPtr< AlgebraLayouts > CreateLocalAlgebraLayouts()
Definition parallelization_util.cpp:936
#define PROFILE_BEGIN(name)
Definition profiler.h:254
#define PROFILE_END()
Definition profiler.h:256
#define PROFILE_FUNC()
Definition profiler.h:257
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836