ug4
Loading...
Searching...
No Matches
slicing.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 * Author: Arne Nägel
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__LIB_DISC__OPERATOR__LINEAR_OPERATOR__SCHUR_SLICING_H_
34#define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__SCHUR_SLICING_H_
35
36
37
38
39#include <iostream>
40#include <sstream>
41#include <string>
42#include <set>
43
44#ifdef UG_PARALLEL
46#include "pcl/pcl.h"
47#endif
48
49#include "common/log.h"
50
51
52namespace ug{
53
54
56extern DebugID SchurDebug;
57
58
59/*
60 *
61 * Allows splitting index sets for vectors/operator into different slices,
62 * e.g., for schur complement, component-wise splitting, etc.
63 *
64 *
65 * **/
66
67template <class TVec, size_t N>
69
70public:
71 typedef std::vector<int> slice_desc_set;
73 typedef typename TVec::value_type slice_desc_type;
74
75protected:
76 bool m_valid;
79
80
81public:
84 {}
85
88 : m_valid(true), m_slice_types(types)
89 {
91 }
92
94 void set_types(const slice_desc_type_vector &types, bool bClear=false)
95 {
96 UG_DLOG(SchurDebug, 5,"SlicingData::set_types:" << bClear);
97 m_slice_types = types;
99 }
100
101 bool is_valid()
102 { return m_valid;}
103
104
105protected:
106
109 {
110 const size_t ntypes = m_slice_types.size();
111 for (size_t i=0; i< ntypes; ++i)
112 {
114 slice(tt).clear();
115 }
116 }
117
119
121 {
122 UG_DLOG(SchurDebug, 5, "SlicingData::fill_set_mappings" << std::endl);
123 const size_t ntypes = m_slice_types.size();
124 for (size_t i=0; i< ntypes; ++i)
125 {
127 slice_desc_set &set= slice(tt);
128 UG_DLOG(SchurDebug, 8, i << "->" << tt << std::endl);
129 set.push_back(i);
130 }
131
132
133 UG_DLOG(SchurDebug, 5, "SlicingData::fill_set_mappings:" << ntypes);
134 for (size_t tt=0; tt<N; ++tt)
135 UG_DLOG(SchurDebug, 6, " " << m_slice_set[tt].size());
136
137 UG_DLOG(SchurDebug, 6, std::endl);
138
139 m_valid = true;
140 }
141
144 {
145 UG_DLOG(SchurDebug, 5,"SlicingData::reset_set_mappings");
148 }
149public:
150
151
152
153
155 template<class VT>
156 void get_vector_slice(const VT &full_src, slice_desc_type desc, VT &small_dst) const
157 {
158 const slice_desc_set &slice_desc = slice(desc);
159
160 // UG_DLOG("get_vector_slice:" << slice_desc.size());
161 small_dst.resize(slice_desc.size());
162 slice_desc_set::const_iterator elem = slice_desc.begin();
163 for (size_t i=0; i<slice_desc.size(); ++i, ++elem)
164 { small_dst[i] = full_src[*elem]; }
165 }
166
167
168
170 template<class VT>
171 void set_vector_slice(const VT &small_src, VT &full_dst, slice_desc_type desc) const
172 {
173 const slice_desc_set &slice_desc = slice(desc);
174 //UG_LOG("get_vector_slice:" << slice_desc.size());
175 slice_desc_set::const_iterator elem = slice_desc.begin();
176 for (size_t i=0; i<slice_desc.size(); ++i, ++elem)
177 {
178 // UG_DLOG(SchurDebug, 8, "Copy: "<< i << "->" << *elem << "v=" << small_src[i] << std::endl)
179 full_dst[*elem] = small_src[i];
180 }
181 }
182
184 template<class VT>
185 void add_vector_slice(const VT &full_src, slice_desc_type desc, VT &small_dst, double sigma=1.0) const
186 {
187 const slice_desc_set &slice_desc = slice(desc);
188 small_dst.resize(slice_desc.size());
189 slice_desc_set::const_iterator elem = slice_desc.begin();
190 for (size_t i=0; i<slice_desc.size(); ++i, ++elem)
191 { small_dst[i] += sigma*full_src[*elem]; }
192 }
193
195 template<class VT>
196 void add_vector_slice(const VT &small_src, VT &full_dst, slice_desc_type desc, double sigma=1.0) const
197 {
198 const slice_desc_set &slice_desc = slice(desc);
199 slice_desc_set::const_iterator elem = slice_desc.begin();
200 for (size_t i=0; i<slice_desc.size(); ++i, ++elem)
201 { full_dst[*elem] += sigma*small_src[i]; }
202 }
203
204
206 template<class VT>
207 void subtract_vector_slice(const VT &full_src, slice_desc_type desc, VT &small_dst) const
208 { add_vector_slice(full_src, desc, small_dst, -1.0); }
209
211 template<class VT>
212 void subtract_vector_slice(const VT &small_src, VT &full_dst, slice_desc_type desc) const
213 { add_vector_slice(small_src, full_dst, desc, -1.0); }
214
215
216 // Extracts a slice from a (full) matrix
217 template<class MT>
218 void get_matrix(const MT &A, slice_desc_type row_type, slice_desc_type col_type, MT &Aslice) const
219 {
220 const slice_desc_set &row_slice = slice(row_type);
221 const slice_desc_set &col_slice = slice(col_type);
222 UG_DLOG(SchurDebug, 5, "SlicingData::get_matrix:" << row_slice.size() << "x" << col_slice.size()<< std::endl)
223 Aslice.resize_and_clear(row_slice.size(), col_slice.size());
224
225 int ii=0;
226 for (slice_desc_set::const_iterator elem = row_slice.begin();
227 elem!=row_slice.end(); ++elem, ++ii)
228 {
229 const int i = *elem; // global index
230 for(typename MT::const_row_iterator it = A.begin_row(i);
231 it != A.end_row(i); ++it)
232
233 {
234 const int j=it.index();
235 int jj;
236 // if (get_type(j)!=col_type) continue;
237 if (find_index(col_type, j, jj))
238 {
239 Aslice(ii, jj) = it.value();
240 }
241 }
242 }
243
244 }
245
248 {return slice(type).size();}
249
250
252 template<class VT>
254 {
255 const slice_desc_set &slice_desc = slice(type);
256
257 SmartPtr<VT> clone(new VT(slice_desc.size()));
258#ifdef UG_PARALLEL
259 // set layout
260 clone->set_layouts(create_slice_layouts(full_src.layouts(), type));
261
262 // set mask
263 uint mask = full_src.get_storage_mask();
264 clone->set_storage_type(mask);
265 //std::cout << "Storage Type:" << full_src.get_storage_type() <<", mask=" << mask << std::endl;
266#endif
267 return clone;
268
269 }
270
272 template<class VT>
273 void setup_slice_like(const VT &full_src, slice_desc_type type, VT & vectorslice) const
274 {
275 const slice_desc_set &slice_desc = slice(type);
276
277 vectorslice.resize(slice_desc.size());
278#ifdef UG_PARALLEL
279 // set layout
280 vectorslice.set_layouts(create_slice_layouts(full_src.layouts(), type));
281
282 // set mask
283 uint mask = full_src.get_storage_mask();
284 vectorslice->set_storage_type(mask);
285 //std::cout << "Storage Type:" << full_src.get_storage_type() <<", mask=" << mask << std::endl;
286#endif
287 }
288
290 template<class VT>
291 SmartPtr<VT> slice_clone(const VT &full_src, slice_desc_type type) const
292 {
293 SmartPtr<VT> clone = slice_clone_without_values(full_src, type);
294 get_vector_slice(full_src, type, *clone);
295 return clone;
296 }
297
298
299protected:
300
303 {return m_slice_types[index];}
304
305
308 {return m_slice_set[type];}
309
312
313
315 bool find_index(slice_desc_type type, int gindex, int &index) const
316 {
317 // WARNING int index < size_t myset.size() WARNING
318 bool found=false;
319
320 const slice_desc_set &myset=slice(type);
321 index = myset.size();
322
323 //slice_desc_set::const_iterator it = myset.find(gindex);
324 slice_desc_set::const_iterator it = lower_bound(myset.begin(), myset.end(), gindex);
325 if (it != myset.end() && *it == gindex) {
326 //index =// *it;
327 index = std::distance(myset.begin(), it);
328 found = true;
329 }
330 // if (found && index >=myset.size()) {
331 UG_ASSERT( (!found || index<(int)myset.size()) , "Invalid index found!");
332 // }
333 return found;
334 }
335#ifdef UG_PARALLEL
336public:
339 {
340 UG_DLOG(SchurDebug, 4, "SlicingData::create_slice_layouts" << std::endl);
341 // Copy layouts.
342 SmartPtr<AlgebraLayouts> slice_layouts = make_sp(new AlgebraLayouts(*fullLayouts));
343
344 // Convert layouts (vector->slice)
345 replace_indices_in_layout(type, slice_layouts->master());
346 replace_indices_in_layout(type, slice_layouts->slave());
347
348 UG_DLOG(SchurDebug, 3, "BEFORE:")
349 UG_DLOG(SchurDebug, 3, *fullLayouts);
350 UG_DLOG(SchurDebug, 3, "AFTER:")
351 UG_DLOG(SchurDebug, 3, *slice_layouts);
352 return slice_layouts;
353 }
354
355protected:
357 {
358
359 UG_DLOG(SchurDebug, 5, "SlicingData::replace_indices_in_layout" << std::endl);
361 for (iter = il.begin(); iter!=il.end(); ++iter)
362 {
363 // iterate over interfaces
364 // i.e., pcl::OrderedInterface<size_t, std::vector>
365 IndexLayout::Interface &interf=il.interface(iter);
366
368 for (eiter = interf.begin(); eiter!=interf.end(); ++eiter)
369 {
370 // replace elements in interface
371 size_t myindex = interf.get_element(eiter);
372 int newindex;
373 bool found=find_index(type, myindex, newindex);
374 // UG_COND_THROW(!found, "SlicingData:: Did not find "<< myindex << "in typelist for type="<< type << std::endl);
375
376 // Replace or remove from IF
377 if (found) {
378 UG_DLOG(SchurDebug, 7, "Replacing:" << myindex << "->" <<newindex << std::endl);
379 interf.get_element(eiter) = newindex;
380 }
381 else {
382 UG_DLOG(SchurDebug, 7, "Deleting:" << myindex << std::endl);
383 interf.erase(eiter);
384 eiter--;
385 }
386
387 }
388 }
389 }
390
391#endif /* UG_PARALLEL */
392
393};
394
395}
396
397
398#endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__SCHUR_SLICING_H_ */
Definition smart_pointer.h:296
Definition smart_pointer.h:108
You may add elements to this interface and iterate over them.
Definition pcl_communication_structs.h:207
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
iterator erase(iterator iter)
Definition pcl_communication_structs.h:286
iterator end(size_t level=0)
returns the iterator to the last interface of the layout.
Definition pcl_communication_structs.h:492
Interface & interface(iterator iter)
returns the interface to the given iterator.
Definition pcl_communication_structs.h:505
iterator begin(size_t level=0)
returns the iterator to the first interface of the layout.
Definition pcl_communication_structs.h:486
InterfaceMap::iterator iterator
An iterator that allows to iterate over the interfaces stored in the layout.
Definition pcl_communication_structs.h:476
Extends the HorizontalAlgebraLayouts by vertical layouts.
Definition algebra_layouts.h:121
Definition slicing.h:68
slice_desc_set & slice(slice_desc_type type)
Definition slicing.h:310
void fill_set_mappings()
Auto fill for sets.
Definition slicing.h:120
void subtract_vector_slice(const VT &small_src, VT &full_dst, slice_desc_type desc) const
substract: small vector -> slice of a vector
Definition slicing.h:212
slice_desc_type get_type(size_t index)
returns type for a global index
Definition slicing.h:302
void add_vector_slice(const VT &small_src, VT &full_dst, slice_desc_type desc, double sigma=1.0) const
Add: small vector -> slice of a vector.
Definition slicing.h:196
SlicingData()
Constructor.
Definition slicing.h:83
SlicingData(const slice_desc_type_vector &types)
Definition slicing.h:87
TVec slice_desc_type_vector
Definition slicing.h:72
void get_vector_slice(const VT &full_src, slice_desc_type desc, VT &small_dst) const
Copy: slice of vector -> small vector.
Definition slicing.h:156
void replace_indices_in_layout(slice_desc_type type, IndexLayout &il) const
Definition slicing.h:356
std::vector< int > slice_desc_set
Definition slicing.h:71
bool m_valid
Definition slicing.h:76
SmartPtr< AlgebraLayouts > create_slice_layouts(ConstSmartPtr< AlgebraLayouts > fullLayouts, slice_desc_type type) const
Create new slice layouts (as a subset from full layouts).
Definition slicing.h:338
size_t get_num_elems(slice_desc_type type) const
Number of elements for each type.
Definition slicing.h:247
void add_vector_slice(const VT &full_src, slice_desc_type desc, VT &small_dst, double sigma=1.0) const
Add: slice of vector -> small vector.
Definition slicing.h:185
void reset_set_mappings()
Clear & auto fill.
Definition slicing.h:143
bool is_valid()
Definition slicing.h:101
TVec::value_type slice_desc_type
Definition slicing.h:73
const slice_desc_set & slice(slice_desc_type type) const
returns the set of global indices for a given type
Definition slicing.h:307
void set_vector_slice(const VT &small_src, VT &full_dst, slice_desc_type desc) const
Copy: small vector -> slice of a vector.
Definition slicing.h:171
slice_desc_type_vector m_slice_types
global vector with mappings iglobal -> type(iglobal)
Definition slicing.h:77
void get_matrix(const MT &A, slice_desc_type row_type, slice_desc_type col_type, MT &Aslice) const
Definition slicing.h:218
void setup_slice_like(const VT &full_src, slice_desc_type type, VT &vectorslice) const
Sets an existing sliced vector up correctly.
Definition slicing.h:273
slice_desc_set m_slice_set[N]
N mappings: islice(type) -> iglobal.
Definition slicing.h:78
SmartPtr< VT > slice_clone_without_values(const VT &full_src, slice_desc_type type) const
Create a (partial) clone of the vector, without copying values.
Definition slicing.h:253
void subtract_vector_slice(const VT &full_src, slice_desc_type desc, VT &small_dst) const
substract: slice of vector -> small vector
Definition slicing.h:207
void clear_set_mappings()
Clear all sets.
Definition slicing.h:108
bool find_index(slice_desc_type type, int gindex, int &index) const
returns: local index for a global index (if found) or undefined else.
Definition slicing.h:315
void set_types(const slice_desc_type_vector &types, bool bClear=false)
Copy types.
Definition slicing.h:94
SmartPtr< VT > slice_clone(const VT &full_src, slice_desc_type type) const
Create a (partial) clone.
Definition slicing.h:291
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_DLOG(__debugID__, level, msg)
Definition log.h:298
unsigned int uint
Definition types.h:114
the ug namespace
DebugID SchurDebug
todo: replace DebugID
Definition schur.h:55
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836