Line data Source code
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
45 : #include "lib_algebra/parallelization/parallelization.h"
46 : #include "pcl/pcl.h"
47 : #endif
48 :
49 : #include "common/log.h"
50 :
51 :
52 : namespace ug{
53 :
54 :
55 : //! todo: replace DebugID
56 : extern 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 :
67 : template <class TVec, size_t N>
68 : class SlicingData{
69 :
70 : public:
71 : typedef std::vector<int> slice_desc_set;
72 : typedef TVec slice_desc_type_vector;
73 : typedef typename TVec::value_type slice_desc_type;
74 :
75 : protected:
76 : bool m_valid;
77 : slice_desc_type_vector m_slice_types; //!< global vector with mappings iglobal -> type(iglobal)
78 : slice_desc_set m_slice_set[N]; //!< N mappings: islice(type) -> iglobal
79 :
80 :
81 : public:
82 : //! Constructor
83 0 : SlicingData() : m_valid(false)
84 : {}
85 :
86 : /*! Builds index mappings based on types */
87 : SlicingData(const slice_desc_type_vector &types)
88 : : m_valid(true), m_slice_types(types)
89 : {
90 : reset_set_mappings();
91 : }
92 :
93 : /// Copy types.
94 : void set_types(const slice_desc_type_vector &types, bool bClear=false)
95 : {
96 : UG_DLOG(SchurDebug, 5,"SlicingData::set_types:" << bClear);
97 0 : m_slice_types = types;
98 : reset_set_mappings();
99 0 : }
100 :
101 : bool is_valid()
102 : { return m_valid;}
103 :
104 :
105 : protected:
106 :
107 : //! Clear all sets.
108 0 : void clear_set_mappings()
109 : {
110 : const size_t ntypes = m_slice_types.size();
111 0 : for (size_t i=0; i< ntypes; ++i)
112 : {
113 : slice_desc_type tt = get_type(i);
114 : slice(tt).clear();
115 : }
116 0 : }
117 :
118 : //! Auto fill for sets.
119 : /*! Assigns every index i=0.. m_slice_types.size()-1 to exactly one set. */
120 0 : void fill_set_mappings()
121 : {
122 : UG_DLOG(SchurDebug, 5, "SlicingData::fill_set_mappings" << std::endl);
123 : const size_t ntypes = m_slice_types.size();
124 0 : for (size_t i=0; i< ntypes; ++i)
125 : {
126 : slice_desc_type tt = get_type(i);
127 : slice_desc_set &set= slice(tt);
128 : UG_DLOG(SchurDebug, 8, i << "->" << tt << std::endl);
129 0 : 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 0 : m_valid = true;
140 0 : }
141 :
142 : //! Clear & auto fill
143 : void reset_set_mappings()
144 : {
145 : UG_DLOG(SchurDebug, 5,"SlicingData::reset_set_mappings");
146 0 : clear_set_mappings();
147 0 : fill_set_mappings();
148 : }
149 : public:
150 :
151 :
152 :
153 :
154 : /// Copy: slice of vector -> small vector.
155 : template<class VT>
156 0 : 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 0 : small_dst.resize(slice_desc.size());
162 : slice_desc_set::const_iterator elem = slice_desc.begin();
163 0 : for (size_t i=0; i<slice_desc.size(); ++i, ++elem)
164 0 : { small_dst[i] = full_src[*elem]; }
165 0 : }
166 :
167 :
168 :
169 : /// Copy: small vector -> slice of a vector.
170 : template<class VT>
171 0 : 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 0 : 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 0 : full_dst[*elem] = small_src[i];
180 : }
181 0 : }
182 :
183 : /// Add: slice of vector -> small vector
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 :
194 : /// Add: small vector -> slice of a vector
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 :
205 : /// substract: slice of vector -> small vector
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 :
210 : /// substract: small vector -> slice of a vector
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 0 : 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 0 : Aslice.resize_and_clear(row_slice.size(), col_slice.size());
224 :
225 : int ii=0;
226 0 : for (slice_desc_set::const_iterator elem = row_slice.begin();
227 0 : elem!=row_slice.end(); ++elem, ++ii)
228 : {
229 0 : const int i = *elem; // global index
230 0 : for(typename MT::const_row_iterator it = A.begin_row(i);
231 0 : 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 0 : if (find_index(col_type, j, jj))
238 : {
239 0 : Aslice(ii, jj) = it.value();
240 : }
241 : }
242 : }
243 :
244 0 : }
245 :
246 : //! Number of elements for each type
247 : size_t get_num_elems(slice_desc_type type) const
248 : {return slice(type).size();}
249 :
250 :
251 : //! Create a (partial) clone of the vector, without copying values
252 : template<class VT>
253 0 : SmartPtr<VT> slice_clone_without_values(const VT &full_src, slice_desc_type type) const
254 : {
255 : const slice_desc_set &slice_desc = slice(type);
256 :
257 0 : 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 0 : return clone;
268 :
269 : }
270 :
271 : //! Sets an existing sliced vector up correctly
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 :
289 : //! Create a (partial) clone
290 : template<class VT>
291 0 : SmartPtr<VT> slice_clone(const VT &full_src, slice_desc_type type) const
292 : {
293 0 : SmartPtr<VT> clone = slice_clone_without_values(full_src, type);
294 0 : get_vector_slice(full_src, type, *clone);
295 0 : return clone;
296 : }
297 :
298 :
299 : protected:
300 :
301 : //! returns type for a global index
302 : slice_desc_type get_type(size_t index)
303 : {return m_slice_types[index];}
304 :
305 :
306 : //! returns the set of global indices for a given type
307 : const slice_desc_set &slice(slice_desc_type type) const
308 0 : {return m_slice_set[type];}
309 :
310 : slice_desc_set &slice(slice_desc_type type)
311 0 : {return m_slice_set[type];}
312 :
313 :
314 : //! returns: local index for a global index (if found) or undefined else.
315 0 : 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 0 : index = myset.size();
322 :
323 : //slice_desc_set::const_iterator it = myset.find(gindex);
324 0 : slice_desc_set::const_iterator it = lower_bound(myset.begin(), myset.end(), gindex);
325 0 : if (it != myset.end() && *it == gindex) {
326 : //index =// *it;
327 0 : 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 0 : return found;
334 : }
335 : #ifdef UG_PARALLEL
336 : public:
337 : //! Create new slice layouts (as a subset from full layouts).
338 : SmartPtr<AlgebraLayouts> create_slice_layouts(ConstSmartPtr<AlgebraLayouts> fullLayouts, slice_desc_type type) const
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 :
355 : protected:
356 : void replace_indices_in_layout(slice_desc_type type, IndexLayout &il) const
357 : {
358 :
359 : UG_DLOG(SchurDebug, 5, "SlicingData::replace_indices_in_layout" << std::endl);
360 : IndexLayout::iterator iter;
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 :
367 : IndexLayout::Interface::iterator eiter;
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_ */
|
