ug4
Loading...
Searching...
No Matches
matrix_overlap_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2017: G-CSC, Goethe University Frankfurt
3 * Author: Sebastian Reiter
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__UG_matrix_overlap_impl
34#define __H__UG_matrix_overlap_impl
35
36#include "algebra_layouts.h"
39
40namespace ug{
41
42template <class TLayout>
43void GetLayoutTargetProcs(std::vector<int>& procsOut, const TLayout& layout)
44{
45 procsOut.clear();
46 typedef typename TLayout::const_iterator const_iter_t;
47 for(const_iter_t iter = layout.begin(); iter != layout.end(); ++iter)
48 {
49 procsOut.push_back(layout.interface(iter).get_target_proc());
50 }
51}
52
53
54template <class TMatrix>
56{
57 public:
58 ComPol_MatCopyDiag(TMatrix& mat): m_mat(mat) {}
59
60 int get_required_buffer_size (const Interface& interface)
61 {
63 return interface.size() * sizeof(typename TMatrix::value_type);
64 else
65 return -1;
66 }
67
68 bool collect (ug::BinaryBuffer& buff, const Interface& interface)
69 {
70 PROFILE_BEGIN_GROUP(ComPol_MatCopyDiag_collect, "algebra parallelization");
71
72 for(typename Interface::const_iterator iter = interface.begin();
73 iter != interface.end(); ++iter)
74 {
75 const size_t index = interface.get_element(iter);
76 Serialize(buff, m_mat(index, index));
77 }
78 return true;
79 }
80
81 bool extract (ug::BinaryBuffer& buff, const Interface& interface)
82 {
83 PROFILE_BEGIN_GROUP(ComPol_MatCopyDiag_extract, "algebra parallelization");
84
85 for(typename Interface::const_iterator iter = interface.begin();
86 iter != interface.end(); ++iter)
87 {
88 const size_t index = interface.get_element(iter);
89 Deserialize(buff, m_mat(index, index));
90 }
91 return true;
92 }
93
94 private:
95 TMatrix& m_mat;
96};
97
98
99
101
111template <class TMatrix>
113 : public pcl::ICommunicationPolicy<IndexLayout>
114{
115 public:
116 typedef typename TMatrix::value_type block_t;
118
121 ComPol_MatCreateOverlap(TMatrix& rMat, AlgebraIDVec& vGlobalID)
122 : m_mat(rMat), m_globalIDs(vGlobalID)
123 {
124 UG_COND_THROW(vGlobalID.size() < m_mat.num_rows(), "Not enough GlobalIDs");
125 // fill the map global->local
127 }
128
130 virtual bool collect(ug::BinaryBuffer& buff, const Interface& interface)
131 {
132 PROFILE_BEGIN_GROUP(ComPol_MatAddRowsOverlap0_collect, "algebra parallelization");
133 typedef typename TMatrix::row_iterator row_iterator;
134
135 // loop interface
136 for(typename Interface::const_iterator iter = interface.begin();
137 iter != interface.end(); ++iter)
138 {
139 // get index
140 const size_t index = interface.get_element(iter);
141
142 // count number of row entries
143 const row_iterator rowEnd = m_mat.end_row(index);
144 size_t numRowEntry = 0;
145 for(row_iterator it_k = m_mat.begin_row(index); it_k != rowEnd; ++it_k)
146 numRowEntry++;
147
148 // write number of row entries to stream
149 Serialize(buff, numRowEntry);
150
151 // write entries and global id to stream
152 for(row_iterator it_k = m_mat.begin_row(index); it_k != rowEnd; ++it_k)
153 {
154 const size_t k = it_k.index();
155 block_t& a_ik = it_k.value();
156
157 // write global entry to buffer
158 Serialize(buff, m_globalIDs[k]);
159
160 // write entry into buffer
161 Serialize(buff, a_ik);
162 }
163 }
164
166 return true;
167 }
168
169
171 virtual bool extract(ug::BinaryBuffer& buff, const Interface& interface)
172 {
173 PROFILE_BEGIN_GROUP(ComPol_MatAddRowsOverlap0_extract, "algebra parallelization");
174
175 // we'll read global ids into this variable
176 AlgebraID gID;
177
178 // we'll read blocks into this var
179 block_t block;
180
181 const int targetProc = interface.get_target_proc ();
182
183 // loop interface
184 for(typename Interface::const_iterator iter = interface.begin();
185 iter != interface.end(); ++iter)
186 {
187 // get index
188 const size_t index = interface.get_element(iter);
189
190 // read the number of connections
191 size_t numConnections = 0;
192 Deserialize(buff, numConnections);
193
194 // read each connection
195 for(size_t i_conn = 0; i_conn < numConnections; ++i_conn){
196 Deserialize(buff, gID);
197 Deserialize(buff, block);
198
199 // if gID exists on this process, then set the connection to
200 // the new value.
201 size_t conInd;
202 if(m_algIDHash.get_entry(conInd, gID)){
203 m_mat(index, conInd) += block;
204 }
205 else {
206 const ExtCon ec(index, m_globalIDs[index], gID, targetProc);
207 m_recvNewCons [ec] = block;
208 m_recvNewIDs.insert (gID);
209 }
210 }
211 }
212
214 return true;
215 }
216
217
219 virtual void post_process()
220 {
221 // create new master overlap
222 using namespace std;
223
224 // we create a new layout, since the old one may be shared between many
225 // different vectors and matrices. H-Master and H-Slave layouts are still
226 // the same.
227 SmartPtr <AlgebraLayouts> newLayout = make_sp(new AlgebraLayouts);
228 *newLayout = *m_mat.layouts();
229 newLayout->enable_overlap(true);
230 m_mat.set_layouts (newLayout);
231
232 const size_t oldSize = m_mat.num_rows();
233 const size_t numNewInds = m_recvNewIDs.size();
234 const size_t newSize = oldSize + numNewInds;
235
236 // add new entries to the algebra hash and to the globalID array
237 {
238 m_globalIDs.reserve(newSize);
239 size_t i = oldSize;
240 for(set<AlgebraID>::iterator iter = m_recvNewIDs.begin();
241 iter != m_recvNewIDs.end(); ++iter, ++i)
242 {
243 m_algIDHash.insert(*iter, i);
244 m_globalIDs.push_back(*iter);
245 }
246 }
247
248 if(newSize != oldSize) {
249 // For each new DoF we set a DirichletRow
250 m_mat.resize_and_keep_values (newSize, newSize);
251 for(size_t i = oldSize; i < newSize; ++i){
253 }
254 }
255
256
257 vector<int> slaveProcs; // process ranks of processes with associated slave interfaces
258 GetLayoutTargetProcs(slaveProcs, newLayout->master());
259
260 // the following buffers and vectors are used to collect globalIDs of
261 // newly created entries sorted by the slave process from which they
262 // were received.
263 BinaryBuffer sendBuf;
264 // size of the message for the i-th process in slaveProcs
265 vector<int> msgSizeForSlaveProcs(slaveProcs.size(), 0);
266 int slaveInd = -1;
267
268 // create new master-overlap and add connections to matrix
269 {
270 IndexLayout::Interface* itfc = NULL;
271 vector<bool> added;// keeps track of whether an index was already added.
272
273 for (typename map<ExtCon, block_t>::iterator iter = m_recvNewCons.begin();
274 iter != m_recvNewCons.end(); ++iter)
275 {
276 const ExtCon& curExtCon = iter->first;
277 const int targetProc = curExtCon.conProc;
278 if (!itfc || targetProc != itfc->get_target_proc()){
279 itfc = &newLayout->master_overlap().interface(targetProc);
280 added.clear();
281 added.resize(numNewInds, false);
282
283 // find the index of the corresponding slave proc
284 slaveInd = -1;
285 for(size_t islave = 0; islave < slaveProcs.size(); ++islave){
286 if(slaveProcs[islave] == targetProc){
287 slaveInd = int(islave);
288 break;
289 }
290 }
291 UG_COND_THROW(slaveInd == -1, "slaveProcs doas not contain referenced slave rank");
292 }
293
294 size_t toInd;
295 UG_COND_THROW(!m_algIDHash.get_entry(toInd, curExtCon.toID),
296 "Expected AlgebraID not found in Hash");
297
298 if(!added[toInd - oldSize]){
299 itfc->push_back(toInd);
300 const size_t oldWritePos = sendBuf.write_pos();
301 Serialize(sendBuf, curExtCon.toID);
302 msgSizeForSlaveProcs[slaveInd] += int(sendBuf.write_pos() - oldWritePos);
303 }
304 added[toInd - oldSize] = true;
305
306 m_mat(curExtCon.fromInd, toInd) += iter->second;
307 }
308 }
309
310 // master processing done!
311 // now find all processes which contain master interfaces to local slave interfaces
312 vector<int> masterProcs;
313 GetLayoutTargetProcs(masterProcs, newLayout->slave());
314
315 vector<int> recvSizes (masterProcs.size());
316 BinaryBuffer recvBuf;
317
318 newLayout->proc_comm().distribute_data(
319 recvBuf, &recvSizes.front(), &masterProcs.front(),
320 (int)masterProcs.size(),
321 sendBuf.buffer(), &msgSizeForSlaveProcs.front(),
322 &slaveProcs.front(), (int)slaveProcs.size());
323
324 // Extract content of recvBuf and create new slave-overlap
325 for(size_t i = 0; i < masterProcs.size(); ++i)
326 {
327 IndexLayout::Interface& itfc = newLayout->slave_overlap()
328 .interface(masterProcs[i]);
329
330 size_t oldReadPos = recvBuf.read_pos();
331 while(recvBuf.read_pos() < oldReadPos + recvSizes[i]){
332 AlgebraID globID;
333 Deserialize(recvBuf, globID);
334 size_t locID;
335 if(m_algIDHash.get_entry(locID, globID)){
336 itfc.push_back(locID);
337 }
338 else{
339 UG_THROW("GlobalID " << globID << " expected on this process "
340 "but not found.");
341 }
342 }
343 }
344
345 {
346 // Make matrix partialy consistent on slave interfaces
348 newLayout->comm().send_data(newLayout->master(), comPolMatCopy);
349 newLayout->comm().receive_data(newLayout->slave(), comPolMatCopy);
350 newLayout->comm().communicate();
351
352 //WARNING: THE CODE BELOW DOESN'T WORK WELL E.G. WITH ILU-OVERLAP
353 // INSTEAD ONLY THE DIAGONAL ENTRIES ARE COPIED TO MASTER_OVERLAP NODES
354 // (SEE BELOW)
355 // MAYBE SOME INVESTIGATIONS WOULD BE NICE REGARDING THIS BEHAVIOR...
356 // Make matrix partialy consistent on h-master-overlap
357 // newLayout->comm().send_data(newLayout->slave_overlap(), comPolMatCopy);
358 // newLayout->comm().receive_data(newLayout->master_overlap(), comPolMatCopy);
359 // newLayout->comm().communicate();
360 }
361
362 // {
363 //NOTE: THIS CODE IS DISABLED BECAUSE THE MATRIX IS ALREADY MADE
364 // PARTIALLY CONSISTENT ABOVE. IT CURRENTLY REMAINS FOR DEBUG PURPOSES
365 // // set h-slave rows to diagonal-only
366 // IndexLayout& slaveLayout = newLayout->slave();
367 // for(size_t lvl = 0; lvl < slaveLayout.num_levels(); ++lvl){
368 // for(IndexLayout::const_iterator interfaceIter = slaveLayout.begin(lvl);
369 // interfaceIter != slaveLayout.end(lvl); ++interfaceIter)
370 // {
371 // const IndexLayout::Interface& interface = slaveLayout.interface(interfaceIter);
372 // for(IndexLayout::Interface::const_iterator iter = interface.begin();
373 // iter != interface.end(); ++iter)
374 // {
375 // const size_t index = interface.get_element(iter);
376 // typename TMatrix::value_type diag = m_mat(index, index);
377 // SetDirichletRow(m_mat, index);
378 // m_mat(index, index) = diag;
379 // }
380 // }
381 // }
382
383 {
384 // Copy diagonal entries to master-overlap entries
385 ComPol_MatCopyDiag<TMatrix> comPolMatCopyDiag(m_mat);
386 newLayout->comm().send_data(newLayout->slave_overlap(), comPolMatCopyDiag);
387 newLayout->comm().receive_data(newLayout->master_overlap(), comPolMatCopyDiag);
388 newLayout->comm().communicate();
389 }
390 }
391
392 private:
393 TMatrix& m_mat;
394
397
400
401 struct ExtCon {
402 ExtCon (size_t _fromInd, const AlgebraID& _fromID,
403 const AlgebraID& _toID, int _conProc) :
404 fromInd (_fromInd),
405 fromID (_fromID),
406 toID (_toID),
407 conProc (_conProc)
408 {}
409
410 size_t fromInd;
414
415 bool operator < (const ExtCon& ec) const
416 {
417 if(conProc < ec.conProc) return true;
418 else if(conProc > ec.conProc) return false;
419 if(toID < ec.toID) return true;
420 else if(toID > ec.toID) return false;
421 return (fromID < ec.fromID);
422 }
423 };
424
426 std::map<ExtCon, block_t> m_recvNewCons;
428 std::set<AlgebraID> m_recvNewIDs;
429};
430
431
432
433template <class TMatrix>
434void CreateOverlap (TMatrix& mat)
435{
436 using namespace std;
437 PROFILE_FUNC_GROUP("algebra parallelization");
438
439 vector<AlgebraID> globalIDs;
440
441 GenerateGlobalAlgebraIDs(mat.layouts()->comm(),
442 globalIDs,
443 mat.num_rows(),
444 mat.layouts()->master(),
445 mat.layouts()->slave());
446
447 ComPol_MatCreateOverlap<TMatrix> comPolOverlap(mat, globalIDs);
448 mat.layouts()->comm().send_data(mat.layouts()->slave(), comPolOverlap);
449 mat.layouts()->comm().receive_data(mat.layouts()->master(), comPolOverlap);
450 mat.layouts()->comm().communicate();
451
452 comPolOverlap.post_process();
453
454// todo: Remove this once Overlap creation is stable. Check for redistributed grids
455// which have h-masters and h-slaves on one process!
456 // TestHorizontalAlgebraLayouts(mat, NULL, false);
457}
458
459}// end of namespace
460
461#endif //__H__UG_matrix_overlap_impl
specializations are responsible to pack and unpack interface data during communication.
Definition pcl_communication_structs.h:790
You may add elements to this interface and iterate over them.
Definition pcl_communication_structs.h:207
iterator push_back(const Element &elem)
Definition pcl_communication_structs.h:245
ElemContainer::const_iterator const_iterator
Definition pcl_communication_structs.h:237
iterator end()
Definition pcl_communication_structs.h:293
iterator begin()
Definition pcl_communication_structs.h:292
Element & get_element(iterator iter)
Definition pcl_communication_structs.h:298
int get_target_proc() const
Definition pcl_communication_structs.h:309
size_t size() const
returns the number of elements that are stored in the interface.
Definition pcl_communication_structs.h:306
Extends the HorizontalAlgebraLayouts by vertical layouts.
Definition algebra_layouts.h:121
A Buffer for binary data.
Definition binary_buffer.h:56
size_t read_pos() const
returns the current read-pos (in bytes)
Definition binary_buffer_impl.h:48
char * buffer()
returns the raw buffer pointer or NULL if the buffer is empty (capacity() == 0)
Definition binary_buffer_impl.h:94
size_t write_pos() const
returns the current write-pos (in bytes)
Definition binary_buffer_impl.h:53
Definition matrix_overlap_impl.h:56
int get_required_buffer_size(const Interface &interface)
returns the size of the buffer in bytes, that will be required for interface-communication.
Definition matrix_overlap_impl.h:60
ComPol_MatCopyDiag(TMatrix &mat)
Definition matrix_overlap_impl.h:58
TMatrix & m_mat
Definition matrix_overlap_impl.h:95
bool collect(ug::BinaryBuffer &buff, const Interface &interface)
should write data which is associated with the interface elements to the buffer.
Definition matrix_overlap_impl.h:68
bool extract(ug::BinaryBuffer &buff, const Interface &interface)
extract data from the buffer and assigns it to the interface-elements.
Definition matrix_overlap_impl.h:81
Communication Policy to copy couplings between interfaces.
Definition communication_policies.h:1273
Highly specialized communication policy for matrix overlap creation.
Definition matrix_overlap_impl.h:114
TMatrix & m_mat
Definition matrix_overlap_impl.h:393
virtual void post_process()
After communication is done, this method should be called to create the overlap.
Definition matrix_overlap_impl.h:219
TMatrix::value_type block_t
Definition matrix_overlap_impl.h:116
AlgebraIDVec & m_globalIDs
map localID->globalID
Definition matrix_overlap_impl.h:396
std::set< AlgebraID > m_recvNewIDs
Here we just collect the new IDs which are added to the local matrix.
Definition matrix_overlap_impl.h:428
AlgebraIDHashList m_algIDHash
map globalID->localID
Definition matrix_overlap_impl.h:399
std::map< ExtCon, block_t > m_recvNewCons
New connections received from other processes.
Definition matrix_overlap_impl.h:426
ComPol_MatCreateOverlap(TMatrix &rMat, AlgebraIDVec &vGlobalID)
Constructor setting the vector.
Definition matrix_overlap_impl.h:121
virtual bool collect(ug::BinaryBuffer &buff, const Interface &interface)
writes the interface values into a buffer that will be sent
Definition matrix_overlap_impl.h:130
virtual bool extract(ug::BinaryBuffer &buff, const Interface &interface)
writes values from a buffer into the interface
Definition matrix_overlap_impl.h:171
void insert(const key_t &key, const value_t &val)
Definition hash_impl.hpp:199
value_t & get_entry(const key_t &key)
Definition hash_impl.hpp:158
void GenerateGlobalAlgebraIDs(pcl::InterfaceCommunicator< TLayout > &communicator, std::vector< AlgebraID > &idsOut, size_t numIDs, const TLayout &masterLayout, const TLayout &slaveLayout)
Generates a set of unique global algebra ids.
Definition parallelization_util.h:81
void SetDirichletRow(TSparseMatrix &A, size_t i, size_t alpha)
Definition sparsematrix_util.h:796
#define UG_THROW(msg)
Definition error.h:57
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
Definition smart_pointer.h:814
the ug namespace
void GetLayoutTargetProcs(std::vector< int > &procsOut, const TLayout &layout)
Definition matrix_overlap_impl.h:43
void CreateOverlap(TMatrix &matInOut)
Adds coefficients and connections to create an overlapping matrix.
Definition matrix_overlap_impl.h:434
void Deserialize(TIStream &buf, ParallelVector< T > &v)
Deerialize for ParallelVector<T>
Definition restart_bridge.cpp:112
std::vector< AlgebraID > AlgebraIDVec
Definition algebra_id.h:63
void GenerateAlgebraIDHashList(AlgebraIDHashList &hash, AlgebraIDVec &algebraIDs)
Creates a hash which allows a algebraID->localIndex mapping.
Definition algebra_id.h:70
void Serialize(TOStream &buf, const ParallelVector< T > &v)
Serialize for ParallelVector<T>
Definition restart_bridge.cpp:103
#define PROFILE_FUNC_GROUP(groups)
Definition profiler.h:258
#define PROFILE_BEGIN_GROUP(name, groups)
Definition profiler.h:255
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
this type is used to identify distributed objects.
Definition algebra_id.h:46
Definition matrix_overlap_impl.h:401
AlgebraID fromID
Definition matrix_overlap_impl.h:411
int conProc
Definition matrix_overlap_impl.h:413
size_t fromInd
Definition matrix_overlap_impl.h:410
ExtCon(size_t _fromInd, const AlgebraID &_fromID, const AlgebraID &_toID, int _conProc)
Definition matrix_overlap_impl.h:402
bool operator<(const ExtCon &ec) const
Definition matrix_overlap_impl.h:415
AlgebraID toID
Definition matrix_overlap_impl.h:412
Definition communication_policies.h:58