Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_OverlappingPartitioner_def.hpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#ifndef IFPACK2_OVERLAPPINGPARTITIONER_DEF_HPP
44#define IFPACK2_OVERLAPPINGPARTITIONER_DEF_HPP
45#include <vector>
46#include <string>
47#include <algorithm>
48#include "Ifpack2_ConfigDefs.hpp"
49#include "Ifpack2_OverlappingPartitioner_decl.hpp"
50#include "Teuchos_Array.hpp"
51#include "Teuchos_ArrayRCP.hpp"
52
53namespace Ifpack2 {
54
55template<class GraphType>
57OverlappingPartitioner (const Teuchos::RCP<const row_graph_type>& graph) :
58 NumLocalParts_ (1),
59 Graph_ (graph),
60 OverlappingLevel_ (0),
61 IsComputed_ (false),
62 verbose_ (false),
63 maintainSparsity_(false)
64{}
65
66
67template<class GraphType>
69
70
71template<class GraphType>
72int
74{
75 return NumLocalParts_;
76}
77
78
79template<class GraphType>
81{
82 return OverlappingLevel_;
83}
84
85
86template<class GraphType>
87typename GraphType::local_ordinal_type
89operator () (const local_ordinal_type MyRow) const
90{
91 TEUCHOS_TEST_FOR_EXCEPTION(
92 MyRow < 0 || Teuchos::as<size_t> (MyRow) > Graph_->getLocalNumRows (),
93 std::runtime_error,
94 "Ifpack2::OverlappingPartitioner::operator(): "
95 "Invalid local row index " << MyRow << ".");
96
97 return Partition_[MyRow];
98}
99
100
101//==============================================================================
102template<class GraphType>
103typename GraphType::local_ordinal_type
105operator() (const local_ordinal_type i, const local_ordinal_type j) const
106{
107 TEUCHOS_TEST_FOR_EXCEPTION(
108 i < 0 || i > Teuchos::as<local_ordinal_type> (NumLocalParts_),
109 std::runtime_error,
110 "Ifpack2::OverlappingPartitioner::operator(): "
111 "Invalid local row index i=" << i << ".");
112 TEUCHOS_TEST_FOR_EXCEPTION(
113 j < 0 || j > Teuchos::as<local_ordinal_type> (Parts_[i].size ()),
114 std::runtime_error,
115 "Ifpack2::OverlappingPartitioner::operator(): "
116 "Invalid node index j=" << j << ".");
117 return Parts_[i][j];
118}
119
120//==============================================================================
121template<class GraphType>
122size_t
124numRowsInPart (const local_ordinal_type Part) const
125{
126 TEUCHOS_TEST_FOR_EXCEPTION(
127 Part < 0 || Part > Teuchos::as<local_ordinal_type> (NumLocalParts_),
128 std::runtime_error,
129 "Ifpack2::OverlappingPartitioner::numRowsInPart: "
130 "Invalid partition index Part=" << Part << ".");
131 return Parts_[Part].size ();
132}
133
134//==============================================================================
135template<class GraphType>
136void
138rowsInPart (const local_ordinal_type Part,
139 Teuchos::ArrayRCP<local_ordinal_type>& List) const
140{
141 // Let numRowsInPart do the sanity checking...
142 const size_t numRows = numRowsInPart (Part);
143 for (size_t i = 0; i < numRows; ++i) {
144 List[i] = Parts_[Part][i];
145 }
146}
147
148//==============================================================================
149template<class GraphType>
150Teuchos::ArrayView<const typename GraphType::local_ordinal_type>
152{
153 return Partition_.view (0, Graph_->getLocalNumRows ());
154}
155
156//==============================================================================
157template<class GraphType>
158void
160setParameters (Teuchos::ParameterList& List)
161{
162 NumLocalParts_ = List.get("partitioner: local parts", NumLocalParts_);
163 OverlappingLevel_ = List.get("partitioner: overlap", OverlappingLevel_);
164 verbose_ = List.get("partitioner: print level", verbose_);
165 maintainSparsity_ = List.get("partitioner: maintain sparsity", false);
166 typedef Teuchos::RCP< Tpetra::Map<typename GraphType::local_ordinal_type, typename GraphType::global_ordinal_type, typename GraphType::node_type> const > map_type;
167 typedef Teuchos::RCP<Tpetra::Import<typename GraphType::local_ordinal_type, typename GraphType::global_ordinal_type, typename GraphType::node_type> const > import_type;
168
169 // when using overlapping schwarz wth combineMode ADD, we need import and
170 // overlap row map to adjust weights for ith dof. Specifically, we sum
171 // all blocks (including off processor ones) that contain i.
172 import_type theImport = Graph_->getImporter();
173 if (theImport != Teuchos::null) List.set< import_type >("theImport",theImport);
174 List.set< map_type >("OverlapRowMap",Graph_->getRowMap());
175
176 if (NumLocalParts_ < 0) {
177 NumLocalParts_ = Graph_->getLocalNumRows() / (-NumLocalParts_);
178 }
179 if (NumLocalParts_ == 0) {
180 NumLocalParts_ = 1;
181 }
182
183 // Sanity checking
184 TEUCHOS_TEST_FOR_EXCEPTION(
185 NumLocalParts_ < 0 ||
186 Teuchos::as<size_t> (NumLocalParts_) > Graph_->getLocalNumRows(),
187 std::runtime_error,
188 "Ifpack2::OverlappingPartitioner::setParameters: "
189 "Invalid NumLocalParts_ = " << NumLocalParts_ << ".");
190 TEUCHOS_TEST_FOR_EXCEPTION(
191 OverlappingLevel_ < 0, std::runtime_error,
192 "Ifpack2::OverlappingPartitioner::setParameters: "
193 "Invalid OverlappingLevel_ = " << OverlappingLevel_ << ".");
194
195 setPartitionParameters(List);
196}
197
198//==============================================================================
199template<class GraphType>
201{
202 using std::cout;
203 using std::endl;
204
205 TEUCHOS_TEST_FOR_EXCEPTION(
206 NumLocalParts_ < 1 || OverlappingLevel_ < 0,
207 std::runtime_error,
208 "Ifpack2::OverlappingPartitioner::compute: "
209 "Invalid NumLocalParts_ or OverlappingLevel_.");
210
211 // std::string's constructor has some overhead, so it's better to
212 // use const char[] for local constant strings.
213 const char printMsg[] = "OverlappingPartitioner: ";
214
215 if (verbose_ && (Graph_->getComm()->getRank() == 0)) {
216 cout << printMsg << "Number of local parts = "
217 << NumLocalParts_ << endl;
218 cout << printMsg << "Approx. Number of global parts = "
219 << NumLocalParts_ * Graph_->getComm ()->getSize () << endl;
220 cout << printMsg << "Amount of overlap = "
221 << OverlappingLevel_ << endl;
222 }
223
224 // 1.- allocate memory
225 Partition_.resize (Graph_->getLocalNumRows ());
226 //Parts_ is allocated in computeOverlappingPartitions_, where it is used
227
228 // 2.- sanity checks on input graph
229 TEUCHOS_TEST_FOR_EXCEPTION(
230 ! Graph_->isFillComplete (), std::runtime_error,
231 "Ifpack2::OverlappingPartitioner::compute: "
232 "The input graph must be fill complete.");
233
234 TEUCHOS_TEST_FOR_EXCEPTION(
235 Graph_->getGlobalNumRows () != Graph_->getGlobalNumCols (),
236 std::runtime_error,
237 "Ifpack2::OverlappingPartitioner::compute: "
238 "The input graph must be (globally) square.");
239
240 // 3.- perform non-overlapping partition
241 computePartitions ();
242
243 // 4.- compute the partitions with overlapping
244 computeOverlappingPartitions ();
245
246 // 5.- mark as computed
247 IsComputed_ = true;
248}
249
250//==============================================================================
251template<class GraphType>
253{
254 //If user has explicitly specified parts, then Partition_ size has been set to 0.
255 //In this case, there is no need to compute Parts_.
256 if (Partition_.size() == 0)
257 return;
258
259 const local_ordinal_type invalid =
260 Teuchos::OrdinalTraits<local_ordinal_type>::invalid();
261
262 // Old FIXME from Ifpack: the first part of this function should be elsewhere
263
264 // start defining the subgraphs for no overlap
265
266 std::vector<size_t> sizes;
267 sizes.resize (NumLocalParts_);
268
269 // 1.- compute how many rows are in each subgraph
270 for (int i = 0; i < NumLocalParts_; ++i) {
271 sizes[i] = 0;
272 }
273
274 for (size_t i = 0; i < Graph_->getLocalNumRows (); ++i) {
275 TEUCHOS_TEST_FOR_EXCEPTION(
276 Partition_[i] >= NumLocalParts_, std::runtime_error,
277 "Ifpack2::OverlappingPartitioner::computeOverlappingPartitions: "
278 "Partition_[i] > NumLocalParts_.");
279 // invalid indicates that this unknown is not in a nonoverlapping
280 // partition
281 if (Partition_[i] != invalid) {
282 sizes[Partition_[i]]++;
283 }
284 }
285
286 // 2.- allocate space for each subgraph
287 Parts_.resize (NumLocalParts_);
288 for (int i = 0; i < NumLocalParts_; ++i) {
289 Parts_[i].resize (sizes[i]);
290 }
291
292 // 3.- cycle over all rows and populate the vectors
293 for (int i = 0; i < NumLocalParts_; ++i) {
294 sizes[i] = 0;
295 }
296
297 for (size_t i = 0; i < Graph_->getLocalNumRows (); ++i) {
298 const local_ordinal_type part = Partition_[i];
299 if (part != invalid) {
300 const size_t count = sizes[part];
301 Parts_[part][count] = i;
302 sizes[part]++;
303 }
304 }
305
306 // If there is no overlap, we're done, so return
307 if (OverlappingLevel_ == 0) {
308 return;
309 }
310
311 // wider overlap requires further computations
312 for (int level = 1; level <= OverlappingLevel_; ++level) {
313 std::vector<std::vector<size_t> > tmp;
314 tmp.resize (NumLocalParts_);
315
316 // cycle over all rows in the local graph (that is the overlapping
317 // graph). For each row, all columns will belong to the subgraph
318 // of row `i'.
319
320 int MaxNumEntries_tmp = Graph_->getLocalMaxNumRowEntries();
321 nonconst_local_inds_host_view_type Indices("Indices",MaxNumEntries_tmp);
322 nonconst_local_inds_host_view_type newIndices("newIndices",MaxNumEntries_tmp);
323
324 if (!maintainSparsity_) {
325
326 local_ordinal_type numLocalRows = Graph_->getLocalNumRows();
327 for (int part = 0; part < NumLocalParts_ ; ++part) {
328 for (size_t i = 0; i < Teuchos::as<size_t> (Parts_[part].size ()); ++i) {
329 const local_ordinal_type LRID = Parts_[part][i];
330
331 size_t numIndices;
332 Graph_->getLocalRowCopy (LRID, Indices, numIndices);
333
334 for (size_t j = 0; j < numIndices; ++j) {
335 // use *local* indices only
336 const local_ordinal_type col = Indices[j];
337 if (col >= numLocalRows) {
338 continue;
339 }
340
341 // has this column already been inserted?
342 std::vector<size_t>::iterator where =
343 std::find (tmp[part].begin (), tmp[part].end (), Teuchos::as<size_t> (col));
344
345
346 if (where == tmp[part].end()) {
347 tmp[part].push_back (col);
348 }
349 }
350
351 //10-Jan-2017 JHU : A possible optimization is to avoid the std::find's in the loop above,
352 //but instead sort and make unique afterwards. One would have to be careful to only sort/
353 //make unique the insertions associated with the current row LRID, as for each row, LRID
354 //may occur after all other col indices for row LRID (see comment about zero pivot below).
355
356 // has this column already been inserted?
357 std::vector<size_t>::iterator where =
358 std::find (tmp[part].begin (), tmp[part].end (), Teuchos::as<size_t> (LRID));
359
360 // This happens here b/c Vanka on Stokes with Stabilized elements will have
361 // a zero pivot entry if this gets pushed back first. So... Last.
362 if (where == tmp[part].end ()) {
363 tmp[part].push_back (LRID);
364 }
365 }
366 } //for (int part = 0; ...
367
368 } else {
369 //maintainSparsity_ == true
370
371 for (int part = 0; part < NumLocalParts_ ; ++part) {
372 for (size_t i = 0; i < Teuchos::as<size_t> (Parts_[part].size ()); ++i) {
373 const local_ordinal_type LRID = Parts_[part][i];
374
375 size_t numIndices;
376 Graph_->getLocalRowCopy (LRID, Indices, numIndices);
377 //JJH: the entries in Indices are already sorted. However, the Tpetra documentation states
378 // that we can't count on this always being true, hence we sort. Also note that there are
379 // unused entries at the end of Indices (it's sized to hold any row). This means we can't
380 // just use Indices.end() in sorting and in std::includes
381 Tpetra::sort(Indices,numIndices);
382
383 for (size_t j = 0; j < numIndices; ++j) {
384 // use *local* indices only
385 const local_ordinal_type col = Indices[j];
386 if (Teuchos::as<size_t> (col) >= Graph_->getLocalNumRows ()) {
387 continue;
388 }
389
390 // has this column already been inserted?
391 std::vector<size_t>::iterator where =
392 std::find (tmp[part].begin (), tmp[part].end (), Teuchos::as<size_t> (col));
393
394
395 if (where == tmp[part].end()) {
396 // Check if row associated with "col" increases connectivity already defined by row LRID's stencil.
397 // If it does and maintainSparsity_ is true, do not add "col" to the current partition (block).
398 size_t numNewIndices;
399 Graph_->getLocalRowCopy(col, newIndices, numNewIndices);
400 Tpetra::sort(newIndices,numNewIndices);
401 auto Indices_rcp = Kokkos::Compat::persistingView<nonconst_local_inds_host_view_type>(Indices, 0, numIndices);
402 auto newIndices_rcp = Kokkos::Compat::persistingView<nonconst_local_inds_host_view_type>(newIndices, 0, numNewIndices);
403 bool isSubset = std::includes(Indices_rcp.begin(),Indices_rcp.begin()+numIndices,
404 newIndices_rcp.begin(),newIndices_rcp.begin()+numNewIndices);
405 if (isSubset) {
406 tmp[part].push_back (col);
407 }
408 }
409 }
410
411 // has this column already been inserted?
412 std::vector<size_t>::iterator where =
413 std::find (tmp[part].begin (), tmp[part].end (), Teuchos::as<size_t> (LRID));
414
415 // This happens here b/c Vanka on Stokes with Stabilized elements will have
416 // a zero pivot entry if this gets pushed back first. So... Last.
417 if (where == tmp[part].end ()) {
418 tmp[part].push_back (LRID);
419 }
420 }
421 } //for (int part = 0;
422
423 } //if (maintainSparsity_) ... else
424
425 // now I convert the STL vectors into Teuchos Array RCP's
426 //
427 // FIXME (mfh 12 July 2013) You could have started with ArrayRCP
428 // in the first place (which implements push_back and iterators)
429 // and avoided the copy.
430 for (int i = 0; i < NumLocalParts_; ++i) {
431 Parts_[i].resize (tmp[i].size ());
432 for (size_t j = 0; j < tmp[i].size (); ++j) {
433 Parts_[i][j] = tmp[i][j];
434 }
435 }
436 }
437}
438
439//==============================================================================
440template<class GraphType>
442{
443 return IsComputed_;
444}
445
446//==============================================================================
447template<class GraphType>
448std::ostream&
450{
451 Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
452 fos.setOutputToRootOnly (0);
453 describe (fos);
454 return os;
455}
456
457//==============================================================================
458template<class GraphType>
460{
461 std::ostringstream oss;
462 oss << Teuchos::Describable::description();
463 if (isComputed()) {
464 oss << "{status = computed";
465 }
466 else {
467 oss << "{status = is not computed";
468 }
469 oss <<"}";
470 return oss.str();
471}
472
473//==============================================================================
474template<class GraphType>
475void OverlappingPartitioner<GraphType>::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const
476{
477 using std::endl;
478 if (verbLevel == Teuchos::VERB_NONE) {
479 return;
480 }
481
482 os << "================================================================================" << endl;
483 os << "Ifpack2::OverlappingPartitioner" << endl;
484 os << "Number of local rows = " << Graph_->getLocalNumRows() << endl;
485 os << "Number of global rows = " << Graph_->getGlobalNumRows() << endl;
486 os << "Number of local parts = " << NumLocalParts_ << endl;
487 os << "Overlapping level = " << OverlappingLevel_ << endl;
488 os << "Is computed = " << IsComputed_ << endl;
489 os << "================================================================================" << endl;
490}
491
492}// namespace Ifpack2
493
494#define IFPACK2_OVERLAPPINGPARTITIONER_INSTANT(LO,GO,N) \
495 template class Ifpack2::OverlappingPartitioner<Tpetra::CrsGraph< LO, GO, N > >; \
496 template class Ifpack2::OverlappingPartitioner<Tpetra::RowGraph< LO, GO, N > >;
497
498#endif // IFPACK2_OVERLAPPINGPARTITIONER_DEF_HPP
size_t numRowsInPart(const local_ordinal_type Part) const
the number of rows contained in the given partition.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:124
OverlappingPartitioner(const Teuchos::RCP< const row_graph_type > &graph)
Constructor.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:57
virtual std::ostream & print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:449
virtual void compute()
Computes the partitions. Returns 0 if successful.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:200
virtual void setParameters(Teuchos::ParameterList &List)
Set all the parameters for the partitioner.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:160
virtual bool isComputed() const
Returns true if partitions have been computed successfully.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:441
virtual void computeOverlappingPartitions()
Computes the partitions. Returns 0 if successful.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:252
int numLocalParts() const
Number of computed local partitions.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:73
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:475
void rowsInPart(const local_ordinal_type Part, Teuchos::ArrayRCP< local_ordinal_type > &List) const
Fill List with the local indices of the rows in the (overlapping) partition Part.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:138
int overlappingLevel() const
The number of levels of overlap.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:80
virtual ~OverlappingPartitioner()
Destructor.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:68
virtual Teuchos::ArrayView< const local_ordinal_type > nonOverlappingPartition() const
A view of the local indices of the nonoverlapping partitions of each local row.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:151
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:459
local_ordinal_type operator()(const local_ordinal_type MyRow) const
Local index of the nonoverlapping partition of the given row.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:89
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74