Tpetra parallel linear algebra  Version of the Day
Tpetra_Experimental_BlockMultiVector_def.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 
42 #ifndef TPETRA_EXPERIMENTAL_BLOCKMULTIVECTOR_DEF_HPP
43 #define TPETRA_EXPERIMENTAL_BLOCKMULTIVECTOR_DEF_HPP
44 
45 #include "Tpetra_Experimental_BlockMultiVector_decl.hpp"
46 
47 namespace { // anonymous
48 
60  template<class MultiVectorType>
61  struct RawPtrFromMultiVector {
62  typedef typename MultiVectorType::impl_scalar_type impl_scalar_type;
63 
64  static impl_scalar_type* getRawPtr (MultiVectorType& X) {
65  typedef typename MultiVectorType::dual_view_type dual_view_type;
66  typedef typename dual_view_type::t_host::memory_space host_memory_space;
67 
68  // We're getting a nonconst View, so mark the MultiVector as
69  // modified on the host. This will throw an exception if the
70  // MultiVector is already modified on the host.
71  X.template modify<host_memory_space> ();
72 
73  dual_view_type X_view = X.getDualView ();
74  impl_scalar_type* X_raw = X_view.h_view.ptr_on_device ();
75  return X_raw;
76  }
77  };
78 
91  template<class S, class LO, class GO, class N>
93  getRawPtrFromMultiVector (Tpetra::MultiVector<S, LO, GO, N>& X) {
95  return RawPtrFromMultiVector<MV>::getRawPtr (X);
96  }
97 
98 } // namespace (anonymous)
99 
100 namespace Tpetra {
101 namespace Experimental {
102 
103 template<class Scalar, class LO, class GO, class Node>
104 typename BlockMultiVector<Scalar, LO, GO, Node>::mv_type
105 BlockMultiVector<Scalar, LO, GO, Node>::
106 getMultiVectorView ()
107 {
108  // Make sure that mv_ has view semantics.
109  mv_.setCopyOrView (Teuchos::View);
110  // Now the one-argument copy constructor will make a shallow copy,
111  // and those view semantics will persist in all of its offspring.
112  return mv_;
113 }
114 
115 template<class Scalar, class LO, class GO, class Node>
116 Teuchos::RCP<const BlockMultiVector<Scalar, LO, GO, Node> >
119 {
121  const BMV* src_bmv = dynamic_cast<const BMV*> (&src);
122  TEUCHOS_TEST_FOR_EXCEPTION(
123  src_bmv == NULL, std::invalid_argument, "Tpetra::Experimental::"
124  "BlockMultiVector: The source object of an Import or Export to a "
125  "BlockMultiVector, must also be a BlockMultiVector.");
126  return Teuchos::rcp (src_bmv, false); // nonowning RCP
127 }
128 
129 template<class Scalar, class LO, class GO, class Node>
131 BlockMultiVector (const map_type& meshMap,
132  const LO blockSize,
133  const LO numVecs) :
134  dist_object_type (Teuchos::rcp (new map_type (meshMap))), // shallow copy
135  meshMap_ (meshMap),
136  pointMap_ (makePointMap (meshMap, blockSize)),
137  mv_ (Teuchos::rcpFromRef (pointMap_), numVecs), // nonowning RCP is OK, since pointMap_ won't go away
138  mvData_ (getRawPtrFromMultiVector (mv_)),
139  blockSize_ (blockSize)
140 {
141  // Make sure that mv_ has view semantics.
142  mv_.setCopyOrView (Teuchos::View);
143 }
144 
145 template<class Scalar, class LO, class GO, class Node>
147 BlockMultiVector (const map_type& meshMap,
148  const map_type& pointMap,
149  const LO blockSize,
150  const LO numVecs) :
151  dist_object_type (Teuchos::rcp (new map_type (meshMap))), // shallow copy
152  meshMap_ (meshMap),
153  pointMap_ (pointMap),
154  mv_ (Teuchos::rcpFromRef (pointMap_), numVecs),
155  mvData_ (getRawPtrFromMultiVector (mv_)),
156  blockSize_ (blockSize)
157 {
158  // Make sure that mv_ has view semantics.
159  mv_.setCopyOrView (Teuchos::View);
160 }
161 
162 template<class Scalar, class LO, class GO, class Node>
165  const map_type& meshMap,
166  const LO blockSize) :
167  dist_object_type (Teuchos::rcp (new map_type (meshMap))), // shallow copy
168  meshMap_ (meshMap),
169  mvData_ (NULL), // just for now
170  blockSize_ (blockSize)
171 {
172  using Teuchos::RCP;
173 
174  if (X_mv.getCopyOrView () == Teuchos::View) {
175  // The input MultiVector has view semantics, so assignment just
176  // does a shallow copy.
177  mv_ = X_mv;
178  }
179  else if (X_mv.getCopyOrView () == Teuchos::Copy) {
180  // The input MultiVector has copy semantics. We can't change
181  // that, but we can make a view of the input MultiVector and
182  // change the view to have view semantics. (That sounds silly;
183  // shouldn't views always have view semantics? but remember that
184  // "view semantics" just governs the default behavior of the copy
185  // constructor and assignment operator.)
186  RCP<const mv_type> X_view_const;
187  const size_t numCols = X_mv.getNumVectors ();
188  if (numCols == 0) {
189  Teuchos::Array<size_t> cols (0); // view will have zero columns
190  X_view_const = X_mv.subView (cols ());
191  } else { // Range1D is an inclusive range
192  X_view_const = X_mv.subView (Teuchos::Range1D (0, numCols-1));
193  }
194  TEUCHOS_TEST_FOR_EXCEPTION(
195  X_view_const.is_null (), std::logic_error, "Tpetra::Experimental::"
196  "BlockMultiVector constructor: X_mv.subView(...) returned null. This "
197  "should never happen. Please report this bug to the Tpetra developers.");
198 
199  // It's perfectly OK to cast away const here. Those view methods
200  // should be marked const anyway, because views can't change the
201  // allocation (just the entries themselves).
202  RCP<mv_type> X_view = Teuchos::rcp_const_cast<mv_type> (X_view_const);
203  X_view->setCopyOrView (Teuchos::View);
204  TEUCHOS_TEST_FOR_EXCEPTION(
205  X_view->getCopyOrView () != Teuchos::View, std::logic_error, "Tpetra::"
206  "Experimental::BlockMultiVector constructor: We just set a MultiVector "
207  "to have view semantics, but it claims that it doesn't have view "
208  "semantics. This should never happen. "
209  "Please report this bug to the Tpetra developers.");
210  mv_ = *X_view; // MultiVector::operator= does a shallow copy here
211  }
212 
213  // At this point, mv_ has been assigned, so we can ignore X_mv.
214  Teuchos::RCP<const map_type> pointMap = mv_.getMap ();
215  if (! pointMap.is_null ()) {
216  pointMap_ = *pointMap; // Map::operator= also does a shallow copy
217  }
218  mvData_ = getRawPtrFromMultiVector (mv_);
219 }
220 
221 template<class Scalar, class LO, class GO, class Node>
224  const map_type& newMeshMap,
225  const map_type& newPointMap,
226  const size_t offset) :
227  dist_object_type (Teuchos::rcp (new map_type (newMeshMap))), // shallow copy
228  meshMap_ (newMeshMap),
229  pointMap_ (newPointMap),
230  mv_ (X.mv_, newPointMap, offset * X.getBlockSize ()), // MV "offset view" constructor
231  mvData_ (getRawPtrFromMultiVector (mv_)),
232  blockSize_ (X.getBlockSize ())
233 {
234  // Make sure that mv_ has view semantics.
235  mv_.setCopyOrView (Teuchos::View);
236 }
237 
238 template<class Scalar, class LO, class GO, class Node>
241  const map_type& newMeshMap,
242  const size_t offset) :
243  dist_object_type (Teuchos::rcp (new map_type (newMeshMap))), // shallow copy
244  meshMap_ (newMeshMap),
245  pointMap_ (makePointMap (newMeshMap, X.getBlockSize ())),
246  mv_ (X.mv_, pointMap_, offset * X.getBlockSize ()), // MV "offset view" constructor
247  mvData_ (getRawPtrFromMultiVector (mv_)),
248  blockSize_ (X.getBlockSize ())
249 {
250  // Make sure that mv_ has view semantics.
251  mv_.setCopyOrView (Teuchos::View);
252 }
253 
254 template<class Scalar, class LO, class GO, class Node>
257  dist_object_type (Teuchos::null),
258  mvData_ (NULL),
259  blockSize_ (0)
260 {
261  // Make sure that mv_ has view semantics.
262  mv_.setCopyOrView (Teuchos::View);
263 }
264 
265 template<class Scalar, class LO, class GO, class Node>
268 makePointMap (const map_type& meshMap, const LO blockSize)
269 {
270  typedef Tpetra::global_size_t GST;
271  typedef typename Teuchos::ArrayView<const GO>::size_type size_type;
272 
273  const GST gblNumMeshMapInds =
274  static_cast<GST> (meshMap.getGlobalNumElements ());
275  const size_t lclNumMeshMapIndices =
276  static_cast<size_t> (meshMap.getNodeNumElements ());
277  const GST gblNumPointMapInds =
278  gblNumMeshMapInds * static_cast<GST> (blockSize);
279  const size_t lclNumPointMapInds =
280  lclNumMeshMapIndices * static_cast<size_t> (blockSize);
281  const GO indexBase = meshMap.getIndexBase ();
282 
283  if (meshMap.isContiguous ()) {
284  return map_type (gblNumPointMapInds, lclNumPointMapInds, indexBase,
285  meshMap.getComm (), meshMap.getNode ());
286  }
287  else {
288  // "Hilbert's Hotel" trick: multiply each process' GIDs by
289  // blockSize, and fill in. That ensures correctness even if the
290  // mesh Map is overlapping.
291  Teuchos::ArrayView<const GO> lclMeshGblInds = meshMap.getNodeElementList ();
292  const size_type lclNumMeshGblInds = lclMeshGblInds.size ();
293  Teuchos::Array<GO> lclPointGblInds (lclNumPointMapInds);
294  for (size_type g = 0; g < lclNumMeshGblInds; ++g) {
295  const GO meshGid = lclMeshGblInds[g];
296  const GO pointGidStart = indexBase +
297  (meshGid - indexBase) * static_cast<GO> (blockSize);
298  const size_type offset = g * static_cast<size_type> (blockSize);
299  for (LO k = 0; k < blockSize; ++k) {
300  const GO pointGid = pointGidStart + static_cast<GO> (k);
301  lclPointGblInds[offset + static_cast<size_type> (k)] = pointGid;
302  }
303  }
304  return map_type (gblNumPointMapInds, lclPointGblInds (), indexBase,
305  meshMap.getComm (), meshMap.getNode ());
306  }
307 }
308 
309 
310 template<class Scalar, class LO, class GO, class Node>
311 void
313 replaceLocalValuesImpl (const LO localRowIndex,
314  const LO colIndex,
315  const Scalar vals[]) const
316 {
317  little_vec_type X_dst = getLocalBlock (localRowIndex, colIndex);
318  const LO strideX = 1;
319  const_little_vec_type X_src (reinterpret_cast<const impl_scalar_type*> (vals),
320  getBlockSize (), strideX);
321  deep_copy (X_dst, X_src);
322 }
323 
324 
325 template<class Scalar, class LO, class GO, class Node>
326 bool
328 replaceLocalValues (const LO localRowIndex,
329  const LO colIndex,
330  const Scalar vals[]) const
331 {
332  if (! meshMap_.isNodeLocalElement (localRowIndex)) {
333  return false;
334  } else {
335  replaceLocalValuesImpl (localRowIndex, colIndex, vals);
336  return true;
337  }
338 }
339 
340 template<class Scalar, class LO, class GO, class Node>
341 bool
343 replaceGlobalValues (const GO globalRowIndex,
344  const LO colIndex,
345  const Scalar vals[]) const
346 {
347  const LO localRowIndex = meshMap_.getLocalElement (globalRowIndex);
348  if (localRowIndex == Teuchos::OrdinalTraits<LO>::invalid ()) {
349  return false;
350  } else {
351  replaceLocalValuesImpl (localRowIndex, colIndex, vals);
352  return true;
353  }
354 }
355 
356 template<class Scalar, class LO, class GO, class Node>
357 void
359 sumIntoLocalValuesImpl (const LO localRowIndex,
360  const LO colIndex,
361  const Scalar vals[]) const
362 {
363  little_vec_type X_dst = getLocalBlock (localRowIndex, colIndex);
364  const LO strideX = 1;
365  const_little_vec_type X_src (reinterpret_cast<const impl_scalar_type*> (vals),
366  getBlockSize (), strideX);
367  AXPY (STS::one (), X_src, X_dst);
368 }
369 
370 template<class Scalar, class LO, class GO, class Node>
371 bool
373 sumIntoLocalValues (const LO localRowIndex,
374  const LO colIndex,
375  const Scalar vals[]) const
376 {
377  if (! meshMap_.isNodeLocalElement (localRowIndex)) {
378  return false;
379  } else {
380  sumIntoLocalValuesImpl (localRowIndex, colIndex, vals);
381  return true;
382  }
383 }
384 
385 template<class Scalar, class LO, class GO, class Node>
386 bool
388 sumIntoGlobalValues (const GO globalRowIndex,
389  const LO colIndex,
390  const Scalar vals[]) const
391 {
392  const LO localRowIndex = meshMap_.getLocalElement (globalRowIndex);
393  if (localRowIndex == Teuchos::OrdinalTraits<LO>::invalid ()) {
394  return false;
395  } else {
396  sumIntoLocalValuesImpl (localRowIndex, colIndex, vals);
397  return true;
398  }
399 }
400 
401 template<class Scalar, class LO, class GO, class Node>
402 bool
404 getLocalRowView (const LO localRowIndex, const LO colIndex, Scalar*& vals) const
405 {
406  if (! meshMap_.isNodeLocalElement (localRowIndex)) {
407  return false;
408  } else {
409  little_vec_type X_ij = getLocalBlock (localRowIndex, colIndex);
410  vals = reinterpret_cast<Scalar*> (X_ij.ptr_on_device ());
411  return true;
412  }
413 }
414 
415 template<class Scalar, class LO, class GO, class Node>
416 bool
418 getGlobalRowView (const GO globalRowIndex, const LO colIndex, Scalar*& vals) const
419 {
420  const LO localRowIndex = meshMap_.getLocalElement (globalRowIndex);
421  if (localRowIndex == Teuchos::OrdinalTraits<LO>::invalid ()) {
422  return false;
423  } else {
424  little_vec_type X_ij = getLocalBlock (localRowIndex, colIndex);
425  vals = reinterpret_cast<Scalar*> (X_ij.ptr_on_device ());
426  return true;
427  }
428 }
429 
430 template<class Scalar, class LO, class GO, class Node>
433 getLocalBlock (const LO localRowIndex,
434  const LO colIndex) const
435 {
436  if (! isValidLocalMeshIndex (localRowIndex)) {
437  return little_vec_type (NULL, 0, 0);
438  } else {
439  const LO strideX = this->getStrideX ();
440  const size_t blockSize = getBlockSize ();
441  const size_t offset = colIndex * this->getStrideY () +
442  localRowIndex * blockSize * strideX;
443  impl_scalar_type* blockRaw = this->getRawPtr () + offset;
444  return little_vec_type (blockRaw, blockSize, strideX);
445  }
446 }
447 
448 template<class Scalar, class LO, class GO, class Node>
449 Teuchos::RCP<const typename BlockMultiVector<Scalar, LO, GO, Node>::mv_type>
452 {
453  using Teuchos::rcp;
454  using Teuchos::rcpFromRef;
455 
456  // The source object of an Import or Export must be a
457  // BlockMultiVector or MultiVector (a Vector is a MultiVector). Try
458  // them in that order; one must succeed. Note that the target of
459  // the Import or Export calls checkSizes in DistObject's doTransfer.
461  const this_type* srcBlkVec = dynamic_cast<const this_type*> (&src);
462  if (srcBlkVec == NULL) {
463  const mv_type* srcMultiVec = dynamic_cast<const mv_type*> (&src);
464  if (srcMultiVec == NULL) {
465  // FIXME (mfh 05 May 2014) Tpetra::MultiVector's operator=
466  // currently does a shallow copy by default. This is why we
467  // return by RCP, rather than by value.
468  return rcp (new mv_type ());
469  } else { // src is a MultiVector
470  return rcp (srcMultiVec, false); // nonowning RCP
471  }
472  } else { // src is a BlockMultiVector
473  return rcpFromRef (srcBlkVec->mv_); // nonowning RCP
474  }
475 }
476 
477 template<class Scalar, class LO, class GO, class Node>
480 {
481  return ! getMultiVectorFromSrcDistObject (src).is_null ();
482 }
483 
484 template<class Scalar, class LO, class GO, class Node>
487  size_t numSameIDs,
488  const Teuchos::ArrayView<const LO>& permuteToLIDs,
489  const Teuchos::ArrayView<const LO>& permuteFromLIDs)
490 {
491  const char tfecfFuncName[] = "copyAndPermute: ";
492  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
493  permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error,
494  "permuteToLIDs and permuteFromLIDs must have the same size."
495  << std::endl << "permuteToLIDs.size() = " << permuteToLIDs.size ()
496  << " != permuteFromLIDs.size() = " << permuteFromLIDs.size () << ".");
497 
499  Teuchos::RCP<const BMV> srcAsBmvPtr = getBlockMultiVectorFromSrcDistObject (src);
500  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
501  srcAsBmvPtr.is_null (), std::invalid_argument,
502  "The source of an Import or Export to a BlockMultiVector "
503  "must also be a BlockMultiVector.");
504  const BMV& srcAsBmv = *srcAsBmvPtr;
505 
506  // FIXME (mfh 23 Apr 2014) This implementation is sequential and
507  // assumes UVM.
508 
509  const LO numVecs = getNumVectors ();
510  const LO numSame = static_cast<LO> (numSameIDs);
511  for (LO j = 0; j < numVecs; ++j) {
512  for (LO lclRow = 0; lclRow < numSame; ++lclRow) {
513  deep_copy (getLocalBlock (lclRow, j), srcAsBmv.getLocalBlock (lclRow, j));
514  }
515  }
516 
517  // FIXME (mfh 20 June 2012) For an Export with duplicate GIDs on the
518  // same process, this merges their values by replacement of the last
519  // encountered GID, not by the specified merge rule (such as ADD).
520  const LO numPermuteLIDs = static_cast<LO> (permuteToLIDs.size ());
521  for (LO j = 0; j < numVecs; ++j) {
522  for (LO k = numSame; k < numPermuteLIDs; ++k) {
523  deep_copy (getLocalBlock (permuteToLIDs[k], j), srcAsBmv.getLocalBlock (permuteFromLIDs[k], j));
524  }
525  }
526 }
527 
528 template<class Scalar, class LO, class GO, class Node>
531  const Teuchos::ArrayView<const LO>& exportLIDs,
532  Teuchos::Array<impl_scalar_type>& exports,
533  const Teuchos::ArrayView<size_t>& /* numPacketsPerLID */,
534  size_t& constantNumPackets,
535  Tpetra::Distributor& /* distor */)
536 {
538  typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
539  const char tfecfFuncName[] = "packAndPrepare: ";
540 
541  Teuchos::RCP<const BMV> srcAsBmvPtr = getBlockMultiVectorFromSrcDistObject (src);
542  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
543  srcAsBmvPtr.is_null (), std::invalid_argument,
544  "The source of an Import or Export to a BlockMultiVector "
545  "must also be a BlockMultiVector.");
546  const BMV& srcAsBmv = *srcAsBmvPtr;
547 
548  const LO numVecs = getNumVectors ();
549  const LO blockSize = getBlockSize ();
550 
551  // Number of things to pack per LID is the block size, times the
552  // number of columns. Input LIDs correspond to the mesh points, not
553  // the degrees of freedom (DOFs).
554  constantNumPackets =
555  static_cast<size_t> (blockSize) * static_cast<size_t> (numVecs);
556  const size_type numMeshLIDs = exportLIDs.size ();
557 
558  const size_type requiredExportsSize = numMeshLIDs *
559  static_cast<size_type> (blockSize) * static_cast<size_type> (numVecs);
560  exports.resize (requiredExportsSize);
561 
562  try {
563  size_type curExportPos = 0;
564  for (size_type meshLidIndex = 0; meshLidIndex < numMeshLIDs; ++meshLidIndex) {
565  for (LO j = 0; j < numVecs; ++j, curExportPos += blockSize) {
566  const LO meshLid = exportLIDs[meshLidIndex];
567  impl_scalar_type* const curExportPtr = &exports[curExportPos];
568  little_vec_type X_dst (curExportPtr, blockSize, 1);
569  little_vec_type X_src = srcAsBmv.getLocalBlock (meshLid, j);
570 
571  deep_copy (X_dst, X_src);
572  }
573  }
574  } catch (std::exception& e) {
575  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
576  true, std::logic_error, "Oh no! packAndPrepare on Process "
577  << meshMap_.getComm ()->getRank () << " raised the following exception: "
578  << e.what ());
579  }
580 }
581 
582 template<class Scalar, class LO, class GO, class Node>
584 unpackAndCombine (const Teuchos::ArrayView<const LO>& importLIDs,
585  const Teuchos::ArrayView<const impl_scalar_type>& imports,
586  const Teuchos::ArrayView<size_t>& numPacketsPerLID,
587  size_t constantNumPackets,
588  Tpetra::Distributor& distor,
590 {
591  typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
592  const char tfecfFuncName[] = "unpackAndCombine: ";
593 
594  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
595  CM != ADD && CM != REPLACE && CM != INSERT && CM != ABSMAX && CM != ZERO,
596  std::invalid_argument, "Invalid CombineMode: " << CM << ". Valid "
597  "CombineMode values are ADD, REPLACE, INSERT, ABSMAX, and ZERO.");
598 
599  if (CM == ZERO) {
600  return; // Combining does nothing, so we don't have to combine anything.
601  }
602 
603  // Number of things to pack per LID is the block size.
604  // Input LIDs correspond to the mesh points, not the DOFs.
605  const size_type numMeshLIDs = importLIDs.size ();
606  const LO blockSize = getBlockSize ();
607  const LO numVecs = getNumVectors ();
608 
609  const size_type requiredImportsSize = numMeshLIDs *
610  static_cast<size_type> (blockSize) * static_cast<size_type> (numVecs);
611  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
612  imports.size () < requiredImportsSize, std::logic_error,
613  ": imports.size () = " << imports.size ()
614  << " < requiredImportsSize = " << requiredImportsSize << ".");
615 
616  size_type curImportPos = 0;
617  for (size_type meshLidIndex = 0; meshLidIndex < numMeshLIDs; ++meshLidIndex) {
618  for (LO j = 0; j < numVecs; ++j, curImportPos += blockSize) {
619  const LO meshLid = importLIDs[meshLidIndex];
620  const impl_scalar_type* const curImportPtr = &imports[curImportPos];
621 
622  const_little_vec_type X_src (curImportPtr, blockSize, 1);
623  little_vec_type X_dst = getLocalBlock (meshLid, j);
624 
625  if (CM == INSERT || CM == REPLACE) {
626  deep_copy (X_dst, X_src);
627  } else if (CM == ADD) {
628  AXPY (STS::one (), X_src, X_dst);
629  } else if (CM == ABSMAX) {
630  Impl::absMax (X_dst, X_src);
631  }
632  }
633  }
634 }
635 
636 template<class Scalar, class LO, class GO, class Node>
638 putScalar (const Scalar& val)
639 {
640  getMultiVectorView ().putScalar (val);
641 }
642 
643 template<class Scalar, class LO, class GO, class Node>
645 scale (const Scalar& val)
646 {
647  getMultiVectorView ().scale (val);
648 }
649 
650 } // namespace Experimental
651 } // namespace Tpetra
652 
653 //
654 // Explicit instantiation macro
655 //
656 // Must be expanded from within the Tpetra namespace!
657 //
658 #define TPETRA_EXPERIMENTAL_BLOCKMULTIVECTOR_INSTANT(S,LO,GO,NODE) \
659  template class Experimental::BlockMultiVector< S, LO, GO, NODE >;
660 
661 #endif // TPETRA_EXPERIMENTAL_BLOCKMULTIVECTOR_DEF_HPP
void scale(const Scalar &alpha)
Scale in place: this = alpha*this.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
mv_type mv_
The Tpetra::MultiVector used to represent the data.
bool isContiguous() const
True if this Map is distributed contiguously, else false.
void putScalar(const Scalar &val)
Fill all entries with the given value val.
mv_type getMultiVectorView()
Get a Tpetra::MultiVector that views this BlockMultiVector&#39;s data.
Tpetra::Map< LO, GO, Node > map_type
The specialization of Tpetra::Map that this class uses.
bool replaceGlobalValues(const GO globalRowIndex, const LO colIndex, const Scalar vals[]) const
Replace all values at the given mesh point, using a global index.
void setCopyOrView(const Teuchos::DataAccess copyOrView)
Set whether this has copy (copyOrView = Teuchos::Copy) or view (copyOrView = Teuchos::View) semantics...
void AXPY(const CoefficientType &alpha, const ViewType1 &x, const ViewType2 &y)
y := y + alpha * x (dense vector or matrix update)
One or more distributed dense vectors.
bool sumIntoGlobalValues(const GO globalRowIndex, const LO colIndex, const Scalar vals[]) const
Sum into all values at the given mesh point, using a global index.
MultiVector for multiple degrees of freedom per mesh point.
GlobalOrdinal getIndexBase() const
The index base for this Map.
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.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
virtual void unpackAndCombine(const Teuchos::ArrayView< const LO > &importLIDs, const Teuchos::ArrayView< const impl_scalar_type > &imports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t constantNumPackets, Tpetra::Distributor &distor, Tpetra::CombineMode CM)
Perform any unpacking and combining after communication.
void deep_copy(const LittleBlock< ST2, LO > &dst, const LittleBlock< ST1, LO > &src, typename std::enable_if< std::is_convertible< ST1, ST2 >::value &&!std::is_const< ST2 >::value, int >::type *=NULL)
Copy the LittleBlock src into the LittleBlock dst.
void scale(const Scalar &val)
Multiply all entries in place by the given value val.
little_vec_type getLocalBlock(const LO localRowIndex, const LO colIndex) const
Get a view of the degrees of freedom at the given mesh point.
bool isNodeLocalElement(LocalOrdinal localIndex) const
Whether the given local index is valid for this Map on the calling process.
LO getBlockSize() const
Get the number of degrees of freedom per mesh point.
size_t global_size_t
Global size_t object.
Insert new values that don&#39;t currently exist.
impl_scalar_type * getRawPtr() const
Raw pointer to the MultiVector&#39;s data.
Nonowning view of a set of degrees of freedom corresponding to a mesh point in a block vector or mult...
LO getNumVectors() const
Get the number of columns (vectors) in the BlockMultiVector.
Sets up and executes a communication plan for a Tpetra DistObject.
CombineMode
Rule for combining data in an Import or Export.
Sum new values into existing values.
Teuchos::RCP< Node > getNode() const
Get this Map&#39;s Node object.
bool getLocalRowView(const LO localRowIndex, const LO colIndex, Scalar *&vals) const
Get a writeable view of the entries at the given mesh point, using a local index. ...
Replace old value with maximum of magnitudes of old and new values.
Abstract base class for objects that can be the source of an Import or Export operation.
bool isValidLocalMeshIndex(const LO meshLocalIndex) const
True if and only if meshLocalIndex is a valid local index in the mesh Map.
Scalar * ptr_on_device() const
Pointer to the vector&#39;s entries.
Replace existing values with new values.
Replace old values with zero.
Teuchos::DataAccess getCopyOrView() const
Get whether this has copy (copyOrView = Teuchos::Copy) or view (copyOrView = Teuchos::View) semantics...
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
Kokkos::Details::ArithTraits< Scalar >::val_type impl_scalar_type
The type used internally in place of Scalar.
LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const
The local index corresponding to the given global index.
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a view of the global indices owned by this process.
virtual bool checkSizes(const Tpetra::SrcDistObject &source)
Compare the source and target (this) objects for compatibility.
size_t getNodeNumElements() const
The number of elements belonging to the calling process.
LittleVector< impl_scalar_type, LO > little_vec_type
"Block view" of all degrees of freedom at a mesh point, for a single column of the MultiVector...
size_t getStrideY() const
Stride between consecutive local entries in the same row.
virtual void copyAndPermute(const Tpetra::SrcDistObject &source, size_t numSameIDs, const Teuchos::ArrayView< const LO > &permuteToLIDs, const Teuchos::ArrayView< const LO > &permuteFromLIDs)
Perform copies and permutations that are local to this process.
size_t getNumVectors() const
Number of columns in the multivector.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
bool replaceLocalValues(const LO localRowIndex, const LO colIndex, const Scalar vals[]) const
Replace all values at the given mesh point, using local row and column indices.
Teuchos::RCP< const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > subView(const Teuchos::Range1D &colRng) const
Return a const MultiVector with const views of selected columns.
size_t getStrideX() const
Stride between consecutive local entries in the same column.
virtual void packAndPrepare(const Tpetra::SrcDistObject &source, const Teuchos::ArrayView< const LO > &exportLIDs, Teuchos::Array< impl_scalar_type > &exports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t &constantNumPackets, Tpetra::Distributor &distor)
Perform any packing or preparation required for communication.
bool sumIntoLocalValues(const LO localRowIndex, const LO colIndex, const Scalar vals[]) const
Sum into all values at the given mesh point, using a local index.
global_size_t getGlobalNumElements() const
The number of elements in this Map.
bool getGlobalRowView(const GO globalRowIndex, const LO colIndex, Scalar *&vals) const
Get a writeable view of the entries at the given mesh point, using a global index.