Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_BlockCrsMatrix_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Tpetra: Templated Linear Algebra Services Package
4//
5// Copyright 2008 NTESS and the Tpetra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef TPETRA_BLOCKCRSMATRIX_DEF_HPP
11#define TPETRA_BLOCKCRSMATRIX_DEF_HPP
12
15
16#include "Tpetra_Details_Behavior.hpp"
17#include "Tpetra_Details_PackTraits.hpp"
18#include "Tpetra_Details_Profiling.hpp"
19#include "Tpetra_BlockMultiVector.hpp"
20#include "Tpetra_BlockView.hpp"
21
22#include "Teuchos_TimeMonitor.hpp"
23#ifdef HAVE_TPETRA_DEBUG
24#include <set>
25#endif // HAVE_TPETRA_DEBUG
26
27#include "KokkosSparse_spmv.hpp"
28
29namespace Tpetra {
30
31namespace Impl {
32
33template <typename T>
34struct BlockCrsRowStruct {
35 T totalNumEntries, totalNumBytes, maxRowLength;
36 KOKKOS_DEFAULTED_FUNCTION BlockCrsRowStruct() = default;
37 KOKKOS_DEFAULTED_FUNCTION BlockCrsRowStruct(const BlockCrsRowStruct& b) = default;
38 KOKKOS_INLINE_FUNCTION BlockCrsRowStruct(const T& numEnt, const T& numBytes, const T& rowLength)
39 : totalNumEntries(numEnt)
40 , totalNumBytes(numBytes)
41 , maxRowLength(rowLength) {}
42};
43
44template <typename T>
45static KOKKOS_INLINE_FUNCTION void operator+=(BlockCrsRowStruct<T>& a, const BlockCrsRowStruct<T>& b) {
46 a.totalNumEntries += b.totalNumEntries;
47 a.totalNumBytes += b.totalNumBytes;
48 a.maxRowLength = a.maxRowLength > b.maxRowLength ? a.maxRowLength : b.maxRowLength;
49}
50
51template <typename T, typename ExecSpace>
52struct BlockCrsReducer {
53 typedef BlockCrsReducer reducer;
54 typedef T value_type;
55 typedef Kokkos::View<value_type, ExecSpace, Kokkos::MemoryTraits<Kokkos::Unmanaged> > result_view_type;
56 value_type* value;
57
58 KOKKOS_INLINE_FUNCTION
59 BlockCrsReducer(value_type& val)
60 : value(&val) {}
61
62 KOKKOS_INLINE_FUNCTION void join(value_type& dst, value_type& src) const { dst += src; }
63 KOKKOS_INLINE_FUNCTION void join(value_type& dst, const value_type& src) const { dst += src; }
64 KOKKOS_INLINE_FUNCTION void init(value_type& val) const { val = value_type(); }
65 KOKKOS_INLINE_FUNCTION value_type& reference() { return *value; }
66 KOKKOS_INLINE_FUNCTION result_view_type view() const { return result_view_type(value); }
67};
68
69} // namespace Impl
70
71namespace { // (anonymous)
72
73// Implementation of BlockCrsMatrix::getLocalDiagCopy (non-deprecated
74// version that takes two Kokkos::View arguments).
75template <class Scalar, class LO, class GO, class Node>
76class GetLocalDiagCopy {
77 public:
78 typedef typename Node::device_type device_type;
79 typedef size_t diag_offset_type;
80 typedef Kokkos::View<const size_t*, device_type,
81 Kokkos::MemoryUnmanaged>
82 diag_offsets_type;
83 typedef typename ::Tpetra::CrsGraph<LO, GO, Node> global_graph_type;
84 typedef typename global_graph_type::local_graph_device_type local_graph_device_type;
85 typedef typename local_graph_device_type::row_map_type row_offsets_type;
86 typedef typename ::Tpetra::BlockMultiVector<Scalar, LO, GO, Node>::impl_scalar_type IST;
87 typedef Kokkos::View<IST***, device_type, Kokkos::MemoryUnmanaged> diag_type;
88 typedef Kokkos::View<const IST*, device_type, Kokkos::MemoryUnmanaged> values_type;
89
90 // Constructor
91 GetLocalDiagCopy(const diag_type& diag,
92 const values_type& val,
93 const diag_offsets_type& diagOffsets,
94 const row_offsets_type& ptr,
95 const LO blockSize)
96 : diag_(diag)
97 , diagOffsets_(diagOffsets)
98 , ptr_(ptr)
99 , blockSize_(blockSize)
100 , offsetPerBlock_(blockSize_ * blockSize_)
101 , val_(val) {}
102
103 KOKKOS_INLINE_FUNCTION void
104 operator()(const LO& lclRowInd) const {
105 using Kokkos::ALL;
106
107 // Get row offset
108 const size_t absOffset = ptr_[lclRowInd];
109
110 // Get offset relative to start of row
111 const size_t relOffset = diagOffsets_[lclRowInd];
112
113 // Get the total offset
114 const size_t pointOffset = (absOffset + relOffset) * offsetPerBlock_;
115
116 // Get a view of the block. BCRS currently uses LayoutRight
117 // regardless of the device.
118 typedef Kokkos::View<const IST**, Impl::BlockCrsMatrixLittleBlockArrayLayout,
119 device_type, Kokkos::MemoryTraits<Kokkos::Unmanaged> >
120 const_little_block_type;
121 const_little_block_type D_in(val_.data() + pointOffset,
122 blockSize_, blockSize_);
123 auto D_out = Kokkos::subview(diag_, lclRowInd, ALL(), ALL());
124 COPY(D_in, D_out);
125 }
126
127 private:
128 diag_type diag_;
129 diag_offsets_type diagOffsets_;
130 row_offsets_type ptr_;
131 LO blockSize_;
132 LO offsetPerBlock_;
133 values_type val_;
134};
135} // namespace
136
137template <class Scalar, class LO, class GO, class Node>
138std::ostream&
139BlockCrsMatrix<Scalar, LO, GO, Node>::
140 markLocalErrorAndGetStream() {
141 *(this->localError_) = true;
142 if ((*errs_).is_null()) {
143 *errs_ = Teuchos::rcp(new std::ostringstream());
144 }
145 return **errs_;
146}
147
148template <class Scalar, class LO, class GO, class Node>
149BlockCrsMatrix<Scalar, LO, GO, Node>::
150 BlockCrsMatrix()
151 : dist_object_type(Teuchos::rcp(new map_type()))
152 , // nonnull, so DistObject doesn't throw
153 graph_(Teuchos::rcp(new map_type()), 0)
154 , // FIXME (mfh 16 May 2014) no empty ctor yet
155 blockSize_(static_cast<LO>(0))
156 , X_colMap_(new Teuchos::RCP<BMV>())
157 , // ptr to a null ptr
158 Y_rowMap_(new Teuchos::RCP<BMV>())
159 , // ptr to a null ptr
160 pointImporter_(new Teuchos::RCP<typename crs_graph_type::import_type>())
161 , offsetPerBlock_(0)
162 , localError_(new bool(false))
163 , errs_(new Teuchos::RCP<std::ostringstream>()) // ptr to a null ptr
164{
165}
166
167template <class Scalar, class LO, class GO, class Node>
168BlockCrsMatrix<Scalar, LO, GO, Node>::
169 BlockCrsMatrix(const crs_graph_type& graph,
170 const LO blockSize)
171 : dist_object_type(graph.getMap())
172 , graph_(graph)
173 , rowMeshMap_(*(graph.getRowMap()))
174 , blockSize_(blockSize)
175 , X_colMap_(new Teuchos::RCP<BMV>())
176 , // ptr to a null ptr
177 Y_rowMap_(new Teuchos::RCP<BMV>())
178 , // ptr to a null ptr
179 pointImporter_(new Teuchos::RCP<typename crs_graph_type::import_type>())
180 , offsetPerBlock_(0)
181 , localError_(new bool(false))
182 , errs_(new Teuchos::RCP<std::ostringstream>()) // ptr to a null ptr
183{
185 TEUCHOS_TEST_FOR_EXCEPTION(
186 !graph_.isSorted(), std::invalid_argument,
187 "Tpetra::"
188 "BlockCrsMatrix constructor: The input CrsGraph does not have sorted "
189 "rows (isSorted() is false). This class assumes sorted rows.");
190
191 graphRCP_ = Teuchos::rcpFromRef(graph_);
192
193 // Trick to test whether LO is nonpositive, without a compiler
194 // warning in case LO is unsigned (which is generally a bad idea
195 // anyway). I don't promise that the trick works, but it
196 // generally does with gcc at least, in my experience.
197 const bool blockSizeIsNonpositive = (blockSize + 1 <= 1);
198 TEUCHOS_TEST_FOR_EXCEPTION(
199 blockSizeIsNonpositive, std::invalid_argument,
200 "Tpetra::"
201 "BlockCrsMatrix constructor: The input blockSize = "
202 << blockSize << " <= 0. The block size must be positive.");
203
204 const size_t blockSizeProd = static_cast<size_t>(blockSize) * static_cast<size_t>(blockSize);
205 TEUCHOS_TEST_FOR_EXCEPTION(
206 blockSizeProd > static_cast<size_t>(std::numeric_limits<LO>::max()),
207 std::overflow_error,
208 "Tpetra::BlockCrsMatrix constructor: blockSize * blockSize overflows LO");
209 offsetPerBlock_ = static_cast<LO>(blockSizeProd);
210
211 domainPointMap_ = BMV::makePointMap(*(graph.getDomainMap()), blockSize);
212 rangePointMap_ = BMV::makePointMap(*(graph.getRangeMap()), blockSize);
213
214 {
215 // These are rcp
216 const auto domainPointMap = getDomainMap();
217 const auto colPointMap = Teuchos::rcp(new typename BMV::map_type(BMV::makePointMap(*graph_.getColMap(), blockSize_)));
218 *pointImporter_ = Teuchos::rcp(new typename crs_graph_type::import_type(domainPointMap, colPointMap));
219 }
220 {
221 auto local_graph_h = graph.getLocalGraphHost();
222 auto ptr_h = local_graph_h.row_map;
223 ptrHost_ = decltype(ptrHost_)(Kokkos::ViewAllocateWithoutInitializing("graph row offset"), ptr_h.extent(0));
224 Kokkos::deep_copy(ptrHost_, ptr_h);
225
226 auto ind_h = local_graph_h.entries;
227 indHost_ = decltype(indHost_)(Kokkos::ViewAllocateWithoutInitializing("graph column indices"), ind_h.extent(0));
228 // DEEP_COPY REVIEW - HOST-TO-HOST
229 Kokkos::deep_copy(indHost_, ind_h);
230
231 const auto numValEnt = ind_h.extent(0) * offsetPerBlock();
232 val_ = decltype(val_)(impl_scalar_type_dualview("val", numValEnt));
233 }
234}
235
236template <class Scalar, class LO, class GO, class Node>
237BlockCrsMatrix<Scalar, LO, GO, Node>::
238 BlockCrsMatrix(const crs_graph_type& graph,
239 const typename local_matrix_device_type::values_type& values,
240 const LO blockSize)
241 : dist_object_type(graph.getMap())
242 , graph_(graph)
243 , rowMeshMap_(*(graph.getRowMap()))
244 , blockSize_(blockSize)
245 , X_colMap_(new Teuchos::RCP<BMV>())
246 , // ptr to a null ptr
247 Y_rowMap_(new Teuchos::RCP<BMV>())
248 , // ptr to a null ptr
249 pointImporter_(new Teuchos::RCP<typename crs_graph_type::import_type>())
250 , offsetPerBlock_(0)
251 , localError_(new bool(false))
252 , errs_(new Teuchos::RCP<std::ostringstream>()) // ptr to a null ptr
253{
255 TEUCHOS_TEST_FOR_EXCEPTION(
256 !graph_.isSorted(), std::invalid_argument,
257 "Tpetra::"
258 "BlockCrsMatrix constructor: The input CrsGraph does not have sorted "
259 "rows (isSorted() is false). This class assumes sorted rows.");
260
261 graphRCP_ = Teuchos::rcpFromRef(graph_);
262
263 // Trick to test whether LO is nonpositive, without a compiler
264 // warning in case LO is unsigned (which is generally a bad idea
265 // anyway). I don't promise that the trick works, but it
266 // generally does with gcc at least, in my experience.
267 const bool blockSizeIsNonpositive = (blockSize + 1 <= 1);
268 TEUCHOS_TEST_FOR_EXCEPTION(
269 blockSizeIsNonpositive, std::invalid_argument,
270 "Tpetra::"
271 "BlockCrsMatrix constructor: The input blockSize = "
272 << blockSize << " <= 0. The block size must be positive.");
273
274 const size_t blockSizeProd = static_cast<size_t>(blockSize) * static_cast<size_t>(blockSize);
275 TEUCHOS_TEST_FOR_EXCEPTION(
276 blockSizeProd > static_cast<size_t>(std::numeric_limits<LO>::max()),
277 std::overflow_error,
278 "Tpetra::BlockCrsMatrix constructor: blockSize * blockSize overflows LO");
279 offsetPerBlock_ = static_cast<LO>(blockSizeProd);
280
281 domainPointMap_ = BMV::makePointMap(*(graph.getDomainMap()), blockSize);
282 rangePointMap_ = BMV::makePointMap(*(graph.getRangeMap()), blockSize);
283
284 {
285 // These are rcp
286 const auto domainPointMap = getDomainMap();
287 const auto colPointMap = Teuchos::rcp(new typename BMV::map_type(BMV::makePointMap(*graph_.getColMap(), blockSize_)));
288 *pointImporter_ = Teuchos::rcp(new typename crs_graph_type::import_type(domainPointMap, colPointMap));
289 }
290 {
291 auto local_graph_h = graph.getLocalGraphHost();
292 auto ptr_h = local_graph_h.row_map;
293 ptrHost_ = decltype(ptrHost_)(Kokkos::ViewAllocateWithoutInitializing("graph row offset"), ptr_h.extent(0));
294 Kokkos::deep_copy(ptrHost_, ptr_h);
295
296 auto ind_h = local_graph_h.entries;
297 indHost_ = decltype(indHost_)(Kokkos::ViewAllocateWithoutInitializing("graph column indices"), ind_h.extent(0));
298 Kokkos::deep_copy(indHost_, ind_h);
299
300 const auto numValEnt = ind_h.extent(0) * offsetPerBlock();
301 TEUCHOS_ASSERT_EQUALITY(numValEnt, values.size());
302 val_ = decltype(val_)(values);
303 }
304}
305
306template <class Scalar, class LO, class GO, class Node>
307BlockCrsMatrix<Scalar, LO, GO, Node>::
308 BlockCrsMatrix(const crs_graph_type& graph,
309 const map_type& domainPointMap,
310 const map_type& rangePointMap,
311 const LO blockSize)
312 : dist_object_type(graph.getMap())
313 , graph_(graph)
314 , rowMeshMap_(*(graph.getRowMap()))
315 , domainPointMap_(domainPointMap)
316 , rangePointMap_(rangePointMap)
317 , blockSize_(blockSize)
318 , X_colMap_(new Teuchos::RCP<BMV>())
319 , // ptr to a null ptr
320 Y_rowMap_(new Teuchos::RCP<BMV>())
321 , // ptr to a null ptr
322 pointImporter_(new Teuchos::RCP<typename crs_graph_type::import_type>())
323 , offsetPerBlock_(0)
324 , localError_(new bool(false))
325 , errs_(new Teuchos::RCP<std::ostringstream>()) // ptr to a null ptr
326{
327 const size_t blockSizeProd = static_cast<size_t>(blockSize) * static_cast<size_t>(blockSize);
328 TEUCHOS_TEST_FOR_EXCEPTION(
329 blockSizeProd > static_cast<size_t>(std::numeric_limits<LO>::max()),
330 std::overflow_error,
331 "Tpetra::BlockCrsMatrix constructor: blockSize * blockSize overflows LO");
332 offsetPerBlock_ = static_cast<LO>(blockSizeProd);
333
334 TEUCHOS_TEST_FOR_EXCEPTION(
335 !graph_.isSorted(), std::invalid_argument,
336 "Tpetra::"
337 "BlockCrsMatrix constructor: The input CrsGraph does not have sorted "
338 "rows (isSorted() is false). This class assumes sorted rows.");
339
340 graphRCP_ = Teuchos::rcpFromRef(graph_);
341
342 // Trick to test whether LO is nonpositive, without a compiler
343 // warning in case LO is unsigned (which is generally a bad idea
344 // anyway). I don't promise that the trick works, but it
345 // generally does with gcc at least, in my experience.
346 const bool blockSizeIsNonpositive = (blockSize + 1 <= 1);
347 TEUCHOS_TEST_FOR_EXCEPTION(
348 blockSizeIsNonpositive, std::invalid_argument,
349 "Tpetra::"
350 "BlockCrsMatrix constructor: The input blockSize = "
351 << blockSize << " <= 0. The block size must be positive.");
352 {
353 // These are rcp
354 const auto rcpDomainPointMap = getDomainMap();
355 const auto colPointMap = Teuchos::rcp(new typename BMV::map_type(BMV::makePointMap(*graph_.getColMap(), blockSize_)));
356 *pointImporter_ = Teuchos::rcp(new typename crs_graph_type::import_type(rcpDomainPointMap, colPointMap));
357 }
358 {
359 auto local_graph_h = graph.getLocalGraphHost();
360 auto ptr_h = local_graph_h.row_map;
361 ptrHost_ = decltype(ptrHost_)(Kokkos::ViewAllocateWithoutInitializing("graph row offset"), ptr_h.extent(0));
362 // DEEP_COPY REVIEW - HOST-TO-HOST
363 Kokkos::deep_copy(ptrHost_, ptr_h);
364
365 auto ind_h = local_graph_h.entries;
366 indHost_ = decltype(indHost_)(Kokkos::ViewAllocateWithoutInitializing("graph column indices"), ind_h.extent(0));
367 // DEEP_COPY REVIEW - HOST-TO-HOST
368 Kokkos::deep_copy(indHost_, ind_h);
369
370 const auto numValEnt = ind_h.extent(0) * offsetPerBlock();
371 val_ = decltype(val_)(impl_scalar_type_dualview("val", numValEnt));
372 }
373}
374
375template <class Scalar, class LO, class GO, class Node>
376Teuchos::RCP<const typename BlockCrsMatrix<Scalar, LO, GO, Node>::map_type>
377BlockCrsMatrix<Scalar, LO, GO, Node>::
378 getDomainMap() const { // Copy constructor of map_type does a shallow copy.
379 // We're only returning by RCP for backwards compatibility.
380 return Teuchos::rcp(new map_type(domainPointMap_));
381}
382
383template <class Scalar, class LO, class GO, class Node>
384Teuchos::RCP<const typename BlockCrsMatrix<Scalar, LO, GO, Node>::map_type>
385BlockCrsMatrix<Scalar, LO, GO, Node>::
386 getRangeMap() const { // Copy constructor of map_type does a shallow copy.
387 // We're only returning by RCP for backwards compatibility.
388 return Teuchos::rcp(new map_type(rangePointMap_));
389}
390
391template <class Scalar, class LO, class GO, class Node>
392Teuchos::RCP<const typename BlockCrsMatrix<Scalar, LO, GO, Node>::map_type>
393BlockCrsMatrix<Scalar, LO, GO, Node>::
394 getRowMap() const {
395 return graph_.getRowMap();
396}
397
398template <class Scalar, class LO, class GO, class Node>
399Teuchos::RCP<const typename BlockCrsMatrix<Scalar, LO, GO, Node>::map_type>
400BlockCrsMatrix<Scalar, LO, GO, Node>::
401 getColMap() const {
402 return graph_.getColMap();
403}
404
405template <class Scalar, class LO, class GO, class Node>
406global_size_t
407BlockCrsMatrix<Scalar, LO, GO, Node>::
408 getGlobalNumRows() const {
409 return graph_.getGlobalNumRows();
410}
411
412template <class Scalar, class LO, class GO, class Node>
413size_t
414BlockCrsMatrix<Scalar, LO, GO, Node>::
415 getLocalNumRows() const {
416 return graph_.getLocalNumRows();
417}
418
419template <class Scalar, class LO, class GO, class Node>
420size_t
421BlockCrsMatrix<Scalar, LO, GO, Node>::
422 getLocalMaxNumRowEntries() const {
423 return graph_.getLocalMaxNumRowEntries();
424}
425
426template <class Scalar, class LO, class GO, class Node>
427void BlockCrsMatrix<Scalar, LO, GO, Node>::
428 apply(const mv_type& X,
429 mv_type& Y,
430 Teuchos::ETransp mode,
431 Scalar alpha,
432 Scalar beta) const {
433 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
434 TEUCHOS_TEST_FOR_EXCEPTION(
435 mode != Teuchos::NO_TRANS && mode != Teuchos::TRANS && mode != Teuchos::CONJ_TRANS,
436 std::invalid_argument,
437 "Tpetra::BlockCrsMatrix::apply: "
438 "Invalid 'mode' argument. Valid values are Teuchos::NO_TRANS, "
439 "Teuchos::TRANS, and Teuchos::CONJ_TRANS.");
440
441 TEUCHOS_TEST_FOR_EXCEPTION(
442 !X.isConstantStride() || !Y.isConstantStride(),
443 std::invalid_argument,
444 "Tpetra::BlockCrsMatrix::apply: "
445 "X and Y must both be constant stride");
446
447 BMV X_view;
448 BMV Y_view;
449 const LO blockSize = getBlockSize();
450 try {
451 X_view = BMV(X, *(graph_.getDomainMap()), blockSize);
452 Y_view = BMV(Y, *(graph_.getRangeMap()), blockSize);
453 } catch (std::invalid_argument& e) {
454 TEUCHOS_TEST_FOR_EXCEPTION(
455 true, std::invalid_argument,
456 "Tpetra::BlockCrsMatrix::"
457 "apply: Either the input MultiVector X or the output MultiVector Y "
458 "cannot be viewed as a BlockMultiVector, given this BlockCrsMatrix's "
459 "graph. BlockMultiVector's constructor threw the following exception: "
460 << e.what());
461 }
462
463 try {
464 // mfh 16 May 2014: Casting 'this' to nonconst is icky, but it's
465 // either that or mark fields of this class as 'mutable'. The
466 // problem is that applyBlock wants to do lazy initialization of
467 // temporary block multivectors.
468 const_cast<this_BCRS_type&>(*this).applyBlock(X_view, Y_view, mode, alpha, beta);
469 } catch (std::invalid_argument& e) {
470 TEUCHOS_TEST_FOR_EXCEPTION(
471 true, std::invalid_argument,
472 "Tpetra::BlockCrsMatrix::"
473 "apply: The implementation method applyBlock complained about having "
474 "an invalid argument. It reported the following: "
475 << e.what());
476 } catch (std::logic_error& e) {
477 TEUCHOS_TEST_FOR_EXCEPTION(
478 true, std::invalid_argument,
479 "Tpetra::BlockCrsMatrix::"
480 "apply: The implementation method applyBlock complained about a "
481 "possible bug in its implementation. It reported the following: "
482 << e.what());
483 } catch (std::exception& e) {
484 TEUCHOS_TEST_FOR_EXCEPTION(
485 true, std::invalid_argument,
486 "Tpetra::BlockCrsMatrix::"
487 "apply: The implementation method applyBlock threw an exception which "
488 "is neither std::invalid_argument nor std::logic_error, but is a "
489 "subclass of std::exception. It reported the following: "
490 << e.what());
491 } catch (...) {
492 TEUCHOS_TEST_FOR_EXCEPTION(
493 true, std::logic_error,
494 "Tpetra::BlockCrsMatrix::"
495 "apply: The implementation method applyBlock threw an exception which "
496 "is not an instance of a subclass of std::exception. This probably "
497 "indicates a bug. Please report this to the Tpetra developers.");
498 }
499}
500
501template <class Scalar, class LO, class GO, class Node>
502void BlockCrsMatrix<Scalar, LO, GO, Node>::
503 applyBlock(const BlockMultiVector<Scalar, LO, GO, Node>& X,
504 BlockMultiVector<Scalar, LO, GO, Node>& Y,
505 Teuchos::ETransp mode,
506 const Scalar alpha,
507 const Scalar beta) {
508 TEUCHOS_TEST_FOR_EXCEPTION(
509 X.getBlockSize() != Y.getBlockSize(), std::invalid_argument,
510 "Tpetra::BlockCrsMatrix::applyBlock: "
511 "X and Y have different block sizes. X.getBlockSize() = "
512 << X.getBlockSize() << " != Y.getBlockSize() = "
513 << Y.getBlockSize() << ".");
514
515 if (mode == Teuchos::NO_TRANS) {
516 applyBlockNoTrans(X, Y, alpha, beta);
517 } else if (mode == Teuchos::TRANS || mode == Teuchos::CONJ_TRANS) {
518 applyBlockTrans(X, Y, mode, alpha, beta);
519 } else {
520 TEUCHOS_TEST_FOR_EXCEPTION(
521 true, std::invalid_argument,
522 "Tpetra::BlockCrsMatrix::"
523 "applyBlock: Invalid 'mode' argument. Valid values are "
524 "Teuchos::NO_TRANS, Teuchos::TRANS, and Teuchos::CONJ_TRANS.");
525 }
526}
527
528template <class Scalar, class LO, class GO, class Node>
529void BlockCrsMatrix<Scalar, LO, GO, Node>::
530 importAndFillComplete(Teuchos::RCP<BlockCrsMatrix<Scalar, LO, GO, Node> >& destMatrix,
531 const Import<LO, GO, Node>& importer) const {
532 using Teuchos::RCP;
533 using Teuchos::rcp;
534 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
535
536 // Right now, we make many assumptions...
537 TEUCHOS_TEST_FOR_EXCEPTION(!destMatrix.is_null(), std::invalid_argument,
538 "destMatrix is required to be null.");
539
540 // BlockCrsMatrix requires a complete graph at construction.
541 // So first step is to import and fill complete the destGraph.
542 RCP<crs_graph_type> srcGraph = rcp(new crs_graph_type(this->getCrsGraph()));
543 RCP<crs_graph_type> destGraph = importAndFillCompleteCrsGraph<crs_graph_type>(srcGraph, importer,
544 srcGraph->getDomainMap(),
545 srcGraph->getRangeMap());
546
547 // Final step, create and import the destMatrix.
548 destMatrix = rcp(new this_BCRS_type(*destGraph, getBlockSize()));
549 destMatrix->doImport(*this, importer, Tpetra::INSERT);
550}
551
552template <class Scalar, class LO, class GO, class Node>
553void BlockCrsMatrix<Scalar, LO, GO, Node>::
554 exportAndFillComplete(Teuchos::RCP<BlockCrsMatrix<Scalar, LO, GO, Node> >& destMatrix,
555 const Export<LO, GO, Node>& exporter) const {
556 using Teuchos::RCP;
557 using Teuchos::rcp;
558 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
559
560 // Right now, we make many assumptions...
561 TEUCHOS_TEST_FOR_EXCEPTION(!destMatrix.is_null(), std::invalid_argument,
562 "destMatrix is required to be null.");
563
564 // BlockCrsMatrix requires a complete graph at construction.
565 // So first step is to import and fill complete the destGraph.
566 RCP<crs_graph_type> srcGraph = rcp(new crs_graph_type(this->getCrsGraph()));
567 RCP<crs_graph_type> destGraph = exportAndFillCompleteCrsGraph<crs_graph_type>(srcGraph, exporter,
568 srcGraph->getDomainMap(),
569 srcGraph->getRangeMap());
570
571 // Final step, create and import the destMatrix.
572 destMatrix = rcp(new this_BCRS_type(*destGraph, getBlockSize()));
573 destMatrix->doExport(*this, exporter, Tpetra::INSERT);
574}
575
576template <class Scalar, class LO, class GO, class Node>
577void BlockCrsMatrix<Scalar, LO, GO, Node>::
578 setAllToScalar(const Scalar& alpha) {
579 auto val_d = val_.getDeviceView(Access::OverwriteAll);
580 Kokkos::deep_copy(val_d, alpha);
581}
582
583template <class Scalar, class LO, class GO, class Node>
584LO BlockCrsMatrix<Scalar, LO, GO, Node>::
585 replaceLocalValues(const LO localRowInd,
586 const LO colInds[],
587 const Scalar vals[],
588 const LO numColInds) const {
589 std::vector<ptrdiff_t> offsets(numColInds);
590 const LO numOffsets = this->getLocalRowOffsets(localRowInd, offsets.data(), colInds, numColInds);
591 const LO validCount = this->replaceLocalValuesByOffsets(localRowInd, offsets.data(), vals, numOffsets);
592 return validCount;
593}
594
595template <class Scalar, class LO, class GO, class Node>
596void BlockCrsMatrix<Scalar, LO, GO, Node>::
597 getLocalDiagOffsets(const Kokkos::View<size_t*, device_type,
598 Kokkos::MemoryUnmanaged>& offsets) const {
599 graph_.getLocalDiagOffsets(offsets);
600}
601
602template <class Scalar, class LO, class GO, class Node>
603void BlockCrsMatrix<Scalar, LO, GO, Node>::
604 getLocalDiagCopy(const Kokkos::View<impl_scalar_type***, device_type,
605 Kokkos::MemoryUnmanaged>& diag,
606 const Kokkos::View<const size_t*, device_type,
607 Kokkos::MemoryUnmanaged>& offsets) const {
608 using Kokkos::parallel_for;
609 const char prefix[] = "Tpetra::BlockCrsMatrix::getLocalDiagCopy (2-arg): ";
610
611 const LO lclNumMeshRows = static_cast<LO>(rowMeshMap_.getLocalNumElements());
612 const LO blockSize = this->getBlockSize();
613 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<LO>(diag.extent(0)) < lclNumMeshRows ||
614 static_cast<LO>(diag.extent(1)) < blockSize ||
615 static_cast<LO>(diag.extent(2)) < blockSize,
616 std::invalid_argument, prefix << "The input Kokkos::View is not big enough to hold all the data.");
617 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<LO>(offsets.size()) < lclNumMeshRows, std::invalid_argument,
618 prefix << "offsets.size() = " << offsets.size() << " < local number of "
619 "diagonal blocks "
620 << lclNumMeshRows << ".");
621
622 typedef Kokkos::RangePolicy<execution_space, LO> policy_type;
623 typedef GetLocalDiagCopy<Scalar, LO, GO, Node> functor_type;
624
625 // FIXME (mfh 26 May 2016) Not really OK to const_cast here, since
626 // we reserve the right to do lazy allocation of device data. (We
627 // don't plan to do lazy allocation for host data; the host
628 // version of the data always exists.)
629 auto val_d = val_.getDeviceView(Access::ReadOnly);
630 parallel_for(policy_type(0, lclNumMeshRows),
631 functor_type(diag, val_d, offsets,
632 graph_.getLocalGraphDevice().row_map, blockSize_));
633}
634
635template <class Scalar, class LO, class GO, class Node>
636LO BlockCrsMatrix<Scalar, LO, GO, Node>::
637 absMaxLocalValues(const LO localRowInd,
638 const LO colInds[],
639 const Scalar vals[],
640 const LO numColInds) const {
641 std::vector<ptrdiff_t> offsets(numColInds);
642 const LO numOffsets = this->getLocalRowOffsets(localRowInd, offsets.data(), colInds, numColInds);
643 const LO validCount = this->absMaxLocalValuesByOffsets(localRowInd, offsets.data(), vals, numOffsets);
644 return validCount;
645}
646
647template <class Scalar, class LO, class GO, class Node>
648LO BlockCrsMatrix<Scalar, LO, GO, Node>::
649 sumIntoLocalValues(const LO localRowInd,
650 const LO colInds[],
651 const Scalar vals[],
652 const LO numColInds) const {
653 std::vector<ptrdiff_t> offsets(numColInds);
654 const LO numOffsets = this->getLocalRowOffsets(localRowInd, offsets.data(), colInds, numColInds);
655 const LO validCount = this->sumIntoLocalValuesByOffsets(localRowInd, offsets.data(), vals, numOffsets);
656 return validCount;
657}
658template <class Scalar, class LO, class GO, class Node>
659void BlockCrsMatrix<Scalar, LO, GO, Node>::
660 getLocalRowCopy(LO LocalRow,
661 nonconst_local_inds_host_view_type& Indices,
662 nonconst_values_host_view_type& Values,
663 size_t& NumEntries) const {
664 auto vals = getValuesHost(LocalRow);
665 const LO numInds = ptrHost_(LocalRow + 1) - ptrHost_(LocalRow);
666 if (numInds > (LO)Indices.extent(0) || numInds * blockSize_ * blockSize_ > (LO)Values.extent(0)) {
667 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
668 "Tpetra::BlockCrsMatrix::getLocalRowCopy : Column and/or values array is not large enough to hold "
669 << numInds << " row entries");
670 }
671 const LO* colInds = indHost_.data() + ptrHost_(LocalRow);
672 for (LO i = 0; i < numInds; ++i) {
673 Indices[i] = colInds[i];
674 }
675 for (LO i = 0; i < numInds * blockSize_ * blockSize_; ++i) {
676 Values[i] = vals[i];
677 }
678 NumEntries = numInds;
679}
680
681template <class Scalar, class LO, class GO, class Node>
682LO BlockCrsMatrix<Scalar, LO, GO, Node>::
683 getLocalRowOffsets(const LO localRowInd,
684 ptrdiff_t offsets[],
685 const LO colInds[],
686 const LO numColInds) const {
687 if (!rowMeshMap_.isNodeLocalElement(localRowInd)) {
688 // We got no offsets, because the input local row index is
689 // invalid on the calling process. That may not be an error, if
690 // numColInds is zero anyway; it doesn't matter. This is the
691 // advantage of returning the number of valid indices.
692 return static_cast<LO>(0);
693 }
694
695 const LO LINV = Teuchos::OrdinalTraits<LO>::invalid();
696 LO hint = 0; // Guess for the relative offset into the current row
697 LO validCount = 0; // number of valid column indices in colInds
698
699 for (LO k = 0; k < numColInds; ++k) {
700 const LO relBlockOffset =
701 this->findRelOffsetOfColumnIndex(localRowInd, colInds[k], hint);
702 if (relBlockOffset != LINV) {
703 offsets[k] = static_cast<ptrdiff_t>(relBlockOffset);
704 hint = relBlockOffset + 1;
705 ++validCount;
706 }
707 }
708 return validCount;
709}
710
711template <class Scalar, class LO, class GO, class Node>
712LO BlockCrsMatrix<Scalar, LO, GO, Node>::
713 replaceLocalValuesByOffsets(const LO localRowInd,
714 const ptrdiff_t offsets[],
715 const Scalar vals[],
716 const LO numOffsets) const {
717 if (!rowMeshMap_.isNodeLocalElement(localRowInd)) {
718 // We modified no values, because the input local row index is
719 // invalid on the calling process. That may not be an error, if
720 // numColInds is zero anyway; it doesn't matter. This is the
721 // advantage of returning the number of valid indices.
722 return static_cast<LO>(0);
723 }
724 const impl_scalar_type* const vIn = reinterpret_cast<const impl_scalar_type*>(vals);
725 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
726 auto val_out = const_cast<this_BCRS_type&>(*this).getValuesHostNonConst(localRowInd);
727 impl_scalar_type* vOut = val_out.data();
728
729 const size_t perBlockSize = offsetPerBlock();
730 const ptrdiff_t STINV = Teuchos::OrdinalTraits<ptrdiff_t>::invalid();
731 size_t pointOffset = 0; // Current offset into input values
732 LO validCount = 0; // number of valid offsets
733
734 for (LO k = 0; k < numOffsets; ++k, pointOffset += perBlockSize) {
735 const size_t blockOffset = offsets[k] * perBlockSize;
736 if (offsets[k] != STINV) {
737 little_block_type A_old = getNonConstLocalBlockFromInput(vOut, blockOffset);
738 const_little_block_type A_new = getConstLocalBlockFromInput(vIn, pointOffset);
739 COPY(A_new, A_old);
740 ++validCount;
741 }
742 }
743 return validCount;
744}
745
746template <class Scalar, class LO, class GO, class Node>
747LO BlockCrsMatrix<Scalar, LO, GO, Node>::
748 absMaxLocalValuesByOffsets(const LO localRowInd,
749 const ptrdiff_t offsets[],
750 const Scalar vals[],
751 const LO numOffsets) const {
752 if (!rowMeshMap_.isNodeLocalElement(localRowInd)) {
753 // We modified no values, because the input local row index is
754 // invalid on the calling process. That may not be an error, if
755 // numColInds is zero anyway; it doesn't matter. This is the
756 // advantage of returning the number of valid indices.
757 return static_cast<LO>(0);
758 }
759 const impl_scalar_type* const vIn = reinterpret_cast<const impl_scalar_type*>(vals);
760 auto val_out = getValuesHost(localRowInd);
761 impl_scalar_type* vOut = const_cast<impl_scalar_type*>(val_out.data());
762
763 const size_t perBlockSize = offsetPerBlock();
764 const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
765 size_t pointOffset = 0; // Current offset into input values
766 LO validCount = 0; // number of valid offsets
767
768 for (LO k = 0; k < numOffsets; ++k, pointOffset += perBlockSize) {
769 const size_t blockOffset = offsets[k] * perBlockSize;
770 if (blockOffset != STINV) {
771 little_block_type A_old = getNonConstLocalBlockFromInput(vOut, blockOffset);
772 const_little_block_type A_new = getConstLocalBlockFromInput(vIn, pointOffset);
773 ::Tpetra::Impl::absMax(A_old, A_new);
774 ++validCount;
775 }
776 }
777 return validCount;
778}
779
780template <class Scalar, class LO, class GO, class Node>
781LO BlockCrsMatrix<Scalar, LO, GO, Node>::
782 sumIntoLocalValuesByOffsets(const LO localRowInd,
783 const ptrdiff_t offsets[],
784 const Scalar vals[],
785 const LO numOffsets) const {
786 if (!rowMeshMap_.isNodeLocalElement(localRowInd)) {
787 // We modified no values, because the input local row index is
788 // invalid on the calling process. That may not be an error, if
789 // numColInds is zero anyway; it doesn't matter. This is the
790 // advantage of returning the number of valid indices.
791 return static_cast<LO>(0);
792 }
793 const impl_scalar_type ONE = static_cast<impl_scalar_type>(1.0);
794 const impl_scalar_type* const vIn = reinterpret_cast<const impl_scalar_type*>(vals);
795 typedef BlockCrsMatrix<Scalar, LO, GO, Node> this_BCRS_type;
796 auto val_out = const_cast<this_BCRS_type&>(*this).getValuesHostNonConst(localRowInd);
797 impl_scalar_type* vOut = val_out.data();
798
799 const size_t perBlockSize = offsetPerBlock();
800 const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
801 size_t pointOffset = 0; // Current offset into input values
802 LO validCount = 0; // number of valid offsets
803
804 for (LO k = 0; k < numOffsets; ++k, pointOffset += perBlockSize) {
805 const size_t blockOffset = offsets[k] * perBlockSize;
806 if (blockOffset != STINV) {
807 little_block_type A_old = getNonConstLocalBlockFromInput(vOut, blockOffset);
808 const_little_block_type A_new = getConstLocalBlockFromInput(vIn, pointOffset);
809 AXPY(ONE, A_new, A_old);
810 ++validCount;
811 }
812 }
813 return validCount;
814}
815
816template <class Scalar, class LO, class GO, class Node>
817typename BlockCrsMatrix<Scalar, LO, GO, Node>::
818 impl_scalar_type_dualview::t_host::const_type
819 BlockCrsMatrix<Scalar, LO, GO, Node>::
820 getValuesHost() const {
821 return val_.getHostView(Access::ReadOnly);
822}
823
824template <class Scalar, class LO, class GO, class Node>
825typename BlockCrsMatrix<Scalar, LO, GO, Node>::
826 impl_scalar_type_dualview::t_dev::const_type
827 BlockCrsMatrix<Scalar, LO, GO, Node>::
828 getValuesDevice() const {
829 return val_.getDeviceView(Access::ReadOnly);
830}
831
832template <class Scalar, class LO, class GO, class Node>
833typename BlockCrsMatrix<Scalar, LO, GO, Node>::
834 impl_scalar_type_dualview::t_host
835 BlockCrsMatrix<Scalar, LO, GO, Node>::
836 getValuesHostNonConst() const {
837 return val_.getHostView(Access::ReadWrite);
838}
839
840template <class Scalar, class LO, class GO, class Node>
841typename BlockCrsMatrix<Scalar, LO, GO, Node>::
842 impl_scalar_type_dualview::t_dev
843 BlockCrsMatrix<Scalar, LO, GO, Node>::
844 getValuesDeviceNonConst() const {
845 return val_.getDeviceView(Access::ReadWrite);
846}
847
848template <class Scalar, class LO, class GO, class Node>
849typename BlockCrsMatrix<Scalar, LO, GO, Node>::
850 impl_scalar_type_dualview::t_host::const_type
851 BlockCrsMatrix<Scalar, LO, GO, Node>::
852 getValuesHost(const LO& lclRow) const {
853 const size_t perBlockSize = offsetPerBlock();
854 auto val = val_.getHostView(Access::ReadOnly);
855 auto r_val = Kokkos::subview(val, Kokkos::pair<size_t, size_t>(ptrHost_(lclRow) * perBlockSize, ptrHost_(lclRow + 1) * perBlockSize));
856 return r_val;
857}
858
859template <class Scalar, class LO, class GO, class Node>
860typename BlockCrsMatrix<Scalar, LO, GO, Node>::
861 impl_scalar_type_dualview::t_dev::const_type
862 BlockCrsMatrix<Scalar, LO, GO, Node>::
863 getValuesDevice(const LO& lclRow) const {
864 const size_t perBlockSize = offsetPerBlock();
865 auto val = val_.getDeviceView(Access::ReadOnly);
866 auto r_val = Kokkos::subview(val, Kokkos::pair<size_t, size_t>(ptrHost_(lclRow) * perBlockSize, ptrHost_(lclRow + 1) * perBlockSize));
867 return r_val;
868}
869
870template <class Scalar, class LO, class GO, class Node>
871typename BlockCrsMatrix<Scalar, LO, GO, Node>::impl_scalar_type_dualview::t_host
872BlockCrsMatrix<Scalar, LO, GO, Node>::
873 getValuesHostNonConst(const LO& lclRow) {
874 const size_t perBlockSize = offsetPerBlock();
875 auto val = val_.getHostView(Access::ReadWrite);
876 auto r_val = Kokkos::subview(val, Kokkos::pair<size_t, size_t>(ptrHost_(lclRow) * perBlockSize, ptrHost_(lclRow + 1) * perBlockSize));
877 return r_val;
878}
879
880template <class Scalar, class LO, class GO, class Node>
881typename BlockCrsMatrix<Scalar, LO, GO, Node>::impl_scalar_type_dualview::t_dev
882BlockCrsMatrix<Scalar, LO, GO, Node>::
883 getValuesDeviceNonConst(const LO& lclRow) {
884 const size_t perBlockSize = offsetPerBlock();
885 auto val = val_.getDeviceView(Access::ReadWrite);
886 auto r_val = Kokkos::subview(val, Kokkos::pair<size_t, size_t>(ptrHost_(lclRow) * perBlockSize, ptrHost_(lclRow + 1) * perBlockSize));
887 return r_val;
888}
889
890template <class Scalar, class LO, class GO, class Node>
891size_t
892BlockCrsMatrix<Scalar, LO, GO, Node>::
893 getNumEntriesInLocalRow(const LO localRowInd) const {
894 const size_t numEntInGraph = graph_.getNumEntriesInLocalRow(localRowInd);
895 if (numEntInGraph == Teuchos::OrdinalTraits<size_t>::invalid()) {
896 return 0; // the calling process doesn't have that row
897 } else {
898 return numEntInGraph;
899 }
900}
901
902template <class Scalar, class LO, class GO, class Node>
903typename BlockCrsMatrix<Scalar, LO, GO, Node>::local_matrix_device_type
904BlockCrsMatrix<Scalar, LO, GO, Node>::
905 getLocalMatrixDevice() const {
906 auto numCols = this->graph_.getColMap()->getLocalNumElements();
907 auto val = val_.getDeviceView(Access::ReadWrite);
908 const LO blockSize = this->getBlockSize();
909 const auto graph = this->graph_.getLocalGraphDevice();
910
911 return local_matrix_device_type("Tpetra::BlockCrsMatrix::lclMatrixDevice",
912 numCols,
913 val,
914 graph,
915 blockSize);
916}
917
918template <class Scalar, class LO, class GO, class Node>
919void BlockCrsMatrix<Scalar, LO, GO, Node>::
920 applyBlockTrans(const BlockMultiVector<Scalar, LO, GO, Node>& X,
921 BlockMultiVector<Scalar, LO, GO, Node>& Y,
922 const Teuchos::ETransp mode,
923 const Scalar alpha,
924 const Scalar beta) {
925 (void)X;
926 (void)Y;
927 (void)mode;
928 (void)alpha;
929 (void)beta;
930
931 TEUCHOS_TEST_FOR_EXCEPTION(
932 true, std::logic_error,
933 "Tpetra::BlockCrsMatrix::apply: "
934 "transpose and conjugate transpose modes are not yet implemented.");
935}
936
937template <class Scalar, class LO, class GO, class Node>
938void BlockCrsMatrix<Scalar, LO, GO, Node>::
939 applyBlockNoTrans(const BlockMultiVector<Scalar, LO, GO, Node>& X,
940 BlockMultiVector<Scalar, LO, GO, Node>& Y,
941 const Scalar alpha,
942 const Scalar beta) {
943 using Teuchos::RCP;
944 using Teuchos::rcp;
945 typedef ::Tpetra::Import<LO, GO, Node> import_type;
946 typedef ::Tpetra::Export<LO, GO, Node> export_type;
947 const Scalar zero = STS::zero();
948 const Scalar one = STS::one();
949 RCP<const import_type> import = graph_.getImporter();
950 // "export" is a reserved C++ keyword, so we can't use it.
951 RCP<const export_type> theExport = graph_.getExporter();
952 const char prefix[] = "Tpetra::BlockCrsMatrix::applyBlockNoTrans: ";
953
954 if (alpha == zero) {
955 if (beta == zero) {
956 Y.putScalar(zero); // replace Inf or NaN (BLAS rules)
957 } else if (beta != one) {
958 Y.scale(beta);
959 }
960 } else { // alpha != 0
961 const BMV* X_colMap = NULL;
962 if (import.is_null()) {
963 try {
964 X_colMap = &X;
965 } catch (std::exception& e) {
966 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, prefix << "Tpetra::MultiVector::"
967 "operator= threw an exception: "
968 << std::endl
969 << e.what());
970 }
971 } else {
972 // X_colMap_ is a pointer to a pointer to BMV. Ditto for
973 // Y_rowMap_ below. This lets us do lazy initialization
974 // correctly with view semantics of BlockCrsMatrix. All views
975 // of this BlockCrsMatrix have the same outer pointer. That
976 // way, we can set the inner pointer in one view, and all
977 // other views will see it.
978 if ((*X_colMap_).is_null() ||
979 (**X_colMap_).getNumVectors() != X.getNumVectors() ||
980 (**X_colMap_).getBlockSize() != X.getBlockSize()) {
981 *X_colMap_ = rcp(new BMV(*(graph_.getColMap()), getBlockSize(),
982 static_cast<LO>(X.getNumVectors())));
983 }
984 (*X_colMap_)->getMultiVectorView().doImport(X.getMultiVectorView(), **pointImporter_, ::Tpetra::REPLACE);
985 try {
986 X_colMap = &(**X_colMap_);
987 } catch (std::exception& e) {
988 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, prefix << "Tpetra::MultiVector::"
989 "operator= threw an exception: "
990 << std::endl
991 << e.what());
992 }
993 }
994
995 BMV* Y_rowMap = NULL;
996 if (theExport.is_null()) {
997 Y_rowMap = &Y;
998 } else if ((*Y_rowMap_).is_null() ||
999 (**Y_rowMap_).getNumVectors() != Y.getNumVectors() ||
1000 (**Y_rowMap_).getBlockSize() != Y.getBlockSize()) {
1001 *Y_rowMap_ = rcp(new BMV(*(graph_.getRowMap()), getBlockSize(),
1002 static_cast<LO>(X.getNumVectors())));
1003 try {
1004 Y_rowMap = &(**Y_rowMap_);
1005 } catch (std::exception& e) {
1006 TEUCHOS_TEST_FOR_EXCEPTION(
1007 true, std::logic_error, prefix << "Tpetra::MultiVector::"
1008 "operator= threw an exception: "
1009 << std::endl
1010 << e.what());
1011 }
1012 }
1013
1014 try {
1015 localApplyBlockNoTrans(*X_colMap, *Y_rowMap, alpha, beta);
1016 } catch (std::exception& e) {
1017 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, prefix << "localApplyBlockNoTrans threw "
1018 "an exception: "
1019 << e.what());
1020 } catch (...) {
1021 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, prefix << "localApplyBlockNoTrans threw "
1022 "an exception not a subclass of std::exception.");
1023 }
1024
1025 if (!theExport.is_null()) {
1026 Y.doExport(*Y_rowMap, *theExport, ::Tpetra::REPLACE);
1027 }
1028 }
1029}
1030
1031template <class Scalar, class LO, class GO, class Node>
1032void BlockCrsMatrix<Scalar, LO, GO, Node>::localApplyBlockNoTrans(
1033 const BlockMultiVector<Scalar, LO, GO, Node>& X,
1034 BlockMultiVector<Scalar, LO, GO, Node>& Y, const Scalar alpha,
1035 const Scalar beta) {
1036 ::Tpetra::Details::ProfilingRegion profile_region(
1037 "Tpetra::BlockCrsMatrix::localApplyBlockNoTrans");
1038 const impl_scalar_type alpha_impl = alpha;
1039 const auto graph = this->graph_.getLocalGraphDevice();
1040
1041 mv_type X_mv = X.getMultiVectorView();
1042 mv_type Y_mv = Y.getMultiVectorView();
1043 auto X_lcl = X_mv.getLocalViewDevice(Access::ReadOnly);
1044 auto Y_lcl = Y_mv.getLocalViewDevice(Access::ReadWrite);
1045
1046 auto A_lcl = getLocalMatrixDevice();
1047 if (!applyHelper.get()) {
1048 // The apply helper does not exist, so create it
1049 applyHelper = std::make_shared<ApplyHelper>(A_lcl.nnz(), A_lcl.graph.row_map);
1050 }
1051 if (applyHelper->shouldUseIntRowptrs()) {
1052 auto A_lcl_int_rowptrs = applyHelper->getIntRowptrMatrix(A_lcl);
1053 KokkosSparse::spmv(
1054 &applyHelper->handle_int, KokkosSparse::NoTranspose,
1055 alpha_impl, A_lcl_int_rowptrs, X_lcl, beta, Y_lcl);
1056 } else {
1057 KokkosSparse::spmv(
1058 &applyHelper->handle, KokkosSparse::NoTranspose,
1059 alpha_impl, A_lcl, X_lcl, beta, Y_lcl);
1060 }
1061}
1062
1063template <class Scalar, class LO, class GO, class Node>
1064LO BlockCrsMatrix<Scalar, LO, GO, Node>::
1065 findRelOffsetOfColumnIndex(const LO localRowIndex,
1066 const LO colIndexToFind,
1067 const LO hint) const {
1068 const size_t absStartOffset = ptrHost_[localRowIndex];
1069 const size_t absEndOffset = ptrHost_[localRowIndex + 1];
1070 const LO numEntriesInRow = static_cast<LO>(absEndOffset - absStartOffset);
1071 // Amortize pointer arithmetic over the search loop.
1072 const LO* const curInd = indHost_.data() + absStartOffset;
1073
1074 // If the hint was correct, then the hint is the offset to return.
1075 if (hint < numEntriesInRow && curInd[hint] == colIndexToFind) {
1076 // Always return the offset relative to the current row.
1077 return hint;
1078 }
1079
1080 // The hint was wrong, so we must search for the given column
1081 // index in the column indices for the given row.
1082 LO relOffset = Teuchos::OrdinalTraits<LO>::invalid();
1083
1084 // We require that the graph have sorted rows. However, binary
1085 // search only pays if the current row is longer than a certain
1086 // amount. We set this to 32, but you might want to tune this.
1087 const LO maxNumEntriesForLinearSearch = 32;
1088 if (numEntriesInRow > maxNumEntriesForLinearSearch) {
1089 // Use binary search. It would probably be better for us to
1090 // roll this loop by hand. If we wrote it right, a smart
1091 // compiler could perhaps use conditional loads and avoid
1092 // branches (according to Jed Brown on May 2014).
1093 const LO* beg = curInd;
1094 const LO* end = curInd + numEntriesInRow;
1095 std::pair<const LO*, const LO*> p =
1096 std::equal_range(beg, end, colIndexToFind);
1097 if (p.first != p.second) {
1098 // offset is relative to the current row
1099 relOffset = static_cast<LO>(p.first - beg);
1100 }
1101 } else { // use linear search
1102 for (LO k = 0; k < numEntriesInRow; ++k) {
1103 if (colIndexToFind == curInd[k]) {
1104 relOffset = k; // offset is relative to the current row
1105 break;
1106 }
1107 }
1108 }
1109
1110 return relOffset;
1111}
1112
1113template <class Scalar, class LO, class GO, class Node>
1114LO BlockCrsMatrix<Scalar, LO, GO, Node>::
1115 offsetPerBlock() const {
1116 return offsetPerBlock_;
1117}
1118
1119template <class Scalar, class LO, class GO, class Node>
1120typename BlockCrsMatrix<Scalar, LO, GO, Node>::const_little_block_type
1121BlockCrsMatrix<Scalar, LO, GO, Node>::
1122 getConstLocalBlockFromInput(const impl_scalar_type* val,
1123 const size_t pointOffset) const {
1124 // Row major blocks
1125 const LO rowStride = blockSize_;
1126 return const_little_block_type(val + pointOffset, blockSize_, rowStride);
1127}
1128
1129template <class Scalar, class LO, class GO, class Node>
1130typename BlockCrsMatrix<Scalar, LO, GO, Node>::little_block_type
1131BlockCrsMatrix<Scalar, LO, GO, Node>::
1132 getNonConstLocalBlockFromInput(impl_scalar_type* val,
1133 const size_t pointOffset) const {
1134 // Row major blocks
1135 const LO rowStride = blockSize_;
1136 return little_block_type(val + pointOffset, blockSize_, rowStride);
1137}
1138
1139template <class Scalar, class LO, class GO, class Node>
1140typename BlockCrsMatrix<Scalar, LO, GO, Node>::little_block_host_type
1141BlockCrsMatrix<Scalar, LO, GO, Node>::
1142 getNonConstLocalBlockFromInputHost(impl_scalar_type* val,
1143 const size_t pointOffset) const {
1144 // Row major blocks
1145 const LO rowStride = blockSize_;
1146 const size_t bs2 = blockSize_ * blockSize_;
1147 return little_block_host_type(val + bs2 * pointOffset, blockSize_, rowStride);
1148}
1149
1150template <class Scalar, class LO, class GO, class Node>
1151typename BlockCrsMatrix<Scalar, LO, GO, Node>::little_block_type
1152BlockCrsMatrix<Scalar, LO, GO, Node>::
1153 getLocalBlockDeviceNonConst(const LO localRowInd, const LO localColInd) const {
1154 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
1155
1156 const size_t absRowBlockOffset = ptrHost_[localRowInd];
1157 const LO relBlockOffset = this->findRelOffsetOfColumnIndex(localRowInd, localColInd);
1158 if (relBlockOffset != Teuchos::OrdinalTraits<LO>::invalid()) {
1159 const size_t absBlockOffset = absRowBlockOffset + relBlockOffset;
1160 auto vals = const_cast<this_BCRS_type&>(*this).getValuesDeviceNonConst();
1161 auto r_val = getNonConstLocalBlockFromInput(vals.data(), absBlockOffset);
1162 return r_val;
1163 } else {
1164 return little_block_type();
1165 }
1166}
1167
1168template <class Scalar, class LO, class GO, class Node>
1169typename BlockCrsMatrix<Scalar, LO, GO, Node>::little_block_host_type
1170BlockCrsMatrix<Scalar, LO, GO, Node>::
1171 getLocalBlockHostNonConst(const LO localRowInd, const LO localColInd) const {
1172 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
1173
1174 const size_t absRowBlockOffset = ptrHost_[localRowInd];
1175 const LO relBlockOffset = this->findRelOffsetOfColumnIndex(localRowInd, localColInd);
1176 if (relBlockOffset != Teuchos::OrdinalTraits<LO>::invalid()) {
1177 const size_t absBlockOffset = absRowBlockOffset + relBlockOffset;
1178 auto vals = const_cast<this_BCRS_type&>(*this).getValuesHostNonConst();
1179 auto r_val = getNonConstLocalBlockFromInputHost(vals.data(), absBlockOffset);
1180 return r_val;
1181 } else {
1182 return little_block_host_type();
1183 }
1184}
1185
1186template <class Scalar, class LO, class GO, class Node>
1187bool BlockCrsMatrix<Scalar, LO, GO, Node>::
1188 checkSizes(const ::Tpetra::SrcDistObject& source) {
1189 using std::endl;
1190 typedef BlockCrsMatrix<Scalar, LO, GO, Node> this_BCRS_type;
1191 const this_BCRS_type* src = dynamic_cast<const this_BCRS_type*>(&source);
1192
1193 if (src == NULL) {
1194 std::ostream& err = markLocalErrorAndGetStream();
1195 err << "checkSizes: The source object of the Import or Export "
1196 "must be a BlockCrsMatrix with the same template parameters as the "
1197 "target object."
1198 << endl;
1199 } else {
1200 // Use a string of ifs, not if-elseifs, because we want to know
1201 // all the errors.
1202 if (src->getBlockSize() != this->getBlockSize()) {
1203 std::ostream& err = markLocalErrorAndGetStream();
1204 err << "checkSizes: The source and target objects of the Import or "
1205 << "Export must have the same block sizes. The source's block "
1206 << "size = " << src->getBlockSize() << " != the target's block "
1207 << "size = " << this->getBlockSize() << "." << endl;
1208 }
1209 if (!src->graph_.isFillComplete()) {
1210 std::ostream& err = markLocalErrorAndGetStream();
1211 err << "checkSizes: The source object of the Import or Export is "
1212 "not fill complete. Both source and target objects must be fill "
1213 "complete."
1214 << endl;
1215 }
1216 if (!this->graph_.isFillComplete()) {
1217 std::ostream& err = markLocalErrorAndGetStream();
1218 err << "checkSizes: The target object of the Import or Export is "
1219 "not fill complete. Both source and target objects must be fill "
1220 "complete."
1221 << endl;
1222 }
1223 if (src->graph_.getColMap().is_null()) {
1224 std::ostream& err = markLocalErrorAndGetStream();
1225 err << "checkSizes: The source object of the Import or Export does "
1226 "not have a column Map. Both source and target objects must have "
1227 "column Maps."
1228 << endl;
1229 }
1230 if (this->graph_.getColMap().is_null()) {
1231 std::ostream& err = markLocalErrorAndGetStream();
1232 err << "checkSizes: The target object of the Import or Export does "
1233 "not have a column Map. Both source and target objects must have "
1234 "column Maps."
1235 << endl;
1236 }
1237 }
1238
1239 return !(*(this->localError_));
1240}
1241
1242template <class Scalar, class LO, class GO, class Node>
1243void BlockCrsMatrix<Scalar, LO, GO, Node>::
1244 copyAndPermute(const ::Tpetra::SrcDistObject& source,
1245 const size_t numSameIDs,
1246 const Kokkos::DualView<const local_ordinal_type*,
1247 buffer_device_type>& permuteToLIDs,
1248 const Kokkos::DualView<const local_ordinal_type*,
1249 buffer_device_type>& permuteFromLIDs,
1250 const CombineMode /*CM*/) {
1251 using std::endl;
1252 using ::Tpetra::Details::Behavior;
1253 using ::Tpetra::Details::dualViewStatusToString;
1254 using ::Tpetra::Details::ProfilingRegion;
1255 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
1256
1257 ProfilingRegion profile_region("Tpetra::BlockCrsMatrix::copyAndPermute");
1258 const bool debug = Behavior::debug();
1259 const bool verbose = Behavior::verbose();
1260
1261 // Define this function prefix
1262 std::string prefix;
1263 {
1264 std::ostringstream os;
1265 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1266 os << "Proc " << myRank << ": BlockCrsMatrix::copyAndPermute : " << endl;
1267 prefix = os.str();
1268 }
1269
1270 // check if this already includes a local error
1271 if (*(this->localError_)) {
1272 std::ostream& err = this->markLocalErrorAndGetStream();
1273 err << prefix
1274 << "The target object of the Import or Export is already in an error state."
1275 << endl;
1276 return;
1277 }
1278
1279 //
1280 // Verbose input dual view status
1281 //
1282 if (verbose) {
1283 std::ostringstream os;
1284 os << prefix << endl
1285 << prefix << " " << dualViewStatusToString(permuteToLIDs, "permuteToLIDs") << endl
1286 << prefix << " " << dualViewStatusToString(permuteFromLIDs, "permuteFromLIDs") << endl;
1287 std::cerr << os.str();
1288 }
1289
1293 if (permuteToLIDs.extent(0) != permuteFromLIDs.extent(0)) {
1294 std::ostream& err = this->markLocalErrorAndGetStream();
1295 err << prefix
1296 << "permuteToLIDs.extent(0) = " << permuteToLIDs.extent(0)
1297 << " != permuteFromLIDs.extent(0) = " << permuteFromLIDs.extent(0)
1298 << "." << endl;
1299 return;
1300 }
1301 if (permuteToLIDs.need_sync_host() || permuteFromLIDs.need_sync_host()) {
1302 std::ostream& err = this->markLocalErrorAndGetStream();
1303 err << prefix
1304 << "Both permuteToLIDs and permuteFromLIDs must be sync'd to host."
1305 << endl;
1306 return;
1307 }
1308
1309 const this_BCRS_type* src = dynamic_cast<const this_BCRS_type*>(&source);
1310 if (src == NULL) {
1311 std::ostream& err = this->markLocalErrorAndGetStream();
1312 err << prefix
1313 << "The source (input) object of the Import or "
1314 "Export is either not a BlockCrsMatrix, or does not have the right "
1315 "template parameters. checkSizes() should have caught this. "
1316 "Please report this bug to the Tpetra developers."
1317 << endl;
1318 return;
1319 }
1320
1321 bool lclErr = false;
1322#ifdef HAVE_TPETRA_DEBUG
1323 std::set<LO> invalidSrcCopyRows;
1324 std::set<LO> invalidDstCopyRows;
1325 std::set<LO> invalidDstCopyCols;
1326 std::set<LO> invalidDstPermuteCols;
1327 std::set<LO> invalidPermuteFromRows;
1328#endif // HAVE_TPETRA_DEBUG
1329
1330 // Copy the initial sequence of rows that are the same.
1331 //
1332 // The two graphs might have different column Maps, so we need to
1333 // do this using global column indices. This is purely local, so
1334 // we only need to check for local sameness of the two column
1335 // Maps.
1336
1337#ifdef HAVE_TPETRA_DEBUG
1338 const map_type& srcRowMap = *(src->graph_.getRowMap());
1339#endif // HAVE_TPETRA_DEBUG
1340 const map_type& dstRowMap = *(this->graph_.getRowMap());
1341 const map_type& srcColMap = *(src->graph_.getColMap());
1342 const map_type& dstColMap = *(this->graph_.getColMap());
1343 const bool canUseLocalColumnIndices = srcColMap.locallySameAs(dstColMap);
1344
1345 const size_t numPermute = static_cast<size_t>(permuteFromLIDs.extent(0));
1346 if (verbose) {
1347 std::ostringstream os;
1348 os << prefix
1349 << "canUseLocalColumnIndices: "
1350 << (canUseLocalColumnIndices ? "true" : "false")
1351 << ", numPermute: " << numPermute
1352 << endl;
1353 std::cerr << os.str();
1354 }
1355
1356 const auto permuteToLIDsHost = permuteToLIDs.view_host();
1357 const auto permuteFromLIDsHost = permuteFromLIDs.view_host();
1358
1359 if (canUseLocalColumnIndices) {
1360 // Copy local rows that are the "same" in both source and target.
1361 for (LO localRow = 0; localRow < static_cast<LO>(numSameIDs); ++localRow) {
1362#ifdef HAVE_TPETRA_DEBUG
1363 if (!srcRowMap.isNodeLocalElement(localRow)) {
1364 lclErr = true;
1365 invalidSrcCopyRows.insert(localRow);
1366 continue; // skip invalid rows
1367 }
1368#endif // HAVE_TPETRA_DEBUG
1369
1370 local_inds_host_view_type lclSrcCols;
1371 values_host_view_type vals;
1372 LO numEntries;
1373 // If this call fails, that means the mesh row local index is
1374 // invalid. That means the Import or Export is invalid somehow.
1375 src->getLocalRowView(localRow, lclSrcCols, vals);
1376 numEntries = lclSrcCols.extent(0);
1377 if (numEntries > 0) {
1378 LO err = this->replaceLocalValues(localRow, lclSrcCols.data(), reinterpret_cast<const scalar_type*>(vals.data()), numEntries);
1379 if (err != numEntries) {
1380 lclErr = true;
1381 if (!dstRowMap.isNodeLocalElement(localRow)) {
1382#ifdef HAVE_TPETRA_DEBUG
1383 invalidDstCopyRows.insert(localRow);
1384#endif // HAVE_TPETRA_DEBUG
1385 } else {
1386 // Once there's an error, there's no sense in saving
1387 // time, so we check whether the column indices were
1388 // invalid. However, only remember which ones were
1389 // invalid in a debug build, because that might take a
1390 // lot of space.
1391 for (LO k = 0; k < numEntries; ++k) {
1392 if (!dstColMap.isNodeLocalElement(lclSrcCols[k])) {
1393 lclErr = true;
1394#ifdef HAVE_TPETRA_DEBUG
1395 (void)invalidDstCopyCols.insert(lclSrcCols[k]);
1396#endif // HAVE_TPETRA_DEBUG
1397 }
1398 }
1399 }
1400 }
1401 }
1402 } // for each "same" local row
1403
1404 // Copy the "permute" local rows.
1405 for (size_t k = 0; k < numPermute; ++k) {
1406 const LO srcLclRow = static_cast<LO>(permuteFromLIDsHost(k));
1407 const LO dstLclRow = static_cast<LO>(permuteToLIDsHost(k));
1408
1409 local_inds_host_view_type lclSrcCols;
1410 values_host_view_type vals;
1411 LO numEntries;
1412 src->getLocalRowView(srcLclRow, lclSrcCols, vals);
1413 numEntries = lclSrcCols.extent(0);
1414 if (numEntries > 0) {
1415 LO err = this->replaceLocalValues(dstLclRow, lclSrcCols.data(), reinterpret_cast<const scalar_type*>(vals.data()), numEntries);
1416 if (err != numEntries) {
1417 lclErr = true;
1418#ifdef HAVE_TPETRA_DEBUG
1419 for (LO c = 0; c < numEntries; ++c) {
1420 if (!dstColMap.isNodeLocalElement(lclSrcCols[c])) {
1421 invalidDstPermuteCols.insert(lclSrcCols[c]);
1422 }
1423 }
1424#endif // HAVE_TPETRA_DEBUG
1425 }
1426 }
1427 }
1428 } else { // must convert column indices to global
1429 // Reserve space to store the destination matrix's local column indices.
1430 const size_t maxNumEnt = src->graph_.getLocalMaxNumRowEntries();
1431 Teuchos::Array<LO> lclDstCols(maxNumEnt);
1432
1433 // Copy local rows that are the "same" in both source and target.
1434 for (LO localRow = 0; localRow < static_cast<LO>(numSameIDs); ++localRow) {
1435 local_inds_host_view_type lclSrcCols;
1436 values_host_view_type vals;
1437 LO numEntries;
1438
1439 // If this call fails, that means the mesh row local index is
1440 // invalid. That means the Import or Export is invalid somehow.
1441 try {
1442 src->getLocalRowView(localRow, lclSrcCols, vals);
1443 numEntries = lclSrcCols.extent(0);
1444 } catch (std::exception& e) {
1445 if (debug) {
1446 std::ostringstream os;
1447 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1448 os << "Proc " << myRank << ": copyAndPermute: At \"same\" localRow "
1449 << localRow << ", src->getLocalRowView() threw an exception: "
1450 << e.what();
1451 std::cerr << os.str();
1452 }
1453 throw e;
1454 }
1455
1456 if (numEntries > 0) {
1457 if (static_cast<size_t>(numEntries) > static_cast<size_t>(lclDstCols.size())) {
1458 lclErr = true;
1459 if (debug) {
1460 std::ostringstream os;
1461 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1462 os << "Proc " << myRank << ": copyAndPermute: At \"same\" localRow "
1463 << localRow << ", numEntries = " << numEntries << " > maxNumEnt = "
1464 << maxNumEnt << endl;
1465 std::cerr << os.str();
1466 }
1467 } else {
1468 // Convert the source matrix's local column indices to the
1469 // destination matrix's local column indices.
1470 Teuchos::ArrayView<LO> lclDstColsView = lclDstCols.view(0, numEntries);
1471 for (LO j = 0; j < numEntries; ++j) {
1472 lclDstColsView[j] = dstColMap.getLocalElement(srcColMap.getGlobalElement(lclSrcCols[j]));
1473 if (lclDstColsView[j] == Teuchos::OrdinalTraits<LO>::invalid()) {
1474 lclErr = true;
1475#ifdef HAVE_TPETRA_DEBUG
1476 invalidDstCopyCols.insert(lclDstColsView[j]);
1477#endif // HAVE_TPETRA_DEBUG
1478 }
1479 }
1480 LO err(0);
1481 try {
1482 err = this->replaceLocalValues(localRow, lclDstColsView.getRawPtr(), reinterpret_cast<const scalar_type*>(vals.data()), numEntries);
1483 } catch (std::exception& e) {
1484 if (debug) {
1485 std::ostringstream os;
1486 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1487 os << "Proc " << myRank << ": copyAndPermute: At \"same\" localRow "
1488 << localRow << ", this->replaceLocalValues() threw an exception: "
1489 << e.what();
1490 std::cerr << os.str();
1491 }
1492 throw e;
1493 }
1494 if (err != numEntries) {
1495 lclErr = true;
1496 if (debug) {
1497 std::ostringstream os;
1498 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1499 os << "Proc " << myRank << ": copyAndPermute: At \"same\" "
1500 "localRow "
1501 << localRow << ", this->replaceLocalValues "
1502 "returned "
1503 << err << " instead of numEntries = "
1504 << numEntries << endl;
1505 std::cerr << os.str();
1506 }
1507 }
1508 }
1509 }
1510 }
1511
1512 // Copy the "permute" local rows.
1513 for (size_t k = 0; k < numPermute; ++k) {
1514 const LO srcLclRow = static_cast<LO>(permuteFromLIDsHost(k));
1515 const LO dstLclRow = static_cast<LO>(permuteToLIDsHost(k));
1516
1517 local_inds_host_view_type lclSrcCols;
1518 values_host_view_type vals;
1519 LO numEntries;
1520
1521 try {
1522 src->getLocalRowView(srcLclRow, lclSrcCols, vals);
1523 numEntries = lclSrcCols.extent(0);
1524 } catch (std::exception& e) {
1525 if (debug) {
1526 std::ostringstream os;
1527 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1528 os << "Proc " << myRank << ": copyAndPermute: At \"permute\" "
1529 "srcLclRow "
1530 << srcLclRow << " and dstLclRow " << dstLclRow
1531 << ", src->getLocalRowView() threw an exception: " << e.what();
1532 std::cerr << os.str();
1533 }
1534 throw e;
1535 }
1536
1537 if (numEntries > 0) {
1538 if (static_cast<size_t>(numEntries) > static_cast<size_t>(lclDstCols.size())) {
1539 lclErr = true;
1540 } else {
1541 // Convert the source matrix's local column indices to the
1542 // destination matrix's local column indices.
1543 Teuchos::ArrayView<LO> lclDstColsView = lclDstCols.view(0, numEntries);
1544 for (LO j = 0; j < numEntries; ++j) {
1545 lclDstColsView[j] = dstColMap.getLocalElement(srcColMap.getGlobalElement(lclSrcCols[j]));
1546 if (lclDstColsView[j] == Teuchos::OrdinalTraits<LO>::invalid()) {
1547 lclErr = true;
1548#ifdef HAVE_TPETRA_DEBUG
1549 invalidDstPermuteCols.insert(lclDstColsView[j]);
1550#endif // HAVE_TPETRA_DEBUG
1551 }
1552 }
1553 LO err = this->replaceLocalValues(dstLclRow, lclDstColsView.getRawPtr(), reinterpret_cast<const scalar_type*>(vals.data()), numEntries);
1554 if (err != numEntries) {
1555 lclErr = true;
1556 }
1557 }
1558 }
1559 }
1560 }
1561
1562 if (lclErr) {
1563 std::ostream& err = this->markLocalErrorAndGetStream();
1564#ifdef HAVE_TPETRA_DEBUG
1565 err << "copyAndPermute: The graph structure of the source of the "
1566 "Import or Export must be a subset of the graph structure of the "
1567 "target. ";
1568 err << "invalidSrcCopyRows = [";
1569 for (typename std::set<LO>::const_iterator it = invalidSrcCopyRows.begin();
1570 it != invalidSrcCopyRows.end(); ++it) {
1571 err << *it;
1572 typename std::set<LO>::const_iterator itp1 = it;
1573 itp1++;
1574 if (itp1 != invalidSrcCopyRows.end()) {
1575 err << ",";
1576 }
1577 }
1578 err << "], invalidDstCopyRows = [";
1579 for (typename std::set<LO>::const_iterator it = invalidDstCopyRows.begin();
1580 it != invalidDstCopyRows.end(); ++it) {
1581 err << *it;
1582 typename std::set<LO>::const_iterator itp1 = it;
1583 itp1++;
1584 if (itp1 != invalidDstCopyRows.end()) {
1585 err << ",";
1586 }
1587 }
1588 err << "], invalidDstCopyCols = [";
1589 for (typename std::set<LO>::const_iterator it = invalidDstCopyCols.begin();
1590 it != invalidDstCopyCols.end(); ++it) {
1591 err << *it;
1592 typename std::set<LO>::const_iterator itp1 = it;
1593 itp1++;
1594 if (itp1 != invalidDstCopyCols.end()) {
1595 err << ",";
1596 }
1597 }
1598 err << "], invalidDstPermuteCols = [";
1599 for (typename std::set<LO>::const_iterator it = invalidDstPermuteCols.begin();
1600 it != invalidDstPermuteCols.end(); ++it) {
1601 err << *it;
1602 typename std::set<LO>::const_iterator itp1 = it;
1603 itp1++;
1604 if (itp1 != invalidDstPermuteCols.end()) {
1605 err << ",";
1606 }
1607 }
1608 err << "], invalidPermuteFromRows = [";
1609 for (typename std::set<LO>::const_iterator it = invalidPermuteFromRows.begin();
1610 it != invalidPermuteFromRows.end(); ++it) {
1611 err << *it;
1612 typename std::set<LO>::const_iterator itp1 = it;
1613 itp1++;
1614 if (itp1 != invalidPermuteFromRows.end()) {
1615 err << ",";
1616 }
1617 }
1618 err << "]" << endl;
1619#else
1620 err << "copyAndPermute: The graph structure of the source of the "
1621 "Import or Export must be a subset of the graph structure of the "
1622 "target."
1623 << endl;
1624#endif // HAVE_TPETRA_DEBUG
1625 }
1626
1627 if (debug) {
1628 std::ostringstream os;
1629 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1630 const bool lclSuccess = !(*(this->localError_));
1631 os << "*** Proc " << myRank << ": copyAndPermute "
1632 << (lclSuccess ? "succeeded" : "FAILED");
1633 if (lclSuccess) {
1634 os << endl;
1635 } else {
1636 os << ": error messages: " << this->errorMessages(); // comes w/ endl
1637 }
1638 std::cerr << os.str();
1639 }
1640}
1641
1642namespace { // (anonymous)
1643
1662template <class LO, class GO>
1663size_t
1664packRowCount(const size_t numEnt,
1665 const size_t numBytesPerValue,
1666 const size_t blkSize) {
1667 using ::Tpetra::Details::PackTraits;
1668
1669 if (numEnt == 0) {
1670 // Empty rows always take zero bytes, to ensure sparsity.
1671 return 0;
1672 } else {
1673 // We store the number of entries as a local index (LO).
1674 LO numEntLO = 0; // packValueCount wants this.
1675 GO gid{};
1676 const size_t numEntLen = PackTraits<LO>::packValueCount(numEntLO);
1677 const size_t gidsLen = numEnt * PackTraits<GO>::packValueCount(gid);
1678 const size_t valsLen = numEnt * numBytesPerValue * blkSize * blkSize;
1679 return numEntLen + gidsLen + valsLen;
1680 }
1681}
1682
1693template <class ST, class LO, class GO>
1694size_t
1695unpackRowCount(const typename ::Tpetra::Details::PackTraits<LO>::input_buffer_type& imports,
1696 const size_t offset,
1697 const size_t numBytes,
1698 const size_t /* numBytesPerValue */) {
1699 using Kokkos::subview;
1700 using ::Tpetra::Details::PackTraits;
1701
1702 if (numBytes == 0) {
1703 // Empty rows always take zero bytes, to ensure sparsity.
1704 return static_cast<size_t>(0);
1705 } else {
1706 LO numEntLO = 0;
1707 const size_t theNumBytes = PackTraits<LO>::packValueCount(numEntLO);
1708 TEUCHOS_TEST_FOR_EXCEPTION(theNumBytes > numBytes, std::logic_error,
1709 "unpackRowCount: "
1710 "theNumBytes = "
1711 << theNumBytes << " < numBytes = " << numBytes
1712 << ".");
1713 const char* const inBuf = imports.data() + offset;
1714 const size_t actualNumBytes = PackTraits<LO>::unpackValue(numEntLO, inBuf);
1715 TEUCHOS_TEST_FOR_EXCEPTION(actualNumBytes > numBytes, std::logic_error,
1716 "unpackRowCount: "
1717 "actualNumBytes = "
1718 << actualNumBytes << " < numBytes = " << numBytes
1719 << ".");
1720 return static_cast<size_t>(numEntLO);
1721 }
1722}
1723
1727template <class ST, class LO, class GO>
1728size_t
1729packRowForBlockCrs(const typename ::Tpetra::Details::PackTraits<LO>::output_buffer_type exports,
1730 const size_t offset,
1731 const size_t numEnt,
1732 const typename ::Tpetra::Details::PackTraits<GO>::input_array_type& gidsIn,
1733 const typename ::Tpetra::Details::PackTraits<ST>::input_array_type& valsIn,
1734 const size_t numBytesPerValue,
1735 const size_t blockSize) {
1736 using Kokkos::subview;
1737 using ::Tpetra::Details::PackTraits;
1738
1739 if (numEnt == 0) {
1740 // Empty rows always take zero bytes, to ensure sparsity.
1741 return 0;
1742 }
1743 const size_t numScalarEnt = numEnt * blockSize * blockSize;
1744
1745 const GO gid = 0; // packValueCount wants this
1746 const LO numEntLO = static_cast<size_t>(numEnt);
1747
1748 const size_t numEntBeg = offset;
1749 const size_t numEntLen = PackTraits<LO>::packValueCount(numEntLO);
1750 const size_t gidsBeg = numEntBeg + numEntLen;
1751 const size_t gidsLen = numEnt * PackTraits<GO>::packValueCount(gid);
1752 const size_t valsBeg = gidsBeg + gidsLen;
1753 const size_t valsLen = numScalarEnt * numBytesPerValue;
1754
1755 char* const numEntOut = exports.data() + numEntBeg;
1756 char* const gidsOut = exports.data() + gidsBeg;
1757 char* const valsOut = exports.data() + valsBeg;
1758
1759 size_t numBytesOut = 0;
1760 int errorCode = 0;
1761 numBytesOut += PackTraits<LO>::packValue(numEntOut, numEntLO);
1762
1763 {
1764 Kokkos::pair<int, size_t> p;
1765 p = PackTraits<GO>::packArray(gidsOut, gidsIn.data(), numEnt);
1766 errorCode += p.first;
1767 numBytesOut += p.second;
1768
1769 p = PackTraits<ST>::packArray(valsOut, valsIn.data(), numScalarEnt);
1770 errorCode += p.first;
1771 numBytesOut += p.second;
1772 }
1773
1774 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
1775 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != expectedNumBytes, std::logic_error,
1776 "packRowForBlockCrs: numBytesOut = " << numBytesOut
1777 << " != expectedNumBytes = " << expectedNumBytes << ".");
1778
1779 TEUCHOS_TEST_FOR_EXCEPTION(errorCode != 0, std::runtime_error,
1780 "packRowForBlockCrs: "
1781 "PackTraits::packArray returned a nonzero error code");
1782
1783 return numBytesOut;
1784}
1785
1786// Return the number of bytes actually read / used.
1787template <class ST, class LO, class GO>
1788size_t
1789unpackRowForBlockCrs(const typename ::Tpetra::Details::PackTraits<GO>::output_array_type& gidsOut,
1790 const typename ::Tpetra::Details::PackTraits<ST>::output_array_type& valsOut,
1791 const typename ::Tpetra::Details::PackTraits<int>::input_buffer_type& imports,
1792 const size_t offset,
1793 const size_t numBytes,
1794 const size_t numEnt,
1795 const size_t numBytesPerValue,
1796 const size_t blockSize) {
1797 using ::Tpetra::Details::PackTraits;
1798
1799 if (numBytes == 0) {
1800 // Rows with zero bytes always have zero entries.
1801 return 0;
1802 }
1803 const size_t numScalarEnt = numEnt * blockSize * blockSize;
1804 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(imports.extent(0)) <= offset,
1805 std::logic_error, "unpackRowForBlockCrs: imports.extent(0) = " << imports.extent(0) << " <= offset = " << offset << ".");
1806 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(imports.extent(0)) < offset + numBytes,
1807 std::logic_error, "unpackRowForBlockCrs: imports.extent(0) = " << imports.extent(0) << " < offset + numBytes = " << (offset + numBytes) << ".");
1808
1809 const GO gid = 0; // packValueCount wants this
1810 const LO lid = 0; // packValueCount wants this
1811
1812 const size_t numEntBeg = offset;
1813 const size_t numEntLen = PackTraits<LO>::packValueCount(lid);
1814 const size_t gidsBeg = numEntBeg + numEntLen;
1815 const size_t gidsLen = numEnt * PackTraits<GO>::packValueCount(gid);
1816 const size_t valsBeg = gidsBeg + gidsLen;
1817 const size_t valsLen = numScalarEnt * numBytesPerValue;
1818
1819 const char* const numEntIn = imports.data() + numEntBeg;
1820 const char* const gidsIn = imports.data() + gidsBeg;
1821 const char* const valsIn = imports.data() + valsBeg;
1822
1823 size_t numBytesOut = 0;
1824 int errorCode = 0;
1825 LO numEntOut;
1826 numBytesOut += PackTraits<LO>::unpackValue(numEntOut, numEntIn);
1827 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(numEntOut) != numEnt, std::logic_error,
1828 "unpackRowForBlockCrs: Expected number of entries " << numEnt
1829 << " != actual number of entries " << numEntOut << ".");
1830
1831 {
1832 Kokkos::pair<int, size_t> p;
1833 p = PackTraits<GO>::unpackArray(gidsOut.data(), gidsIn, numEnt);
1834 errorCode += p.first;
1835 numBytesOut += p.second;
1836
1837 p = PackTraits<ST>::unpackArray(valsOut.data(), valsIn, numScalarEnt);
1838 errorCode += p.first;
1839 numBytesOut += p.second;
1840 }
1841
1842 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != numBytes, std::logic_error,
1843 "unpackRowForBlockCrs: numBytesOut = " << numBytesOut
1844 << " != numBytes = " << numBytes << ".");
1845
1846 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
1847 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != expectedNumBytes, std::logic_error,
1848 "unpackRowForBlockCrs: numBytesOut = " << numBytesOut
1849 << " != expectedNumBytes = " << expectedNumBytes << ".");
1850
1851 TEUCHOS_TEST_FOR_EXCEPTION(errorCode != 0, std::runtime_error,
1852 "unpackRowForBlockCrs: "
1853 "PackTraits::unpackArray returned a nonzero error code");
1854
1855 return numBytesOut;
1856}
1857} // namespace
1858
1859template <class Scalar, class LO, class GO, class Node>
1860void BlockCrsMatrix<Scalar, LO, GO, Node>::
1861 packAndPrepare(const ::Tpetra::SrcDistObject& source,
1862 const Kokkos::DualView<const local_ordinal_type*,
1863 buffer_device_type>& exportLIDs,
1864 Kokkos::DualView<packet_type*,
1865 buffer_device_type>& exports, // output
1866 Kokkos::DualView<size_t*,
1867 buffer_device_type>
1868 numPacketsPerLID, // output
1869 size_t& constantNumPackets) {
1870 using ::Tpetra::Details::Behavior;
1871 using ::Tpetra::Details::dualViewStatusToString;
1872 using ::Tpetra::Details::PackTraits;
1873 using ::Tpetra::Details::ProfilingRegion;
1874
1875 typedef typename Kokkos::View<int*, device_type>::host_mirror_type::execution_space host_exec;
1876
1877 typedef BlockCrsMatrix<Scalar, LO, GO, Node> this_BCRS_type;
1878
1879 ProfilingRegion profile_region("Tpetra::BlockCrsMatrix::packAndPrepare");
1880
1881 const bool debug = Behavior::debug();
1882 const bool verbose = Behavior::verbose();
1883
1884 // Define this function prefix
1885 std::string prefix;
1886 {
1887 std::ostringstream os;
1888 const int myRank = this->graph_.getRowMap()->getComm()->getRank();
1889 os << "Proc " << myRank << ": BlockCrsMatrix::packAndPrepare: " << std::endl;
1890 prefix = os.str();
1891 }
1892
1893 // check if this already includes a local error
1894 if (*(this->localError_)) {
1895 std::ostream& err = this->markLocalErrorAndGetStream();
1896 err << prefix
1897 << "The target object of the Import or Export is already in an error state."
1898 << std::endl;
1899 return;
1900 }
1901
1902 //
1903 // Verbose input dual view status
1904 //
1905 if (verbose) {
1906 std::ostringstream os;
1907 os << prefix << std::endl
1908 << prefix << " " << dualViewStatusToString(exportLIDs, "exportLIDs") << std::endl
1909 << prefix << " " << dualViewStatusToString(exports, "exports") << std::endl
1910 << prefix << " " << dualViewStatusToString(numPacketsPerLID, "numPacketsPerLID") << std::endl;
1911 std::cerr << os.str();
1912 }
1913
1917 if (exportLIDs.extent(0) != numPacketsPerLID.extent(0)) {
1918 std::ostream& err = this->markLocalErrorAndGetStream();
1919 err << prefix
1920 << "exportLIDs.extent(0) = " << exportLIDs.extent(0)
1921 << " != numPacketsPerLID.extent(0) = " << numPacketsPerLID.extent(0)
1922 << "." << std::endl;
1923 return;
1924 }
1925 if (exportLIDs.need_sync_host()) {
1926 std::ostream& err = this->markLocalErrorAndGetStream();
1927 err << prefix << "exportLIDs be sync'd to host." << std::endl;
1928 return;
1929 }
1930
1931 const this_BCRS_type* src = dynamic_cast<const this_BCRS_type*>(&source);
1932 if (src == NULL) {
1933 std::ostream& err = this->markLocalErrorAndGetStream();
1934 err << prefix
1935 << "The source (input) object of the Import or "
1936 "Export is either not a BlockCrsMatrix, or does not have the right "
1937 "template parameters. checkSizes() should have caught this. "
1938 "Please report this bug to the Tpetra developers."
1939 << std::endl;
1940 return;
1941 }
1942
1943 // Graphs and matrices are allowed to have a variable number of
1944 // entries per row. We could test whether all rows have the same
1945 // number of entries, but DistObject can only use this
1946 // optimization if all rows on _all_ processes have the same
1947 // number of entries. Rather than do the all-reduce necessary to
1948 // test for this unlikely case, we tell DistObject (by setting
1949 // constantNumPackets to zero) to assume that different rows may
1950 // have different numbers of entries.
1951 constantNumPackets = 0;
1952
1953 // const values
1954 const crs_graph_type& srcGraph = src->graph_;
1955 const size_t blockSize = static_cast<size_t>(src->getBlockSize());
1956 const size_t numExportLIDs = exportLIDs.extent(0);
1957 size_t numBytesPerValue(0);
1958 {
1959 auto val_host = val_.getHostView(Access::ReadOnly);
1960 numBytesPerValue =
1961 PackTraits<impl_scalar_type>::packValueCount(val_host.extent(0) ? val_host(0) : impl_scalar_type());
1962 }
1963
1964 // Compute the number of bytes ("packets") per row to pack. While
1965 // we're at it, compute the total # of block entries to send, and
1966 // the max # of block entries in any of the rows we're sending.
1967
1968 Impl::BlockCrsRowStruct<size_t> rowReducerStruct;
1969
1970 // Graph information is on host; let's do this on host parallel reduce
1971 auto exportLIDsHost = exportLIDs.view_host();
1972 auto numPacketsPerLIDHost = numPacketsPerLID.view_host(); // we will modify this.
1973 numPacketsPerLID.modify_host();
1974 {
1975 rowReducerStruct = Impl::BlockCrsRowStruct<size_t>();
1976 for (size_t i = 0; i < numExportLIDs; ++i) {
1977 const LO lclRow = exportLIDsHost(i);
1978 size_t numEnt = srcGraph.getNumEntriesInLocalRow(lclRow);
1979 numEnt = (numEnt == Teuchos::OrdinalTraits<size_t>::invalid() ? 0 : numEnt);
1980
1981 const size_t numBytes =
1982 packRowCount<LO, GO>(numEnt, numBytesPerValue, blockSize);
1983 numPacketsPerLIDHost(i) = numBytes;
1984 rowReducerStruct += Impl::BlockCrsRowStruct<size_t>(numEnt, numBytes, numEnt);
1985 }
1986 }
1987
1988 // Compute the number of bytes ("packets") per row to pack. While
1989 // we're at it, compute the total # of block entries to send, and
1990 // the max # of block entries in any of the rows we're sending.
1991 const size_t totalNumBytes = rowReducerStruct.totalNumBytes;
1992 const size_t totalNumEntries = rowReducerStruct.totalNumEntries;
1993 const size_t maxRowLength = rowReducerStruct.maxRowLength;
1994
1995 if (verbose) {
1996 std::ostringstream os;
1997 os << prefix
1998 << "totalNumBytes = " << totalNumBytes << ", totalNumEntries = " << totalNumEntries
1999 << std::endl;
2000 std::cerr << os.str();
2001 }
2002
2003 // We use a "struct of arrays" approach to packing each row's
2004 // entries. All the column indices (as global indices) go first,
2005 // then all their owning process ranks, and then the values.
2006 if (exports.extent(0) != totalNumBytes) {
2007 const std::string oldLabel = exports.view_device().label();
2008 const std::string newLabel = (oldLabel == "") ? "exports" : oldLabel;
2009 exports = Kokkos::DualView<packet_type*, buffer_device_type>(newLabel, totalNumBytes);
2010 }
2011 if (totalNumEntries > 0) {
2012 // Current position (in bytes) in the 'exports' output array.
2013 Kokkos::View<size_t*, host_exec> offset("offset", numExportLIDs + 1);
2014 {
2015 const auto policy = Kokkos::RangePolicy<host_exec>(size_t(0), numExportLIDs + 1);
2016 Kokkos::parallel_scan(policy,
2017 [=](const size_t& i, size_t& update, const bool& final) {
2018 if (final) offset(i) = update;
2019 update += (i == numExportLIDs ? 0 : numPacketsPerLIDHost(i));
2020 });
2021 }
2022 if (offset(numExportLIDs) != totalNumBytes) {
2023 std::ostream& err = this->markLocalErrorAndGetStream();
2024 err << prefix
2025 << "At end of method, the final offset (in bytes) "
2026 << offset(numExportLIDs) << " does not equal the total number of bytes packed "
2027 << totalNumBytes << ". "
2028 << "Please report this bug to the Tpetra developers." << std::endl;
2029 return;
2030 }
2031
2032 // For each block row of the matrix owned by the calling
2033 // process, pack that block row's column indices and values into
2034 // the exports array.
2035
2036 // Source matrix's column Map. We verified in checkSizes() that
2037 // the column Map exists (is not null).
2038 const map_type& srcColMap = *(srcGraph.getColMap());
2039
2040 // Pack the data for each row to send, into the 'exports' buffer.
2041 // exports will be modified on host.
2042 exports.modify_host();
2043 {
2044 typedef Kokkos::TeamPolicy<host_exec> policy_type;
2045 const auto policy =
2046 policy_type(numExportLIDs, 1, 1)
2047 .set_scratch_size(0, Kokkos::PerTeam(sizeof(GO) * maxRowLength));
2048 // The following parallel_for needs const access to the local values of src.
2049 // (the local graph is also accessed on host, but this does not use WDVs).
2050 getValuesHost();
2051 Details::disableWDVTracking();
2052 Kokkos::parallel_for(policy,
2053 [=](const typename policy_type::member_type& member) {
2054 const size_t i = member.league_rank();
2055 Kokkos::View<GO*, typename host_exec::scratch_memory_space>
2056 gblColInds(member.team_scratch(0), maxRowLength);
2057
2058 const LO lclRowInd = exportLIDsHost(i);
2059 local_inds_host_view_type lclColInds;
2060 values_host_view_type vals;
2061
2062 // It's OK to ignore the return value, since if the calling
2063 // process doesn't own that local row, then the number of
2064 // entries in that row on the calling process is zero.
2065 src->getLocalRowView(lclRowInd, lclColInds, vals);
2066 const size_t numEnt = lclColInds.extent(0);
2067
2068 // Convert column indices from local to global.
2069 for (size_t j = 0; j < numEnt; ++j)
2070 gblColInds(j) = srcColMap.getGlobalElement(lclColInds(j));
2071
2072 // Kyungjoo: additional wrapping scratch view is necessary
2073 // the following function interface need the same execution space
2074 // host scratch space somehow is not considered same as the host_exec
2075 // Copy the row's data into the current spot in the exports array.
2076 const size_t numBytes =
2077 packRowForBlockCrs<impl_scalar_type, LO, GO>(exports.view_host(),
2078 offset(i),
2079 numEnt,
2080 Kokkos::View<const GO*, host_exec>(gblColInds.data(), maxRowLength),
2081 vals,
2082 numBytesPerValue,
2083 blockSize);
2084
2085 // numBytes should be same as the difference between offsets
2086 if (debug) {
2087 const size_t offsetDiff = offset(i + 1) - offset(i);
2088 if (numBytes != offsetDiff) {
2089 std::ostringstream os;
2090 os << prefix
2091 << "numBytes computed from packRowForBlockCrs is different from "
2092 << "precomputed offset values, LID = " << i << std::endl;
2093 std::cerr << os.str();
2094 }
2095 }
2096 }); // for each LID (of a row) to send
2097 Details::enableWDVTracking();
2098 }
2099 } // if totalNumEntries > 0
2100
2101 if (debug) {
2102 std::ostringstream os;
2103 const bool lclSuccess = !(*(this->localError_));
2104 os << prefix
2105 << (lclSuccess ? "succeeded" : "FAILED")
2106 << std::endl;
2107 std::cerr << os.str();
2108 }
2109}
2110
2111template <class Scalar, class LO, class GO, class Node>
2112void BlockCrsMatrix<Scalar, LO, GO, Node>::
2113 unpackAndCombine(const Kokkos::DualView<const local_ordinal_type*,
2114 buffer_device_type>& importLIDs,
2115 Kokkos::DualView<packet_type*,
2116 buffer_device_type>
2117 imports,
2118 Kokkos::DualView<size_t*,
2119 buffer_device_type>
2120 numPacketsPerLID,
2121 const size_t /* constantNumPackets */,
2122 const CombineMode combineMode) {
2123 using std::endl;
2124 using ::Tpetra::Details::Behavior;
2125 using ::Tpetra::Details::dualViewStatusToString;
2126 using ::Tpetra::Details::PackTraits;
2127 using ::Tpetra::Details::ProfilingRegion;
2128 using host_exec =
2129 typename Kokkos::View<int*, device_type>::host_mirror_type::execution_space;
2130
2131 ProfilingRegion profile_region("Tpetra::BlockCrsMatrix::unpackAndCombine");
2132 const bool verbose = Behavior::verbose();
2133
2134 // Define this function prefix
2135 std::string prefix;
2136 {
2137 std::ostringstream os;
2138 auto map = this->graph_.getRowMap();
2139 auto comm = map.is_null() ? Teuchos::null : map->getComm();
2140 const int myRank = comm.is_null() ? -1 : comm->getRank();
2141 os << "Proc " << myRank << ": Tpetra::BlockCrsMatrix::unpackAndCombine: ";
2142 prefix = os.str();
2143 if (verbose) {
2144 os << "Start" << endl;
2145 std::cerr << os.str();
2146 }
2147 }
2148
2149 // check if this already includes a local error
2150 if (*(this->localError_)) {
2151 std::ostream& err = this->markLocalErrorAndGetStream();
2152 std::ostringstream os;
2153 os << prefix << "{Im/Ex}port target is already in error." << endl;
2154 if (verbose) {
2155 std::cerr << os.str();
2156 }
2157 err << os.str();
2158 return;
2159 }
2160
2164 if (importLIDs.extent(0) != numPacketsPerLID.extent(0)) {
2165 std::ostream& err = this->markLocalErrorAndGetStream();
2166 std::ostringstream os;
2167 os << prefix << "importLIDs.extent(0) = " << importLIDs.extent(0)
2168 << " != numPacketsPerLID.extent(0) = " << numPacketsPerLID.extent(0)
2169 << "." << endl;
2170 if (verbose) {
2171 std::cerr << os.str();
2172 }
2173 err << os.str();
2174 return;
2175 }
2176
2177 if (combineMode != ADD && combineMode != INSERT &&
2178 combineMode != REPLACE && combineMode != ABSMAX &&
2179 combineMode != ZERO) {
2180 std::ostream& err = this->markLocalErrorAndGetStream();
2181 std::ostringstream os;
2182 os << prefix
2183 << "Invalid CombineMode value " << combineMode << ". Valid "
2184 << "values include ADD, INSERT, REPLACE, ABSMAX, and ZERO."
2185 << std::endl;
2186 if (verbose) {
2187 std::cerr << os.str();
2188 }
2189 err << os.str();
2190 return;
2191 }
2192
2193 if (this->graph_.getColMap().is_null()) {
2194 std::ostream& err = this->markLocalErrorAndGetStream();
2195 std::ostringstream os;
2196 os << prefix << "Target matrix's column Map is null." << endl;
2197 if (verbose) {
2198 std::cerr << os.str();
2199 }
2200 err << os.str();
2201 return;
2202 }
2203
2204 // Target matrix's column Map. Use to convert the global column
2205 // indices in the receive buffer to local indices. We verified in
2206 // checkSizes() that the column Map exists (is not null).
2207 const map_type& tgtColMap = *(this->graph_.getColMap());
2208
2209 // Const values
2210 const size_t blockSize = this->getBlockSize();
2211 const size_t numImportLIDs = importLIDs.extent(0);
2212 // FIXME (mfh 06 Feb 2019) For scalar types with run-time sizes, a
2213 // default-constructed instance could have a different size than
2214 // other instances. (We assume all nominally constructed
2215 // instances have the same size; that's not the issue here.) This
2216 // could be bad if the calling process has no entries, but other
2217 // processes have entries that they want to send to us.
2218 size_t numBytesPerValue(0);
2219 {
2220 auto val_host = val_.getHostView(Access::ReadOnly);
2221 numBytesPerValue =
2222 PackTraits<impl_scalar_type>::packValueCount(val_host.extent(0) ? val_host(0) : impl_scalar_type());
2223 }
2224 const size_t maxRowNumEnt = graph_.getLocalMaxNumRowEntries();
2225 const size_t maxRowNumScalarEnt = maxRowNumEnt * blockSize * blockSize;
2226
2227 if (verbose) {
2228 std::ostringstream os;
2229 os << prefix << "combineMode: "
2230 << ::Tpetra::combineModeToString(combineMode)
2231 << ", blockSize: " << blockSize
2232 << ", numImportLIDs: " << numImportLIDs
2233 << ", numBytesPerValue: " << numBytesPerValue
2234 << ", maxRowNumEnt: " << maxRowNumEnt
2235 << ", maxRowNumScalarEnt: " << maxRowNumScalarEnt << endl;
2236 std::cerr << os.str();
2237 }
2238
2239 if (combineMode == ZERO || numImportLIDs == 0) {
2240 if (verbose) {
2241 std::ostringstream os;
2242 os << prefix << "Nothing to unpack. Done!" << std::endl;
2243 std::cerr << os.str();
2244 }
2245 return;
2246 }
2247
2248 // NOTE (mfh 07 Feb 2019) If we ever implement unpack on device,
2249 // we can remove this sync.
2250 if (imports.need_sync_host()) {
2251 imports.sync_host();
2252 }
2253
2254 // NOTE (mfh 07 Feb 2019) DistObject::doTransferNew has already
2255 // sync'd numPacketsPerLID to host, since it needs to do that in
2256 // order to post MPI_Irecv messages with the correct lengths. A
2257 // hypothetical device-specific boundary exchange implementation
2258 // could possibly receive data without sync'ing lengths to host,
2259 // but we don't need to design for that nonexistent thing yet.
2260
2261 if (imports.need_sync_host() || numPacketsPerLID.need_sync_host() ||
2262 importLIDs.need_sync_host()) {
2263 std::ostream& err = this->markLocalErrorAndGetStream();
2264 std::ostringstream os;
2265 os << prefix << "All input DualViews must be sync'd to host by now. "
2266 << "imports_nc.need_sync_host()="
2267 << (imports.need_sync_host() ? "true" : "false")
2268 << ", numPacketsPerLID.need_sync_host()="
2269 << (numPacketsPerLID.need_sync_host() ? "true" : "false")
2270 << ", importLIDs.need_sync_host()="
2271 << (importLIDs.need_sync_host() ? "true" : "false")
2272 << "." << endl;
2273 if (verbose) {
2274 std::cerr << os.str();
2275 }
2276 err << os.str();
2277 return;
2278 }
2279
2280 const auto importLIDsHost = importLIDs.view_host();
2281 const auto numPacketsPerLIDHost = numPacketsPerLID.view_host();
2282
2283 // FIXME (mfh 06 Feb 2019) We could fuse the scan with the unpack
2284 // loop, by only unpacking on final==true (when we know the
2285 // current offset's value).
2286
2287 Kokkos::View<size_t*, host_exec> offset("offset", numImportLIDs + 1);
2288 {
2289 const auto policy = Kokkos::RangePolicy<host_exec>(size_t(0), numImportLIDs + 1);
2290 Kokkos::parallel_scan("Tpetra::BlockCrsMatrix::unpackAndCombine: offsets", policy,
2291 [=](const size_t& i, size_t& update, const bool& final) {
2292 if (final) offset(i) = update;
2293 update += (i == numImportLIDs ? 0 : numPacketsPerLIDHost(i));
2294 });
2295 }
2296
2297 // this variable does not matter with a race condition (just error flag)
2298 //
2299 // NOTE (mfh 06 Feb 2019) CUDA doesn't necessarily like reductions
2300 // or atomics on bool, so we use int instead. (I know we're not
2301 // launching a CUDA loop here, but we could in the future, and I'd
2302 // like to avoid potential trouble.)
2303 Kokkos::View<int, host_exec, Kokkos::MemoryTraits<Kokkos::Atomic> >
2304 errorDuringUnpack("errorDuringUnpack");
2305 errorDuringUnpack() = 0;
2306 {
2307 using policy_type = Kokkos::TeamPolicy<host_exec>;
2308 size_t scratch_per_row = sizeof(GO) * maxRowNumEnt + sizeof(LO) * maxRowNumEnt + numBytesPerValue * maxRowNumScalarEnt + 2 * sizeof(GO); // Yeah, this is a fudge factor
2309
2310 const auto policy = policy_type(numImportLIDs, 1, 1)
2311 .set_scratch_size(0, Kokkos::PerTeam(scratch_per_row));
2312 using host_scratch_space = typename host_exec::scratch_memory_space;
2313
2314 using pair_type = Kokkos::pair<size_t, size_t>;
2315
2316 // The following parallel_for modifies values on host while unpacking.
2317 getValuesHostNonConst();
2318 Details::disableWDVTracking();
2319 Kokkos::parallel_for("Tpetra::BlockCrsMatrix::unpackAndCombine: unpack", policy,
2320 [=, *this](const typename policy_type::member_type& member) {
2321 const size_t i = member.league_rank();
2322 Kokkos::View<GO*, host_scratch_space> gblColInds(member.team_scratch(0), maxRowNumEnt);
2323 Kokkos::View<LO*, host_scratch_space> lclColInds(member.team_scratch(0), maxRowNumEnt);
2324 Kokkos::View<impl_scalar_type*, host_scratch_space> vals(member.team_scratch(0), maxRowNumScalarEnt);
2325
2326 const size_t offval = offset(i);
2327 const LO lclRow = importLIDsHost(i);
2328 const size_t numBytes = numPacketsPerLIDHost(i);
2329 const size_t numEnt =
2330 unpackRowCount<impl_scalar_type, LO, GO>(imports.view_host(), offval, numBytes, numBytesPerValue);
2331
2332 if (numBytes > 0) {
2333 if (numEnt > maxRowNumEnt) {
2334 errorDuringUnpack() = 1;
2335 if (verbose) {
2336 std::ostringstream os;
2337 os << prefix
2338 << "At i = " << i << ", numEnt = " << numEnt
2339 << " > maxRowNumEnt = " << maxRowNumEnt
2340 << std::endl;
2341 std::cerr << os.str();
2342 }
2343 return;
2344 }
2345 }
2346 const size_t numScalarEnt = numEnt * blockSize * blockSize;
2347 auto gidsOut = Kokkos::subview(gblColInds, pair_type(0, numEnt));
2348 auto lidsOut = Kokkos::subview(lclColInds, pair_type(0, numEnt));
2349 auto valsOut = Kokkos::subview(vals, pair_type(0, numScalarEnt));
2350
2351 // Kyungjoo: additional wrapping scratch view is necessary
2352 // the following function interface need the same execution space
2353 // host scratch space somehow is not considered same as the host_exec
2354 size_t numBytesOut = 0;
2355 try {
2356 numBytesOut =
2357 unpackRowForBlockCrs<impl_scalar_type, LO, GO>(Kokkos::View<GO*, host_exec>(gidsOut.data(), numEnt),
2358 Kokkos::View<impl_scalar_type*, host_exec>(valsOut.data(), numScalarEnt),
2359 imports.view_host(),
2360 offval, numBytes, numEnt,
2361 numBytesPerValue, blockSize);
2362 } catch (std::exception& e) {
2363 errorDuringUnpack() = 1;
2364 if (verbose) {
2365 std::ostringstream os;
2366 os << prefix << "At i = " << i << ", unpackRowForBlockCrs threw: "
2367 << e.what() << endl;
2368 std::cerr << os.str();
2369 }
2370 return;
2371 }
2372
2373 if (numBytes != numBytesOut) {
2374 errorDuringUnpack() = 1;
2375 if (verbose) {
2376 std::ostringstream os;
2377 os << prefix
2378 << "At i = " << i << ", numBytes = " << numBytes
2379 << " != numBytesOut = " << numBytesOut << "."
2380 << std::endl;
2381 std::cerr << os.str();
2382 }
2383 return;
2384 }
2385
2386 // Convert incoming global indices to local indices.
2387 for (size_t k = 0; k < numEnt; ++k) {
2388 lidsOut(k) = tgtColMap.getLocalElement(gidsOut(k));
2389 if (lidsOut(k) == Teuchos::OrdinalTraits<LO>::invalid()) {
2390 if (verbose) {
2391 std::ostringstream os;
2392 os << prefix
2393 << "At i = " << i << ", GID " << gidsOut(k)
2394 << " is not owned by the calling process."
2395 << std::endl;
2396 std::cerr << os.str();
2397 }
2398 return;
2399 }
2400 }
2401
2402 // Combine the incoming data with the matrix's current data.
2403 LO numCombd = 0;
2404 const LO* const lidsRaw = const_cast<const LO*>(lidsOut.data());
2405 const Scalar* const valsRaw = reinterpret_cast<const Scalar*>(const_cast<const impl_scalar_type*>(valsOut.data()));
2406 if (combineMode == ADD) {
2407 numCombd = this->sumIntoLocalValues(lclRow, lidsRaw, valsRaw, numEnt);
2408 } else if (combineMode == INSERT || combineMode == REPLACE) {
2409 numCombd = this->replaceLocalValues(lclRow, lidsRaw, valsRaw, numEnt);
2410 } else if (combineMode == ABSMAX) {
2411 numCombd = this->absMaxLocalValues(lclRow, lidsRaw, valsRaw, numEnt);
2412 }
2413
2414 if (static_cast<LO>(numEnt) != numCombd) {
2415 errorDuringUnpack() = 1;
2416 if (verbose) {
2417 std::ostringstream os;
2418 os << prefix << "At i = " << i << ", numEnt = " << numEnt
2419 << " != numCombd = " << numCombd << "."
2420 << endl;
2421 std::cerr << os.str();
2422 }
2423 return;
2424 }
2425 }); // for each import LID i
2426 Details::enableWDVTracking();
2427 }
2428
2429 if (errorDuringUnpack() != 0) {
2430 std::ostream& err = this->markLocalErrorAndGetStream();
2431 err << prefix << "Unpacking failed.";
2432 if (!verbose) {
2433 err << " Please run again with the environment variable setting "
2434 "TPETRA_VERBOSE=1 to get more verbose diagnostic output.";
2435 }
2436 err << endl;
2437 }
2438
2439 if (verbose) {
2440 std::ostringstream os;
2441 const bool lclSuccess = !(*(this->localError_));
2442 os << prefix
2443 << (lclSuccess ? "succeeded" : "FAILED")
2444 << std::endl;
2445 std::cerr << os.str();
2446 }
2447}
2448
2449template <class Scalar, class LO, class GO, class Node>
2450std::string
2451BlockCrsMatrix<Scalar, LO, GO, Node>::description() const {
2452 using Teuchos::TypeNameTraits;
2453 std::ostringstream os;
2454 os << "\"Tpetra::BlockCrsMatrix\": { "
2455 << "Template parameters: { "
2456 << "Scalar: " << TypeNameTraits<Scalar>::name()
2457 << "LO: " << TypeNameTraits<LO>::name()
2458 << "GO: " << TypeNameTraits<GO>::name()
2459 << "Node: " << TypeNameTraits<Node>::name()
2460 << " }"
2461 << ", Label: \"" << this->getObjectLabel() << "\""
2462 << ", Global dimensions: ["
2463 << graph_.getDomainMap()->getGlobalNumElements() << ", "
2464 << graph_.getRangeMap()->getGlobalNumElements() << "]"
2465 << ", Block size: " << getBlockSize()
2466 << " }";
2467 return os.str();
2468}
2469
2470template <class Scalar, class LO, class GO, class Node>
2471void BlockCrsMatrix<Scalar, LO, GO, Node>::
2472 describe(Teuchos::FancyOStream& out,
2473 const Teuchos::EVerbosityLevel verbLevel) const {
2474 using Teuchos::ArrayRCP;
2475 using Teuchos::CommRequest;
2476 using Teuchos::FancyOStream;
2477 using Teuchos::getFancyOStream;
2478 using Teuchos::ireceive;
2479 using Teuchos::isend;
2480 using Teuchos::outArg;
2481 using Teuchos::TypeNameTraits;
2482 using Teuchos::VERB_DEFAULT;
2483 using Teuchos::VERB_LOW;
2484 using Teuchos::VERB_MEDIUM;
2485 using Teuchos::VERB_NONE;
2486 // using Teuchos::VERB_HIGH;
2487 using std::endl;
2488 using Teuchos::RCP;
2489 using Teuchos::VERB_EXTREME;
2490 using Teuchos::wait;
2491
2492 // Set default verbosity if applicable.
2493 const Teuchos::EVerbosityLevel vl =
2494 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
2495
2496 if (vl == VERB_NONE) {
2497 return; // print nothing
2498 }
2499
2500 // describe() always starts with a tab before it prints anything.
2501 Teuchos::OSTab tab0(out);
2502
2503 out << "\"Tpetra::BlockCrsMatrix\":" << endl;
2504 Teuchos::OSTab tab1(out);
2505
2506 out << "Template parameters:" << endl;
2507 {
2508 Teuchos::OSTab tab2(out);
2509 out << "Scalar: " << TypeNameTraits<Scalar>::name() << endl
2510 << "LO: " << TypeNameTraits<LO>::name() << endl
2511 << "GO: " << TypeNameTraits<GO>::name() << endl
2512 << "Node: " << TypeNameTraits<Node>::name() << endl;
2513 }
2514 out << "Label: \"" << this->getObjectLabel() << "\"" << endl
2515 << "Global dimensions: ["
2516 << graph_.getDomainMap()->getGlobalNumElements() << ", "
2517 << graph_.getRangeMap()->getGlobalNumElements() << "]" << endl;
2518
2519 const LO blockSize = getBlockSize();
2520 out << "Block size: " << blockSize << endl;
2521
2522 // constituent objects
2523 if (vl >= VERB_MEDIUM) {
2524 const Teuchos::Comm<int>& comm = *(graph_.getMap()->getComm());
2525 const int myRank = comm.getRank();
2526 if (myRank == 0) {
2527 out << "Row Map:" << endl;
2528 }
2529 getRowMap()->describe(out, vl);
2530
2531 if (!getColMap().is_null()) {
2532 if (getColMap() == getRowMap()) {
2533 if (myRank == 0) {
2534 out << "Column Map: same as row Map" << endl;
2535 }
2536 } else {
2537 if (myRank == 0) {
2538 out << "Column Map:" << endl;
2539 }
2540 getColMap()->describe(out, vl);
2541 }
2542 }
2543 if (!getDomainMap().is_null()) {
2544 if (getDomainMap() == getRowMap()) {
2545 if (myRank == 0) {
2546 out << "Domain Map: same as row Map" << endl;
2547 }
2548 } else if (getDomainMap() == getColMap()) {
2549 if (myRank == 0) {
2550 out << "Domain Map: same as column Map" << endl;
2551 }
2552 } else {
2553 if (myRank == 0) {
2554 out << "Domain Map:" << endl;
2555 }
2556 getDomainMap()->describe(out, vl);
2557 }
2558 }
2559 if (!getRangeMap().is_null()) {
2560 if (getRangeMap() == getDomainMap()) {
2561 if (myRank == 0) {
2562 out << "Range Map: same as domain Map" << endl;
2563 }
2564 } else if (getRangeMap() == getRowMap()) {
2565 if (myRank == 0) {
2566 out << "Range Map: same as row Map" << endl;
2567 }
2568 } else {
2569 if (myRank == 0) {
2570 out << "Range Map: " << endl;
2571 }
2572 getRangeMap()->describe(out, vl);
2573 }
2574 }
2575 }
2576
2577 if (vl >= VERB_EXTREME) {
2578 const Teuchos::Comm<int>& comm = *(graph_.getMap()->getComm());
2579 const int myRank = comm.getRank();
2580 const int numProcs = comm.getSize();
2581
2582 // Print the calling process' data to the given output stream.
2583 RCP<std::ostringstream> lclOutStrPtr(new std::ostringstream());
2584 RCP<FancyOStream> osPtr = getFancyOStream(lclOutStrPtr);
2585 FancyOStream& os = *osPtr;
2586 os << "Process " << myRank << ":" << endl;
2587 Teuchos::OSTab tab2(os);
2588
2589 const map_type& meshRowMap = *(graph_.getRowMap());
2590 const map_type& meshColMap = *(graph_.getColMap());
2591 for (LO meshLclRow = meshRowMap.getMinLocalIndex();
2592 meshLclRow <= meshRowMap.getMaxLocalIndex();
2593 ++meshLclRow) {
2594 const GO meshGblRow = meshRowMap.getGlobalElement(meshLclRow);
2595 os << "Row " << meshGblRow << ": {";
2596
2597 local_inds_host_view_type lclColInds;
2598 values_host_view_type vals;
2599 LO numInds = 0;
2600 this->getLocalRowView(meshLclRow, lclColInds, vals);
2601 numInds = lclColInds.extent(0);
2602
2603 for (LO k = 0; k < numInds; ++k) {
2604 const GO gblCol = meshColMap.getGlobalElement(lclColInds[k]);
2605
2606 os << "Col " << gblCol << ": [";
2607 for (LO i = 0; i < blockSize; ++i) {
2608 for (LO j = 0; j < blockSize; ++j) {
2609 os << vals[blockSize * blockSize * k + i * blockSize + j];
2610 if (j + 1 < blockSize) {
2611 os << ", ";
2612 }
2613 }
2614 if (i + 1 < blockSize) {
2615 os << "; ";
2616 }
2617 }
2618 os << "]";
2619 if (k + 1 < numInds) {
2620 os << ", ";
2621 }
2622 }
2623 os << "}" << endl;
2624 }
2625
2626 // Print data on Process 0. This will automatically respect the
2627 // current indentation level.
2628 if (myRank == 0) {
2629 out << lclOutStrPtr->str();
2630 lclOutStrPtr = Teuchos::null; // clear it to save space
2631 }
2632
2633 const int sizeTag = 1337;
2634 const int dataTag = 1338;
2635
2636 ArrayRCP<char> recvDataBuf; // only used on Process 0
2637
2638 // Send string sizes and data from each process in turn to
2639 // Process 0, and print on that process.
2640 for (int p = 1; p < numProcs; ++p) {
2641 if (myRank == 0) {
2642 // Receive the incoming string's length.
2643 ArrayRCP<size_t> recvSize(1);
2644 recvSize[0] = 0;
2645 RCP<CommRequest<int> > recvSizeReq =
2646 ireceive<int, size_t>(recvSize, p, sizeTag, comm);
2647 wait<int>(comm, outArg(recvSizeReq));
2648 const size_t numCharsToRecv = recvSize[0];
2649
2650 // Allocate space for the string to receive. Reuse receive
2651 // buffer space if possible. We can do this because in the
2652 // current implementation, we only have one receive in
2653 // flight at a time. Leave space for the '\0' at the end,
2654 // in case the sender doesn't send it.
2655 if (static_cast<size_t>(recvDataBuf.size()) < numCharsToRecv + 1) {
2656 recvDataBuf.resize(numCharsToRecv + 1);
2657 }
2658 ArrayRCP<char> recvData = recvDataBuf.persistingView(0, numCharsToRecv);
2659 // Post the receive of the actual string data.
2660 RCP<CommRequest<int> > recvDataReq =
2661 ireceive<int, char>(recvData, p, dataTag, comm);
2662 wait<int>(comm, outArg(recvDataReq));
2663
2664 // Print the received data. This will respect the current
2665 // indentation level. Make sure that the string is
2666 // null-terminated.
2667 recvDataBuf[numCharsToRecv] = '\0';
2668 out << recvDataBuf.getRawPtr();
2669 } else if (myRank == p) { // if I am not Process 0, and my rank is p
2670 // This deep-copies the string at most twice, depending on
2671 // whether std::string reference counts internally (it
2672 // generally does, so this won't deep-copy at all).
2673 const std::string stringToSend = lclOutStrPtr->str();
2674 lclOutStrPtr = Teuchos::null; // clear original to save space
2675
2676 // Send the string's length to Process 0.
2677 const size_t numCharsToSend = stringToSend.size();
2678 ArrayRCP<size_t> sendSize(1);
2679 sendSize[0] = numCharsToSend;
2680 RCP<CommRequest<int> > sendSizeReq =
2681 isend<int, size_t>(sendSize, 0, sizeTag, comm);
2682 wait<int>(comm, outArg(sendSizeReq));
2683
2684 // Send the actual string to Process 0. We know that the
2685 // string has length > 0, so it's save to take the address
2686 // of the first entry. Make a nonowning ArrayRCP to hold
2687 // the string. Process 0 will add a null termination
2688 // character at the end of the string, after it receives the
2689 // message.
2690 ArrayRCP<const char> sendData(&stringToSend[0], 0, numCharsToSend, false);
2691 RCP<CommRequest<int> > sendDataReq =
2692 isend<int, char>(sendData, 0, dataTag, comm);
2693 wait<int>(comm, outArg(sendDataReq));
2694 }
2695 } // for each process rank p other than 0
2696 } // extreme verbosity level (print the whole matrix)
2697}
2698
2699template <class Scalar, class LO, class GO, class Node>
2700Teuchos::RCP<const Teuchos::Comm<int> >
2701BlockCrsMatrix<Scalar, LO, GO, Node>::
2702 getComm() const {
2703 return graph_.getComm();
2704}
2705
2706template <class Scalar, class LO, class GO, class Node>
2707global_size_t
2708BlockCrsMatrix<Scalar, LO, GO, Node>::
2709 getGlobalNumCols() const {
2710 return graph_.getGlobalNumCols();
2711}
2712
2713template <class Scalar, class LO, class GO, class Node>
2714size_t
2715BlockCrsMatrix<Scalar, LO, GO, Node>::
2716 getLocalNumCols() const {
2717 return graph_.getLocalNumCols();
2718}
2719
2720template <class Scalar, class LO, class GO, class Node>
2721GO BlockCrsMatrix<Scalar, LO, GO, Node>::
2722 getIndexBase() const {
2723 return graph_.getIndexBase();
2724}
2725
2726template <class Scalar, class LO, class GO, class Node>
2727global_size_t
2728BlockCrsMatrix<Scalar, LO, GO, Node>::
2729 getGlobalNumEntries() const {
2730 return graph_.getGlobalNumEntries();
2731}
2732
2733template <class Scalar, class LO, class GO, class Node>
2734size_t
2735BlockCrsMatrix<Scalar, LO, GO, Node>::
2736 getLocalNumEntries() const {
2737 return graph_.getLocalNumEntries();
2738}
2739
2740template <class Scalar, class LO, class GO, class Node>
2741size_t
2742BlockCrsMatrix<Scalar, LO, GO, Node>::
2743 getNumEntriesInGlobalRow(GO globalRow) const {
2744 return graph_.getNumEntriesInGlobalRow(globalRow);
2745}
2746
2747template <class Scalar, class LO, class GO, class Node>
2748size_t
2749BlockCrsMatrix<Scalar, LO, GO, Node>::
2750 getGlobalMaxNumRowEntries() const {
2751 return graph_.getGlobalMaxNumRowEntries();
2752}
2753
2754template <class Scalar, class LO, class GO, class Node>
2755bool BlockCrsMatrix<Scalar, LO, GO, Node>::
2756 hasColMap() const {
2757 return graph_.hasColMap();
2758}
2759
2760template <class Scalar, class LO, class GO, class Node>
2761bool BlockCrsMatrix<Scalar, LO, GO, Node>::
2762 isLocallyIndexed() const {
2763 return graph_.isLocallyIndexed();
2764}
2765
2766template <class Scalar, class LO, class GO, class Node>
2767bool BlockCrsMatrix<Scalar, LO, GO, Node>::
2768 isGloballyIndexed() const {
2769 return graph_.isGloballyIndexed();
2770}
2771
2772template <class Scalar, class LO, class GO, class Node>
2773bool BlockCrsMatrix<Scalar, LO, GO, Node>::
2774 isFillComplete() const {
2775 return graph_.isFillComplete();
2776}
2777
2778template <class Scalar, class LO, class GO, class Node>
2779bool BlockCrsMatrix<Scalar, LO, GO, Node>::
2780 supportsRowViews() const {
2781 return false;
2782}
2783
2784template <class Scalar, class LO, class GO, class Node>
2785void BlockCrsMatrix<Scalar, LO, GO, Node>::
2786 getGlobalRowCopy(GO /*GlobalRow*/,
2787 nonconst_global_inds_host_view_type& /*Indices*/,
2788 nonconst_values_host_view_type& /*Values*/,
2789 size_t& /*NumEntries*/) const {
2790 TEUCHOS_TEST_FOR_EXCEPTION(
2791 true, std::logic_error,
2792 "Tpetra::BlockCrsMatrix::getGlobalRowCopy: "
2793 "This class doesn't support global matrix indexing.");
2794}
2795
2796template <class Scalar, class LO, class GO, class Node>
2797void BlockCrsMatrix<Scalar, LO, GO, Node>::
2798 getGlobalRowView(GO /* GlobalRow */,
2799 global_inds_host_view_type& /* indices */,
2800 values_host_view_type& /* values */) const {
2801 TEUCHOS_TEST_FOR_EXCEPTION(
2802 true, std::logic_error,
2803 "Tpetra::BlockCrsMatrix::getGlobalRowView: "
2804 "This class doesn't support global matrix indexing.");
2805}
2806
2807template <class Scalar, class LO, class GO, class Node>
2808void BlockCrsMatrix<Scalar, LO, GO, Node>::
2809 getLocalRowView(LO localRowInd,
2810 local_inds_host_view_type& colInds,
2811 values_host_view_type& vals) const {
2812 if (!rowMeshMap_.isNodeLocalElement(localRowInd)) {
2813 colInds = local_inds_host_view_type();
2814 vals = values_host_view_type();
2815 } else {
2816 const size_t absBlockOffsetStart = ptrHost_[localRowInd];
2817 const size_t numInds = ptrHost_[localRowInd + 1] - absBlockOffsetStart;
2818 colInds = local_inds_host_view_type(indHost_.data() + absBlockOffsetStart, numInds);
2819
2820 vals = getValuesHost(localRowInd);
2821 }
2822}
2823
2824template <class Scalar, class LO, class GO, class Node>
2825void BlockCrsMatrix<Scalar, LO, GO, Node>::
2826 getLocalRowViewNonConst(LO localRowInd,
2827 local_inds_host_view_type& colInds,
2828 nonconst_values_host_view_type& vals) const {
2829 if (!rowMeshMap_.isNodeLocalElement(localRowInd)) {
2830 colInds = nonconst_local_inds_host_view_type();
2831 vals = nonconst_values_host_view_type();
2832 } else {
2833 const size_t absBlockOffsetStart = ptrHost_[localRowInd];
2834 const size_t numInds = ptrHost_[localRowInd + 1] - absBlockOffsetStart;
2835 colInds = local_inds_host_view_type(indHost_.data() + absBlockOffsetStart, numInds);
2836
2837 using this_BCRS_type = BlockCrsMatrix<Scalar, LO, GO, Node>;
2838 vals = const_cast<this_BCRS_type&>(*this).getValuesHostNonConst(localRowInd);
2839 }
2840}
2841
2842template <class Scalar, class LO, class GO, class Node>
2843void BlockCrsMatrix<Scalar, LO, GO, Node>::
2844 getLocalDiagCopy(::Tpetra::Vector<Scalar, LO, GO, Node>& diag) const {
2845 const size_t lclNumMeshRows = graph_.getLocalNumRows();
2846
2847 Kokkos::View<size_t*, device_type> diagOffsets("diagOffsets", lclNumMeshRows);
2848 graph_.getLocalDiagOffsets(diagOffsets);
2849
2850 // The code below works on host, so use a host View.
2851 auto diagOffsetsHost = Kokkos::create_mirror_view(diagOffsets);
2852 // DEEP_COPY REVIEW - DEVICE-TO-HOSTMIRROR
2853 Kokkos::deep_copy(execution_space(), diagOffsetsHost, diagOffsets);
2854
2855 auto vals_host_out = val_.getHostView(Access::ReadOnly);
2856 const impl_scalar_type* vals_host_out_raw = vals_host_out.data();
2857
2858 // TODO amk: This is a temporary measure to make the code run with Ifpack2
2859 size_t rowOffset = 0;
2860 size_t offset = 0;
2861 LO bs = getBlockSize();
2862 execution_space().fence(); // ensure deep_copy to diagOffsetsHost is done
2863 for (size_t r = 0; r < getLocalNumRows(); r++) {
2864 // move pointer to start of diagonal block
2865 offset = rowOffset + diagOffsetsHost(r) * bs * bs;
2866 for (int b = 0; b < bs; b++) {
2867 diag.replaceLocalValue(r * bs + b, vals_host_out_raw[offset + b * (bs + 1)]);
2868 }
2869 // move pointer to start of next block row
2870 rowOffset += getNumEntriesInLocalRow(r) * bs * bs;
2871 }
2872
2873 // diag.template sync<memory_space> (); // sync vec of diag entries back to dev
2874}
2875
2876template <class Scalar, class LO, class GO, class Node>
2877void BlockCrsMatrix<Scalar, LO, GO, Node>::
2878 leftScale(const ::Tpetra::Vector<Scalar, LO, GO, Node>& /* x */) {
2879 TEUCHOS_TEST_FOR_EXCEPTION(
2880 true, std::logic_error,
2881 "Tpetra::BlockCrsMatrix::leftScale: "
2882 "not implemented.");
2883}
2884
2885template <class Scalar, class LO, class GO, class Node>
2886void BlockCrsMatrix<Scalar, LO, GO, Node>::
2887 rightScale(const ::Tpetra::Vector<Scalar, LO, GO, Node>& /* x */) {
2888 TEUCHOS_TEST_FOR_EXCEPTION(
2889 true, std::logic_error,
2890 "Tpetra::BlockCrsMatrix::rightScale: "
2891 "not implemented.");
2892}
2893
2894template <class Scalar, class LO, class GO, class Node>
2895Teuchos::RCP<const ::Tpetra::RowGraph<LO, GO, Node> >
2896BlockCrsMatrix<Scalar, LO, GO, Node>::
2897 getGraph() const {
2898 return graphRCP_;
2899}
2900
2901template <class Scalar, class LO, class GO, class Node>
2902typename ::Tpetra::RowMatrix<Scalar, LO, GO, Node>::mag_type
2903BlockCrsMatrix<Scalar, LO, GO, Node>::
2904 getFrobeniusNorm() const {
2905 TEUCHOS_TEST_FOR_EXCEPTION(
2906 true, std::logic_error,
2907 "Tpetra::BlockCrsMatrix::getFrobeniusNorm: "
2908 "not implemented.");
2909}
2910
2911} // namespace Tpetra
2912
2913//
2914// Explicit instantiation macro
2915//
2916// Must be expanded from within the Tpetra namespace!
2917//
2918#define TPETRA_BLOCKCRSMATRIX_INSTANT(S, LO, GO, NODE) \
2919 template class BlockCrsMatrix<S, LO, GO, NODE>;
2920
2921#endif // TPETRA_BLOCKCRSMATRIX_DEF_HPP
Tpetra_BlockView.hpp
Linear algebra kernels for small dense matrices and vectors.
Tpetra_Details_Behavior.hpp
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra's behavior.
Tpetra_Details_PackTraits.hpp
Declaration and generic definition of traits class that tells Tpetra::CrsMatrix how to pack and unpac...
Tpetra_Details_Profiling.hpp
Declaration of Tpetra::Details::Profiling, a scope guard for Kokkos Profiling.
Tpetra::BlockCrsMatrix
Sparse matrix whose entries are small dense square blocks, all of the same dimensions.
Definition Tpetra_BlockCrsMatrix_decl.hpp:132
Tpetra::BlockCrsMatrix::sumIntoLocalValues
LO sumIntoLocalValues(const LO localRowInd, const LO colInds[], const Scalar vals[], const LO numColInds) const
Sum into values at the given (mesh, i.e., block) column indices, in the given (mesh,...
Definition Tpetra_BlockCrsMatrix_def.hpp:649
Tpetra::BlockCrsMatrix::getBlockSize
virtual LO getBlockSize() const override
The number of degrees of freedom per mesh point.
Definition Tpetra_BlockCrsMatrix_decl.hpp:405
Tpetra::BlockCrsMatrix::exportAndFillComplete
void exportAndFillComplete(Teuchos::RCP< BlockCrsMatrix< Scalar, LO, GO, Node > > &destMatrix, const Export< LO, GO, Node > &exporter) const
Import from this to the given destination matrix, and make the result fill complete.
Definition Tpetra_BlockCrsMatrix_def.hpp:554
Tpetra::BlockCrsMatrix::isLocallyIndexed
virtual bool isLocallyIndexed() const override
Whether matrix indices are locally indexed.
Definition Tpetra_BlockCrsMatrix_def.hpp:2762
Tpetra::BlockCrsMatrix::isFillComplete
virtual bool isFillComplete() const override
Whether fillComplete() has been called.
Definition Tpetra_BlockCrsMatrix_def.hpp:2774
Tpetra::BlockCrsMatrix::getGlobalMaxNumRowEntries
virtual size_t getGlobalMaxNumRowEntries() const override
The maximum number of entries in any row over all processes in the matrix's communicator.
Definition Tpetra_BlockCrsMatrix_def.hpp:2750
Tpetra::BlockCrsMatrix::applyBlock
void applyBlock(const BlockMultiVector< Scalar, LO, GO, Node > &X, BlockMultiVector< Scalar, LO, GO, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, const Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), const Scalar beta=Teuchos::ScalarTraits< Scalar >::zero())
Version of apply() that takes BlockMultiVector input and output.
Definition Tpetra_BlockCrsMatrix_def.hpp:503
Tpetra::BlockCrsMatrix::replaceLocalValuesByOffsets
LO replaceLocalValuesByOffsets(const LO localRowInd, const ptrdiff_t offsets[], const Scalar vals[], const LO numOffsets) const
Like replaceLocalValues, but avoids computing row offsets.
Definition Tpetra_BlockCrsMatrix_def.hpp:713
Tpetra::BlockCrsMatrix::getLocalRowCopy
virtual void getLocalRowCopy(LO LocalRow, nonconst_local_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const override
Not implemented.
Definition Tpetra_BlockCrsMatrix_def.hpp:660
Tpetra::BlockCrsMatrix::errorMessages
std::string errorMessages() const
The current stream of error messages.
Definition Tpetra_BlockCrsMatrix_decl.hpp:650
Tpetra::BlockCrsMatrix::absMaxLocalValues
LO absMaxLocalValues(const LO localRowInd, const LO colInds[], const Scalar vals[], const LO numColInds) const
Variant of getLocalDiagCopy() that uses precomputed offsets and puts diagonal blocks in a 3-D Kokkos:...
Definition Tpetra_BlockCrsMatrix_def.hpp:637
Tpetra::BlockCrsMatrix::scalar_type
Scalar scalar_type
The type of entries in the matrix (that is, of each entry in each block).
Definition Tpetra_BlockCrsMatrix_decl.hpp:147
Tpetra::BlockCrsMatrix::getGlobalRowView
virtual void getGlobalRowView(GO GlobalRow, global_inds_host_view_type &indices, values_host_view_type &values) const override
Get a constant, nonpersisting, globally indexed view of the given row of the matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:2798
Tpetra::BlockCrsMatrix::describe
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Print a description of this object to the given output stream.
Definition Tpetra_BlockCrsMatrix_def.hpp:2472
Tpetra::BlockCrsMatrix::getLocalMaxNumRowEntries
size_t getLocalMaxNumRowEntries() const override
Maximum number of entries in any row of the matrix, on this process.
Definition Tpetra_BlockCrsMatrix_def.hpp:422
Tpetra::BlockCrsMatrix::local_ordinal_type
LO local_ordinal_type
The type of local indices.
Definition Tpetra_BlockCrsMatrix_decl.hpp:158
Tpetra::BlockCrsMatrix::getGlobalRowCopy
virtual void getGlobalRowCopy(GO GlobalRow, nonconst_global_inds_host_view_type &Indices, nonconst_values_host_view_type &Values, size_t &NumEntries) const override
Get a copy of the given global row's entries.
Definition Tpetra_BlockCrsMatrix_def.hpp:2786
Tpetra::BlockCrsMatrix::getLocalRowOffsets
LO getLocalRowOffsets(const LO localRowInd, ptrdiff_t offsets[], const LO colInds[], const LO numColInds) const
Get relative offsets corresponding to the given rows, given by local row index.
Definition Tpetra_BlockCrsMatrix_def.hpp:683
Tpetra::BlockCrsMatrix::getLocalNumCols
virtual size_t getLocalNumCols() const override
The number of columns needed to apply the forward operator on this node.
Definition Tpetra_BlockCrsMatrix_def.hpp:2716
Tpetra::BlockCrsMatrix::copyAndPermute
virtual void copyAndPermute(const SrcDistObject &sourceObj, const size_t numSameIDs, const Kokkos::DualView< const local_ordinal_type *, buffer_device_type > &permuteToLIDs, const Kokkos::DualView< const local_ordinal_type *, buffer_device_type > &permuteFromLIDs, const CombineMode CM) override
Definition Tpetra_BlockCrsMatrix_def.hpp:1244
Tpetra::BlockCrsMatrix::mv_type
::Tpetra::MultiVector< Scalar, LO, GO, node_type > mv_type
The implementation of MultiVector that this class uses.
Definition Tpetra_BlockCrsMatrix_decl.hpp:177
Tpetra::BlockCrsMatrix::getRangeMap
Teuchos::RCP< const map_type > getRangeMap() const override
Get the (point) range Map of this matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:386
Tpetra::BlockCrsMatrix::sumIntoLocalValuesByOffsets
LO sumIntoLocalValuesByOffsets(const LO localRowInd, const ptrdiff_t offsets[], const Scalar vals[], const LO numOffsets) const
Like sumIntoLocalValues, but avoids computing row offsets.
Definition Tpetra_BlockCrsMatrix_def.hpp:782
Tpetra::BlockCrsMatrix::hasColMap
virtual bool hasColMap() const override
Whether this matrix has a well-defined column Map.
Definition Tpetra_BlockCrsMatrix_def.hpp:2756
Tpetra::BlockCrsMatrix::execution_space
device_type::execution_space execution_space
The Kokkos execution space that this class uses.
Definition Tpetra_BlockCrsMatrix_decl.hpp:170
Tpetra::BlockCrsMatrix::getLocalNumRows
size_t getLocalNumRows() const override
get the local number of block rows
Definition Tpetra_BlockCrsMatrix_def.hpp:415
Tpetra::BlockCrsMatrix::getLocalNumEntries
virtual size_t getLocalNumEntries() const override
The local number of stored (structurally nonzero) entries.
Definition Tpetra_BlockCrsMatrix_def.hpp:2736
Tpetra::BlockCrsMatrix::getValuesHostNonConst
impl_scalar_type_dualview::t_host getValuesHostNonConst() const
Get the host or device View of the matrix's values (val_).
Definition Tpetra_BlockCrsMatrix_def.hpp:836
Tpetra::BlockCrsMatrix::little_block_type
Kokkos::View< impl_scalar_type **, Impl::BlockCrsMatrixLittleBlockArrayLayout, device_type, Kokkos::MemoryTraits< Kokkos::Unmanaged > > little_block_type
The type used to access nonconst matrix blocks.
Definition Tpetra_BlockCrsMatrix_decl.hpp:186
Tpetra::BlockCrsMatrix::getColMap
Teuchos::RCP< const map_type > getColMap() const override
get the (mesh) map for the columns of this block matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:401
Tpetra::BlockCrsMatrix::getFrobeniusNorm
virtual typename::Tpetra::RowMatrix< Scalar, LO, GO, Node >::mag_type getFrobeniusNorm() const override
The Frobenius norm of the matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:2904
Tpetra::BlockCrsMatrix::unpackAndCombine
virtual void unpackAndCombine(const Kokkos::DualView< const local_ordinal_type *, buffer_device_type > &importLIDs, Kokkos::DualView< packet_type *, buffer_device_type > imports, Kokkos::DualView< size_t *, buffer_device_type > numPacketsPerLID, const size_t constantNumPackets, const CombineMode combineMode) override
Definition Tpetra_BlockCrsMatrix_def.hpp:2113
Tpetra::BlockCrsMatrix::getLocalMatrixDevice
local_matrix_device_type getLocalMatrixDevice() const
Definition Tpetra_BlockCrsMatrix_def.hpp:905
Tpetra::BlockCrsMatrix::apply
void apply(const mv_type &X, mv_type &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const override
For this matrix A, compute Y := beta * Y + alpha * Op(A) * X.
Definition Tpetra_BlockCrsMatrix_def.hpp:428
Tpetra::BlockCrsMatrix::getGlobalNumCols
virtual global_size_t getGlobalNumCols() const override
The global number of columns of this matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:2709
Tpetra::BlockCrsMatrix::getLocalDiagOffsets
void getLocalDiagOffsets(const Kokkos::View< size_t *, device_type, Kokkos::MemoryUnmanaged > &offsets) const
Get offsets of the diagonal entries in the matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:597
Tpetra::BlockCrsMatrix::getRowMap
Teuchos::RCP< const map_type > getRowMap() const override
get the (mesh) map for the rows of this block matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:394
Tpetra::BlockCrsMatrix::description
std::string description() const override
One-line description of this object.
Definition Tpetra_BlockCrsMatrix_def.hpp:2451
Tpetra::BlockCrsMatrix::importAndFillComplete
void importAndFillComplete(Teuchos::RCP< BlockCrsMatrix< Scalar, LO, GO, Node > > &destMatrix, const Import< LO, GO, Node > &importer) const
Import from this to the given destination matrix, and make the result fill complete.
Definition Tpetra_BlockCrsMatrix_def.hpp:530
Tpetra::BlockCrsMatrix::getLocalRowView
void getLocalRowView(LO LocalRow, local_inds_host_view_type &indices, values_host_view_type &values) const override
Get a view of the (mesh, i.e., block) row, using local (mesh, i.e., block) indices.
Definition Tpetra_BlockCrsMatrix_def.hpp:2809
Tpetra::BlockCrsMatrix::getLocalDiagCopy
void getLocalDiagCopy(const Kokkos::View< impl_scalar_type ***, device_type, Kokkos::MemoryUnmanaged > &diag, const Kokkos::View< const size_t *, device_type, Kokkos::MemoryUnmanaged > &offsets) const
Variant of getLocalDiagCopy() that uses precomputed offsets and puts diagonal blocks in a 3-D Kokkos:...
Definition Tpetra_BlockCrsMatrix_def.hpp:604
Tpetra::BlockCrsMatrix::getNumEntriesInLocalRow
size_t getNumEntriesInLocalRow(const LO localRowInd) const override
Return the number of entries in the given row on the calling process.
Definition Tpetra_BlockCrsMatrix_def.hpp:893
Tpetra::BlockCrsMatrix::supportsRowViews
virtual bool supportsRowViews() const override
Whether this object implements getLocalRowView() and getGlobalRowView().
Definition Tpetra_BlockCrsMatrix_def.hpp:2780
Tpetra::BlockCrsMatrix::getGraph
virtual Teuchos::RCP< const ::Tpetra::RowGraph< LO, GO, Node > > getGraph() const override
Get the (mesh) graph.
Definition Tpetra_BlockCrsMatrix_def.hpp:2897
Tpetra::BlockCrsMatrix::replaceLocalValues
LO replaceLocalValues(const LO localRowInd, const LO colInds[], const Scalar vals[], const LO numColInds) const
Replace values at the given (mesh, i.e., block) column indices, in the given (mesh,...
Definition Tpetra_BlockCrsMatrix_def.hpp:585
Tpetra::BlockCrsMatrix::getNumEntriesInGlobalRow
virtual size_t getNumEntriesInGlobalRow(GO globalRow) const override
The current number of entries on the calling process in the specified global row.
Definition Tpetra_BlockCrsMatrix_def.hpp:2743
Tpetra::BlockCrsMatrix::buffer_device_type
typename DistObject< Scalar, LO, GO, Node >::buffer_device_type buffer_device_type
Kokkos::Device specialization for communication buffers.
Definition Tpetra_BlockCrsMatrix_decl.hpp:746
Tpetra::BlockCrsMatrix::device_type
Node::device_type device_type
The Kokkos::Device specialization that this class uses.
Definition Tpetra_BlockCrsMatrix_decl.hpp:168
Tpetra::BlockCrsMatrix::setAllToScalar
void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to alpha.
Definition Tpetra_BlockCrsMatrix_def.hpp:578
Tpetra::BlockCrsMatrix::const_little_block_type
Kokkos::View< const impl_scalar_type **, Impl::BlockCrsMatrixLittleBlockArrayLayout, device_type, Kokkos::MemoryTraits< Kokkos::Unmanaged > > const_little_block_type
The type used to access const matrix blocks.
Definition Tpetra_BlockCrsMatrix_decl.hpp:194
Tpetra::BlockCrsMatrix::getLocalRowViewNonConst
void getLocalRowViewNonConst(LO LocalRow, local_inds_host_view_type &indices, nonconst_values_host_view_type &values) const
Definition Tpetra_BlockCrsMatrix_def.hpp:2826
Tpetra::BlockCrsMatrix::packet_type
char packet_type
Implementation detail; tells.
Definition Tpetra_BlockCrsMatrix_decl.hpp:140
Tpetra::BlockCrsMatrix::packAndPrepare
virtual void packAndPrepare(const SrcDistObject &sourceObj, const Kokkos::DualView< const local_ordinal_type *, buffer_device_type > &exportLIDs, Kokkos::DualView< packet_type *, buffer_device_type > &exports, Kokkos::DualView< size_t *, buffer_device_type > numPacketsPerLID, size_t &constantNumPackets) override
Definition Tpetra_BlockCrsMatrix_def.hpp:1861
Tpetra::BlockCrsMatrix::map_type
::Tpetra::Map< LO, GO, node_type > map_type
The implementation of Map that this class uses.
Definition Tpetra_BlockCrsMatrix_decl.hpp:175
Tpetra::BlockCrsMatrix::rightScale
virtual void rightScale(const ::Tpetra::Vector< Scalar, LO, GO, Node > &x) override
Scale the RowMatrix on the right with the given Vector x.
Definition Tpetra_BlockCrsMatrix_def.hpp:2887
Tpetra::BlockCrsMatrix::getIndexBase
virtual GO getIndexBase() const override
The index base for global indices in this matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:2722
Tpetra::BlockCrsMatrix::getGlobalNumEntries
virtual global_size_t getGlobalNumEntries() const override
The global number of stored (structurally nonzero) entries.
Definition Tpetra_BlockCrsMatrix_def.hpp:2729
Tpetra::BlockCrsMatrix::impl_scalar_type
typename BMV::impl_scalar_type impl_scalar_type
The implementation type of entries in the matrix.
Definition Tpetra_BlockCrsMatrix_decl.hpp:155
Tpetra::BlockCrsMatrix::getComm
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const override
The communicator over which this matrix is distributed.
Definition Tpetra_BlockCrsMatrix_def.hpp:2702
Tpetra::BlockCrsMatrix::crs_graph_type
::Tpetra::CrsGraph< LO, GO, node_type > crs_graph_type
The implementation of CrsGraph that this class uses.
Definition Tpetra_BlockCrsMatrix_decl.hpp:179
Tpetra::BlockCrsMatrix::getGlobalNumRows
global_size_t getGlobalNumRows() const override
get the global number of block rows
Definition Tpetra_BlockCrsMatrix_def.hpp:408
Tpetra::BlockCrsMatrix::leftScale
virtual void leftScale(const ::Tpetra::Vector< Scalar, LO, GO, Node > &x) override
Scale the RowMatrix on the left with the given Vector x.
Definition Tpetra_BlockCrsMatrix_def.hpp:2878
Tpetra::BlockCrsMatrix::getDomainMap
Teuchos::RCP< const map_type > getDomainMap() const override
Get the (point) domain Map of this matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:378
Tpetra::BlockCrsMatrix::isGloballyIndexed
virtual bool isGloballyIndexed() const override
Whether matrix indices are globally indexed.
Definition Tpetra_BlockCrsMatrix_def.hpp:2768
Tpetra::BlockCrsMatrix::BlockCrsMatrix
BlockCrsMatrix()
Default constructor: Makes an empty block matrix.
Definition Tpetra_BlockCrsMatrix_def.hpp:150
Tpetra::BlockMultiVector
MultiVector for multiple degrees of freedom per mesh point.
Definition Tpetra_BlockMultiVector_decl.hpp:109
Tpetra::BlockMultiVector::getMultiVectorView
const mv_type & getMultiVectorView() const
Get a Tpetra::MultiVector that views this BlockMultiVector's data.
Definition Tpetra_BlockMultiVector_def.hpp:23
Tpetra::BlockMultiVector::makePointMap
static map_type makePointMap(const map_type &meshMap, const LO blockSize)
Create and return the point Map corresponding to the given mesh Map and block size.
Definition Tpetra_BlockMultiVector_def.hpp:175
Tpetra::BlockMultiVector::map_type
Tpetra::Map< LO, GO, Node > map_type
The specialization of Tpetra::Map that this class uses.
Definition Tpetra_BlockMultiVector_decl.hpp:120
Tpetra::BlockMultiVector::getBlockSize
LO getBlockSize() const
Get the number of degrees of freedom per mesh point.
Definition Tpetra_BlockMultiVector_decl.hpp:312
Tpetra::CrsGraph::getNumEntriesInLocalRow
size_t getNumEntriesInLocalRow(local_ordinal_type localRow) const override
Get the number of entries in the given row (local index).
Definition Tpetra_CrsGraph_def.hpp:2091
Tpetra::CrsGraph::getColMap
Teuchos::RCP< const map_type > getColMap() const override
Returns the Map that describes the column distribution in this graph.
Definition Tpetra_CrsGraph_def.hpp:898
Tpetra::CrsGraph< LO, GO, node_type >::import_type
::Tpetra::Import< LO, GO, node_type > import_type
Definition Tpetra_CrsGraph_decl.hpp:234
Tpetra::CrsGraph< LO, GO, Node >::local_graph_device_type
KokkosSparse::StaticCrsGraph< local_ordinal_type, Kokkos::LayoutLeft, device_type, void, size_t > local_graph_device_type
Definition Tpetra_CrsGraph_decl.hpp:224
Tpetra::CrsGraph::getDomainMap
Teuchos::RCP< const map_type > getDomainMap() const override
Returns the Map associated with the domain of this graph.
Definition Tpetra_CrsGraph_def.hpp:905
Tpetra::CrsGraph::getRangeMap
Teuchos::RCP< const map_type > getRangeMap() const override
Returns the Map associated with the domain of this graph.
Definition Tpetra_CrsGraph_def.hpp:912
Tpetra::CrsMatrix::getColMap
Teuchos::RCP< const map_type > getColMap() const override
The Map that describes the column distribution in this matrix.
Definition Tpetra_CrsMatrix_def.hpp:947
Tpetra::Details::Behavior
Description of Tpetra's behavior.
Definition Tpetra_Details_Behavior.hpp:92
Tpetra::DistObject< char, LO, GO, Node >::getMap
virtual Teuchos::RCP< const map_type > getMap() const
Definition Tpetra_DistObject_decl.hpp:550
Tpetra::Export
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Definition Tpetra_Export_decl.hpp:86
Tpetra::Import
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Definition Tpetra_Import_decl.hpp:77
Tpetra::Map::getGlobalElement
global_ordinal_type getGlobalElement(local_ordinal_type localIndex) const
The global index corresponding to the given local index.
Definition Tpetra_Map_def.hpp:1212
Tpetra::Map::isNodeLocalElement
bool isNodeLocalElement(local_ordinal_type localIndex) const
Whether the given local index is valid for this Map on the calling process.
Definition Tpetra_Map_def.hpp:1261
Tpetra::Map::locallySameAs
bool locallySameAs(const Map< local_ordinal_type, global_ordinal_type, node_type > &map) const
Is this Map locally the same as the input Map?
Definition Tpetra_Map_def.hpp:1358
Tpetra::Map::getLocalElement
local_ordinal_type getLocalElement(global_ordinal_type globalIndex) const
The local index corresponding to the given global index.
Definition Tpetra_Map_def.hpp:1190
Tpetra::Map::getMinLocalIndex
local_ordinal_type getMinLocalIndex() const
The minimum local index.
Definition Tpetra_Map_decl.hpp:586
Tpetra::Map::getMaxLocalIndex
local_ordinal_type getMaxLocalIndex() const
The maximum local index on the calling process.
Definition Tpetra_Map_decl.hpp:600
Tpetra::MultiVector::getLocalViewDevice
dual_view_type::t_dev::const_type getLocalViewDevice(Access::ReadOnlyStruct) const
Return a read-only, up-to-date view of this MultiVector's local data on device. This requires that th...
Definition Tpetra_MultiVector_def.hpp:3734
Tpetra::MultiVector::isConstantStride
bool isConstantStride() const
Whether this multivector has constant stride between columns.
Definition Tpetra_MultiVector_def.hpp:960
Tpetra::Vector
A distributed dense vector.
Definition Tpetra_Vector_decl.hpp:44
Tpetra::Vector::replaceLocalValue
void replaceLocalValue(const LocalOrdinal myRow, const Scalar &value)
Replace current value at the specified location with specified values.
Definition Tpetra_Vector_def.hpp:95
Tpetra::Details::disableWDVTracking
void disableWDVTracking()
Disable WrappedDualView reference-count tracking and syncing. Call this before entering a host-parall...
Definition Tpetra_Details_WrappedDualView.cpp:23
Tpetra::Details::enableWDVTracking
void enableWDVTracking()
Enable WrappedDualView reference-count tracking and syncing. Call this after exiting a host-parallel ...
Definition Tpetra_Details_WrappedDualView.cpp:17
Tpetra::Details::dualViewStatusToString
std::string dualViewStatusToString(const DualViewType &dv, const char name[])
Return the status of the given Kokkos::DualView, as a human-readable string.
Definition Tpetra_Util.hpp:1007
Tpetra::Impl
Namespace for new Tpetra features that are not ready for public release, but are ready for evaluation...
Definition Tpetra_Assembly_Helpers.hpp:15
Tpetra::Impl::BlockCrsMatrixLittleBlockArrayLayout
Kokkos::LayoutRight BlockCrsMatrixLittleBlockArrayLayout
give an option to use layoutleft
Definition Tpetra_BlockCrsMatrix_decl.hpp:123
Tpetra::Impl::absMax
KOKKOS_INLINE_FUNCTION void absMax(const ViewType2 &Y, const ViewType1 &X)
Implementation of Tpetra's ABSMAX CombineMode for the small dense blocks in BlockCrsMatrix,...
Definition Tpetra_BlockView.hpp:128
Tpetra
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Tpetra::exportAndFillCompleteCrsGraph
Teuchos::RCP< CrsGraphType > exportAndFillCompleteCrsGraph(const Teuchos::RCP< const CrsGraphType > &sourceGraph, const Export< typename CrsGraphType::local_ordinal_type, typename CrsGraphType::global_ordinal_type, typename CrsGraphType::node_type > &exporter, const Teuchos::RCP< const Map< typename CrsGraphType::local_ordinal_type, typename CrsGraphType::global_ordinal_type, typename CrsGraphType::node_type > > &domainMap=Teuchos::null, const Teuchos::RCP< const Map< typename CrsGraphType::local_ordinal_type, typename CrsGraphType::global_ordinal_type, typename CrsGraphType::node_type > > &rangeMap=Teuchos::null, const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Nonmember CrsGraph constructor that fuses Export and fillComplete().
Definition Tpetra_CrsGraph_decl.hpp:2766
Tpetra::AXPY
KOKKOS_INLINE_FUNCTION void AXPY(const CoefficientType &alpha, const ViewType1 &x, const ViewType2 &y)
y := y + alpha * x (dense vector or matrix update)
Definition Tpetra_BlockView.hpp:582
Tpetra::importAndFillCompleteCrsGraph
Teuchos::RCP< CrsGraphType > importAndFillCompleteCrsGraph(const Teuchos::RCP< const CrsGraphType > &sourceGraph, const Import< typename CrsGraphType::local_ordinal_type, typename CrsGraphType::global_ordinal_type, typename CrsGraphType::node_type > &importer, const Teuchos::RCP< const Map< typename CrsGraphType::local_ordinal_type, typename CrsGraphType::global_ordinal_type, typename CrsGraphType::node_type > > &domainMap=Teuchos::null, const Teuchos::RCP< const Map< typename CrsGraphType::local_ordinal_type, typename CrsGraphType::global_ordinal_type, typename CrsGraphType::node_type > > &rangeMap=Teuchos::null, const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Nonmember CrsGraph constructor that fuses Import and fillComplete().
Definition Tpetra_CrsGraph_decl.hpp:2644
Tpetra::global_size_t
size_t global_size_t
Global size_t object.
Definition Tpetra_ConfigDefs.hpp:77
Tpetra::combineModeToString
std::string combineModeToString(const CombineMode combineMode)
Human-readable string representation of the given CombineMode.
Definition Tpetra_CombineMode.cpp:61
Tpetra::COPY
KOKKOS_INLINE_FUNCTION void COPY(const ViewType1 &x, const ViewType2 &y)
Deep copy x into y, where x and y are either rank 1 (vectors) or rank 2 (matrices) with the same dime...
Definition Tpetra_BlockView.hpp:610
Tpetra::CombineMode
CombineMode
Rule for combining data in an Import or Export.
Definition Tpetra_CombineMode.hpp:65
Tpetra::REPLACE
@ REPLACE
Replace existing values with new values.
Definition Tpetra_CombineMode.hpp:68
Tpetra::ADD
@ ADD
Sum new values.
Definition Tpetra_CombineMode.hpp:66
Tpetra::ABSMAX
@ ABSMAX
Replace old value with maximum of magnitudes of old and new values.
Definition Tpetra_CombineMode.hpp:69
Tpetra::INSERT
@ INSERT
Insert new values that don't currently exist.
Definition Tpetra_CombineMode.hpp:67
Tpetra::ZERO
@ ZERO
Replace old values with zero.
Definition Tpetra_CombineMode.hpp:70
Tpetra::Details::PackTraits
Traits class for packing / unpacking data of type T.
Definition Tpetra_Details_PackTraits.hpp:29
Generated by doxygen 1.15.0