Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Belos_TpetraAdapter_MP_Vector.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef BELOS_TPETRA_ADAPTER_MP_VECTOR_HPP
43#define BELOS_TPETRA_ADAPTER_MP_VECTOR_HPP
44
45#include "BelosTpetraAdapter.hpp"
48#include "KokkosBlas.hpp"
49
50#ifdef HAVE_BELOS_TSQR
52#endif // HAVE_BELOS_TSQR
53
54namespace Belos {
55
57 //
58 // Implementation of Belos::MultiVecTraits for Tpetra::MultiVector.
59 //
61
72 template<class Storage, class LO, class GO, class Node>
73 class MultiVecTraits<typename Storage::value_type,
74 Tpetra::MultiVector< Sacado::MP::Vector<Storage>,
75 LO, GO, Node > > {
76 public:
77 typedef typename Storage::ordinal_type s_ordinal;
78 typedef typename Storage::value_type BaseScalar;
80 typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
81 public:
82 typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::dot_type dot_type;
83 typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::mag_type mag_type;
84
90 static Teuchos::RCP<MV> Clone (const MV& X, const int numVecs) {
91 Teuchos::RCP<MV> Y (new MV (X.getMap (), numVecs, false));
92 Y->setCopyOrView (Teuchos::View);
93 return Y;
94 }
95
97 static Teuchos::RCP<MV> CloneCopy (const MV& X)
98 {
99 // Make a deep copy of X. The one-argument copy constructor
100 // does a shallow copy by default; the second argument tells it
101 // to do a deep copy.
102 Teuchos::RCP<MV> X_copy (new MV (X, Teuchos::Copy));
103 // Make Tpetra::MultiVector use the new view semantics. This is
104 // a no-op for the Kokkos refactor version of Tpetra; it only
105 // does something for the "classic" version of Tpetra. This
106 // shouldn't matter because Belos only handles MV through RCP
107 // and through this interface anyway, but it doesn't hurt to set
108 // it and make sure that it works.
109 X_copy->setCopyOrView (Teuchos::View);
110 return X_copy;
111 }
112
124 static Teuchos::RCP<MV>
125 CloneCopy (const MV& mv, const std::vector<int>& index)
126 {
127#ifdef HAVE_TPETRA_DEBUG
128 const char fnName[] = "Belos::MultiVecTraits::CloneCopy(mv,index)";
129 const size_t inNumVecs = mv.getNumVectors ();
130 TEUCHOS_TEST_FOR_EXCEPTION(
131 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
132 std::runtime_error, fnName << ": All indices must be nonnegative.");
133 TEUCHOS_TEST_FOR_EXCEPTION(
134 index.size () > 0 &&
135 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
136 std::runtime_error,
137 fnName << ": All indices must be strictly less than the number of "
138 "columns " << inNumVecs << " of the input multivector mv.");
139#endif // HAVE_TPETRA_DEBUG
140
141 // Tpetra wants an array of size_t, not of int.
142 Teuchos::Array<size_t> columns (index.size ());
143 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
144 columns[j] = index[j];
145 }
146 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
147 // continuous column index range in MultiVector::subCopy, so we
148 // don't have to check here.
149 Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
150 X_copy->setCopyOrView (Teuchos::View);
151 return X_copy;
152 }
153
160 static Teuchos::RCP<MV>
161 CloneCopy (const MV& mv, const Teuchos::Range1D& index)
162 {
163 const bool validRange = index.size() > 0 &&
164 index.lbound() >= 0 &&
165 index.ubound() < GetNumberVecs(mv);
166 if (! validRange) { // invalid range; generate error message
167 std::ostringstream os;
168 os << "Belos::MultiVecTraits::CloneCopy(mv,index=["
169 << index.lbound() << "," << index.ubound() << "]): ";
170 TEUCHOS_TEST_FOR_EXCEPTION(
171 index.size() == 0, std::invalid_argument,
172 os.str() << "Empty index range is not allowed.");
173 TEUCHOS_TEST_FOR_EXCEPTION(
174 index.lbound() < 0, std::invalid_argument,
175 os.str() << "Index range includes negative index/ices, which is not "
176 "allowed.");
177 TEUCHOS_TEST_FOR_EXCEPTION(
178 index.ubound() >= GetNumberVecs(mv), std::invalid_argument,
179 os.str() << "Index range exceeds number of vectors "
180 << mv.getNumVectors() << " in the input multivector.");
181 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
182 os.str() << "Should never get here!");
183 }
184 Teuchos::RCP<MV> X_copy = mv.subCopy (index);
185 X_copy->setCopyOrView (Teuchos::View);
186 return X_copy;
187 }
188
189
190 static Teuchos::RCP<MV>
191 CloneViewNonConst (MV& mv, const std::vector<int>& index)
192 {
193#ifdef HAVE_TPETRA_DEBUG
194 const char fnName[] = "Belos::MultiVecTraits::CloneViewNonConst(mv,index)";
195 const size_t numVecs = mv.getNumVectors ();
196 TEUCHOS_TEST_FOR_EXCEPTION(
197 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
198 std::invalid_argument,
199 fnName << ": All indices must be nonnegative.");
200 TEUCHOS_TEST_FOR_EXCEPTION(
201 index.size () > 0 &&
202 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
203 std::invalid_argument,
204 fnName << ": All indices must be strictly less than the number of "
205 "columns " << numVecs << " in the input MultiVector mv.");
206#endif // HAVE_TPETRA_DEBUG
207
208 // Tpetra wants an array of size_t, not of int.
209 Teuchos::Array<size_t> columns (index.size ());
210 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
211 columns[j] = index[j];
212 }
213 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
214 // continuous column index range in
215 // MultiVector::subViewNonConst, so we don't have to check here.
216 Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
217 X_view->setCopyOrView (Teuchos::View);
218 return X_view;
219 }
220
221 static Teuchos::RCP<MV>
222 CloneViewNonConst (MV& mv, const Teuchos::Range1D& index)
223 {
224 // NOTE (mfh 11 Jan 2011) We really should check for possible
225 // overflow of int here. However, the number of columns in a
226 // multivector typically fits in an int.
227 const int numCols = static_cast<int> (mv.getNumVectors());
228 const bool validRange = index.size() > 0 &&
229 index.lbound() >= 0 && index.ubound() < numCols;
230 if (! validRange) {
231 std::ostringstream os;
232 os << "Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
233 << index.lbound() << ", " << index.ubound() << "]): ";
234 TEUCHOS_TEST_FOR_EXCEPTION(
235 index.size() == 0, std::invalid_argument,
236 os.str() << "Empty index range is not allowed.");
237 TEUCHOS_TEST_FOR_EXCEPTION(
238 index.lbound() < 0, std::invalid_argument,
239 os.str() << "Index range includes negative inde{x,ices}, which is "
240 "not allowed.");
241 TEUCHOS_TEST_FOR_EXCEPTION(
242 index.ubound() >= numCols, std::invalid_argument,
243 os.str() << "Index range exceeds number of vectors " << numCols
244 << " in the input multivector.");
245 TEUCHOS_TEST_FOR_EXCEPTION(
246 true, std::logic_error,
247 os.str() << "Should never get here!");
248 }
249 Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
250 X_view->setCopyOrView (Teuchos::View);
251 return X_view;
252 }
253
254 static Teuchos::RCP<const MV>
255 CloneView (const MV& mv, const std::vector<int>& index)
256 {
257#ifdef HAVE_TPETRA_DEBUG
258 const char fnName[] = "Belos::MultiVecTraits<Scalar, "
259 "Tpetra::MultiVector<...> >::CloneView(mv,index)";
260 const size_t numVecs = mv.getNumVectors ();
261 TEUCHOS_TEST_FOR_EXCEPTION(
262 *std::min_element (index.begin (), index.end ()) < 0,
263 std::invalid_argument,
264 fnName << ": All indices must be nonnegative.");
265 TEUCHOS_TEST_FOR_EXCEPTION(
266 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
267 std::invalid_argument,
268 fnName << ": All indices must be strictly less than the number of "
269 "columns " << numVecs << " in the input MultiVector mv.");
270#endif // HAVE_TPETRA_DEBUG
271
272 // Tpetra wants an array of size_t, not of int.
273 Teuchos::Array<size_t> columns (index.size ());
274 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
275 columns[j] = index[j];
276 }
277 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
278 // continuous column index range in MultiVector::subView, so we
279 // don't have to check here.
280 Teuchos::RCP<const MV> X_view = mv.subView (columns);
281 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
282 return X_view;
283 }
284
285 static Teuchos::RCP<const MV>
286 CloneView (const MV& mv, const Teuchos::Range1D& index)
287 {
288 // NOTE (mfh 11 Jan 2011) We really should check for possible
289 // overflow of int here. However, the number of columns in a
290 // multivector typically fits in an int.
291 const int numCols = static_cast<int> (mv.getNumVectors());
292 const bool validRange = index.size() > 0 &&
293 index.lbound() >= 0 && index.ubound() < numCols;
294 if (! validRange) {
295 std::ostringstream os;
296 os << "Belos::MultiVecTraits::CloneView(mv, index=["
297 << index.lbound () << ", " << index.ubound() << "]): ";
298 TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
299 os.str() << "Empty index range is not allowed.");
300 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
301 os.str() << "Index range includes negative index/ices, which is not "
302 "allowed.");
303 TEUCHOS_TEST_FOR_EXCEPTION(
304 index.ubound() >= numCols, std::invalid_argument,
305 os.str() << "Index range exceeds number of vectors " << numCols
306 << " in the input multivector.");
307 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
308 os.str() << "Should never get here!");
309 }
310 Teuchos::RCP<const MV> X_view = mv.subView (index);
311 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
312 return X_view;
313 }
314
315 static ptrdiff_t GetGlobalLength (const MV& mv) {
316 return static_cast<ptrdiff_t> (mv.getGlobalLength ());
317 }
318
319 static int GetNumberVecs (const MV& mv) {
320 return static_cast<int> (mv.getNumVectors ());
321 }
322
323 static bool HasConstantStride (const MV& mv) {
324 return mv.isConstantStride ();
325 }
326
327 template <class DstType, class SrcType>
329 const DstType& dst,
330 const SrcType& src
331 )
332 {
333 if (std::is_same<typename SrcType::memory_space, Kokkos::HostSpace>::value)
334 {
335 auto dst_i = Kokkos::create_mirror(dst);
336 for (size_t i = 0; i < src.extent(0); i++)
337 {
338 for (size_t j = 0; j < src.extent(1); j++)
339 {
340 dst_i(i, j) = src(i, j);
341 }
342 }
343 Kokkos::deep_copy(dst, dst_i);
344 }
345 else
346 {
347 auto src_i = Kokkos::create_mirror(src);
348 Kokkos::deep_copy(src_i, src);
349 for (size_t i = 0; i < src.extent(0); i++)
350 {
351 for (size_t j = 0; j < src.extent(1); j++)
352 {
353 dst(i, j) = src_i(i, j);
354 }
355 }
356 }
357 }
358
359 static void
361 const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
362 const Teuchos::SerialDenseMatrix<int,dot_type>& B,
363 const dot_type& beta,
364 Tpetra::MultiVector<Scalar,LO,GO,Node>& C)
365 {
366 using Teuchos::RCP;
367 using Teuchos::rcp;
368
369 // Check if numRowsB == numColsB == 1, in which case we can call update()
370 const int numRowsB = B.numRows ();
371 const int numColsB = B.numCols ();
372 const int strideB = B.stride ();
373 if (numRowsB == 1 && numColsB == 1) {
374 C.update (alpha*B(0,0), A, beta);
375 return;
376 }
377
378 // Ensure A and C have constant stride
379 RCP<const MV> Atmp;
380 RCP< MV> Ctmp;
381 if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
382 else Atmp = rcp(&A,false);
383
384 if (C.isConstantStride() == false) Ctmp = rcp (new MV (C, Teuchos::Copy));
385 else Ctmp = rcp(&C,false);
386
387 // Create flattened view's
388 typedef Tpetra::MultiVector<dot_type,LO,GO,Node> FMV;
389 typedef typename FMV::dual_view_type::t_dev flat_view_type;
390 typedef typename flat_view_type::execution_space execution_space;
391 typename flat_view_type::const_type flat_A_view = Atmp->getLocalViewDevice(Tpetra::Access::ReadOnly);
392 flat_view_type flat_C_view = Ctmp->getLocalViewDevice(Tpetra::Access::OverwriteAll);
393
394 // Create a view for B on the host
395 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> b_host_view_type;
396 b_host_view_type B_view_host_input( B.values(), strideB, numColsB);
397 auto B_view_host = Kokkos::subview(B_view_host_input,
398 Kokkos::pair<int,int>(0,numRowsB),
399 Kokkos::pair<int,int>(0,numColsB));
400
401 // Create view for B on the device -- need to be careful to get the
402 // right stride to match B
403 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> b_view_type;
404 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> b_1d_view_type;
405 b_1d_view_type B_1d_view_dev(Kokkos::ViewAllocateWithoutInitializing("B"), numRowsB*numColsB);
406 b_view_type B_view_dev( B_1d_view_dev.data(), numRowsB, numColsB);
407
408 if (Kokkos::SpaceAccessibility<Kokkos::HostSpace, typename execution_space::memory_space>::accessible)
409 {
410 Kokkos::deep_copy(B_view_dev, B_view_host);
411 }
412 else
413 {
414 deep_copy_2d_view_with_intercessory_space(B_view_dev, B_view_host);
415 }
416
417 // Do local multiply
418 {
419 const char ctransA = 'N', ctransB = 'N';
421 &ctransA, &ctransB,
422 alpha, flat_A_view, B_view_dev, beta, flat_C_view);
423 }
424 // Copy back to C if we made a copy
425 if (C.isConstantStride() == false)
426 C.assign(*Ctmp);
427 }
428
436 static void
438 const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
439 Scalar beta,
440 const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
441 Tpetra::MultiVector<Scalar,LO,GO,Node>& mv)
442 {
443 mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
444 }
445
446 static void
447 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
448 const Scalar& alpha)
449 {
450 mv.scale (alpha);
451 }
452
453 static void
454 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
455 const std::vector<BaseScalar>& alphas)
456 {
457 std::vector<Scalar> alphas_mp(alphas.size());
458 const size_t sz = alphas.size();
459 for (size_t i=0; i<sz; ++i)
460 alphas_mp[i] = alphas[i];
461 mv.scale (alphas_mp);
462 }
463
464 static void
465 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
466 const std::vector<Scalar>& alphas)
467 {
468 mv.scale (alphas);
469 }
470
471 static void
473 const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
474 const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
475 Teuchos::SerialDenseMatrix<int,dot_type>& C)
476 {
477 using Teuchos::Comm;
478 using Teuchos::RCP;
479 using Teuchos::rcp;
480 using Teuchos::REDUCE_SUM;
481 using Teuchos::reduceAll;
482
483 // Check if numRowsC == numColsC == 1, in which case we can call dot()
484 const int numRowsC = C.numRows ();
485 const int numColsC = C.numCols ();
486 const int strideC = C.stride ();
487 if (numRowsC == 1 && numColsC == 1) {
488 if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
489 // Short-circuit, as required by BLAS semantics.
490 C(0,0) = alpha;
491 return;
492 }
493 A.dot (B, Teuchos::ArrayView<dot_type> (C.values (), 1));
494 if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
495 C(0,0) *= alpha;
496 }
497 return;
498 }
499
500 // Ensure A and B have constant stride
501 RCP<const MV> Atmp, Btmp;
502 if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
503 else Atmp = rcp(&A,false);
504
505 if (B.isConstantStride() == false) Btmp = rcp (new MV (B, Teuchos::Copy));
506 else Btmp = rcp(&B,false);
507
508 // Create flattened Kokkos::MultiVector's
509 typedef Tpetra::MultiVector<dot_type,LO,GO,Node> FMV;
510 typedef typename FMV::dual_view_type::t_dev flat_view_type;
511 typedef typename flat_view_type::execution_space execution_space;
512 typename flat_view_type::const_type flat_A_view = Atmp->getLocalViewDevice(Tpetra::Access::ReadOnly);
513 typename flat_view_type::const_type flat_B_view = Btmp->getLocalViewDevice(Tpetra::Access::ReadOnly);
514
515 // Create a view for C on the host
516 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> c_host_view_type;
517 c_host_view_type C_view_host_input( C.values(), strideC, numColsC);
518 auto C_view_host = Kokkos::subview(C_view_host_input,
519 Kokkos::pair<int,int>(0,numRowsC),
520 Kokkos::pair<int,int>(0,numColsC));
521
522 // Create view for C on the device -- need to be careful to get the
523 // right stride to match C (allow setting to 0 for first-touch)
524 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> c_view_type;
525 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> c_1d_view_type;
526 c_1d_view_type C_1d_view_dev("C", numRowsC*numColsC);
527 c_view_type C_view_dev( C_1d_view_dev.data(), numRowsC, numColsC);
528
529 // Do local multiply
530 {
531 const char ctransA = 'C', ctransB = 'N';
533 &ctransA, &ctransB,
534 alpha, flat_A_view, flat_B_view,
535 Kokkos::Details::ArithTraits<dot_type>::zero(),
536 C_view_dev);
537 }
538
539 // reduce across processors -- could check for RDMA
540 RCP<const Comm<int> > pcomm = A.getMap()->getComm ();
541 if (pcomm->getSize () == 1)
542 {
543 if (Kokkos::SpaceAccessibility<execution_space, Kokkos::HostSpace>::accessible)
544 {
545 Kokkos::deep_copy(C_view_host, C_view_dev);
546 }
547 else
548 {
549 deep_copy_2d_view_with_intercessory_space(C_view_host, C_view_dev);
550 }
551 }
552 else
553 {
554 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, Kokkos::HostSpace> c_1d_host_view_type;
555 c_1d_host_view_type C_1d_view_tmp(Kokkos::ViewAllocateWithoutInitializing("C_tmp"), strideC*numColsC);
556 c_host_view_type C_view_tmp_input( C_1d_view_tmp.data(), strideC, numColsC);
557 auto C_view_tmp = Kokkos::subview(C_view_tmp_input,
558 Kokkos::pair<int,int>(0,numRowsC),
559 Kokkos::pair<int,int>(0,numColsC));
560
561 if (Kokkos::SpaceAccessibility<execution_space, Kokkos::HostSpace>::accessible)
562 {
563 Kokkos::deep_copy(C_view_tmp, C_view_dev);
564 }
565 else
566 {
567 deep_copy_2d_view_with_intercessory_space(C_view_tmp, C_view_dev);
568 }
569
570 reduceAll<int> (*pcomm, REDUCE_SUM, strideC*numColsC,
571 C_view_tmp.data(),
572 C_view_host.data());
573 }
574 }
575
577 static void
578 MvDot (const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
579 const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
580 std::vector<dot_type>& dots)
581 {
582 const size_t numVecs = A.getNumVectors ();
583
584 TEUCHOS_TEST_FOR_EXCEPTION(
585 numVecs != B.getNumVectors (), std::invalid_argument,
586 "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
587 "A and B must have the same number of columns. "
588 "A has " << numVecs << " column(s), "
589 "but B has " << B.getNumVectors () << " column(s).");
590#ifdef HAVE_TPETRA_DEBUG
591 TEUCHOS_TEST_FOR_EXCEPTION(
592 dots.size() < numVecs, std::invalid_argument,
593 "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
594 "The output array 'dots' must have room for all dot products. "
595 "A and B each have " << numVecs << " column(s), "
596 "but 'dots' only has " << dots.size() << " entry(/ies).");
597#endif // HAVE_TPETRA_DEBUG
598
599 Teuchos::ArrayView<dot_type> av (dots);
600 A.dot (B, av (0, numVecs));
601 }
602
604 static void
605 MvNorm (const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
606 std::vector<mag_type> &normvec,
607 NormType type=TwoNorm)
608 {
609
610#ifdef HAVE_TPETRA_DEBUG
611 typedef std::vector<int>::size_type size_type;
612 TEUCHOS_TEST_FOR_EXCEPTION(
613 normvec.size () < static_cast<size_type> (mv.getNumVectors ()),
614 std::invalid_argument,
615 "Belos::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
616 "argument must have at least as many entries as the number of vectors "
617 "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
618 << " < mv.getNumVectors() = " << mv.getNumVectors () << ".");
619#endif
620
621 Teuchos::ArrayView<mag_type> av(normvec);
622 switch (type) {
623 case OneNorm:
624 mv.norm1(av(0,mv.getNumVectors()));
625 break;
626 case TwoNorm:
627 mv.norm2(av(0,mv.getNumVectors()));
628 break;
629 case InfNorm:
630 mv.normInf(av(0,mv.getNumVectors()));
631 break;
632 default:
633 // Throw logic_error rather than invalid_argument, because if
634 // we get here, it's probably the fault of a Belos solver,
635 // rather than a user giving Belos an invalid input.
636 TEUCHOS_TEST_FOR_EXCEPTION(
637 true, std::logic_error,
638 "Belos::MultiVecTraits::MvNorm: Invalid NormType value " << type
639 << ". Valid values are OneNorm=" << OneNorm << ", TwoNorm="
640 << TwoNorm <<", and InfNorm=" << InfNorm << ". If you are a Belos "
641 "user and have not modified Belos in any way, and you get this "
642 "message, then this is probably a bug in the Belos solver you were "
643 "using. Please report this to the Belos developers.");
644 }
645 }
646
647 static void
648 SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
649 {
650 using Teuchos::Range1D;
651 using Teuchos::RCP;
652 const size_t inNumVecs = A.getNumVectors ();
653#ifdef HAVE_TPETRA_DEBUG
654 TEUCHOS_TEST_FOR_EXCEPTION(
655 inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
656 "Belos::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
657 "have no more entries as the number of columns in the input MultiVector"
658 " A. A.getNumVectors() = " << inNumVecs << " < index.size () = "
659 << index.size () << ".");
660#endif // HAVE_TPETRA_DEBUG
661 RCP<MV> mvsub = CloneViewNonConst (mv, index);
662 if (inNumVecs > static_cast<size_t> (index.size ())) {
663 RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
664 ::Tpetra::deep_copy (*mvsub, *Asub);
665 } else {
666 ::Tpetra::deep_copy (*mvsub, A);
667 }
668 }
669
670 static void
671 SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
672 {
673 // Range1D bounds are signed; size_t is unsigned.
674 // Assignment of Tpetra::MultiVector is a deep copy.
675
676 // Tpetra::MultiVector::getNumVectors() returns size_t. It's
677 // fair to assume that the number of vectors won't overflow int,
678 // since the typical use case of multivectors involves few
679 // columns, but it's friendly to check just in case.
680 const size_t maxInt =
681 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
682 const bool overflow =
683 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
684 if (overflow) {
685 std::ostringstream os;
686 os << "Belos::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
687 << ", " << index.ubound () << "], mv): ";
688 TEUCHOS_TEST_FOR_EXCEPTION(
689 maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
690 "of columns (size_t) in the input MultiVector 'A' overflows int.");
691 TEUCHOS_TEST_FOR_EXCEPTION(
692 maxInt < mv.getNumVectors (), std::range_error, os.str () << "Number "
693 "of columns (size_t) in the output MultiVector 'mv' overflows int.");
694 }
695 // We've already validated the static casts above.
696 const int numColsA = static_cast<int> (A.getNumVectors ());
697 const int numColsMv = static_cast<int> (mv.getNumVectors ());
698 // 'index' indexes into mv; it's the index set of the target.
699 const bool validIndex =
700 index.lbound () >= 0 && index.ubound () < numColsMv;
701 // We can't take more columns out of A than A has.
702 const bool validSource = index.size () <= numColsA;
703
704 if (! validIndex || ! validSource) {
705 std::ostringstream os;
706 os << "Belos::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
707 << ", " << index.ubound () << "], mv): ";
708 TEUCHOS_TEST_FOR_EXCEPTION(
709 index.lbound() < 0, std::invalid_argument,
710 os.str() << "Range lower bound must be nonnegative.");
711 TEUCHOS_TEST_FOR_EXCEPTION(
712 index.ubound() >= numColsMv, std::invalid_argument,
713 os.str() << "Range upper bound must be less than the number of "
714 "columns " << numColsA << " in the 'mv' output argument.");
715 TEUCHOS_TEST_FOR_EXCEPTION(
716 index.size() > numColsA, std::invalid_argument,
717 os.str() << "Range must have no more elements than the number of "
718 "columns " << numColsA << " in the 'A' input argument.");
719 TEUCHOS_TEST_FOR_EXCEPTION(
720 true, std::logic_error, "Should never get here!");
721 }
722
723 // View of the relevant column(s) of the target multivector mv.
724 // We avoid view creation overhead by only creating a view if
725 // the index range is different than [0, (# columns in mv) - 1].
726 Teuchos::RCP<MV> mv_view;
727 if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
728 mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
729 } else {
730 mv_view = CloneViewNonConst (mv, index);
731 }
732
733 // View of the relevant column(s) of the source multivector A.
734 // If A has fewer columns than mv_view, then create a view of
735 // the first index.size() columns of A.
736 Teuchos::RCP<const MV> A_view;
737 if (index.size () == numColsA) {
738 A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
739 } else {
740 A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
741 }
742
743 ::Tpetra::deep_copy (*mv_view, *A_view);
744 }
745
746 static void Assign (const MV& A, MV& mv)
747 {
748 const char errPrefix[] = "Belos::MultiVecTraits::Assign(A, mv): ";
749
750 // Range1D bounds are signed; size_t is unsigned.
751 // Assignment of Tpetra::MultiVector is a deep copy.
752
753 // Tpetra::MultiVector::getNumVectors() returns size_t. It's
754 // fair to assume that the number of vectors won't overflow int,
755 // since the typical use case of multivectors involves few
756 // columns, but it's friendly to check just in case.
757 const size_t maxInt =
758 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
759 const bool overflow =
760 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
761 if (overflow) {
762 TEUCHOS_TEST_FOR_EXCEPTION(
763 maxInt < A.getNumVectors(), std::range_error,
764 errPrefix << "Number of columns in the input multivector 'A' "
765 "(a size_t) overflows int.");
766 TEUCHOS_TEST_FOR_EXCEPTION(
767 maxInt < mv.getNumVectors(), std::range_error,
768 errPrefix << "Number of columns in the output multivector 'mv' "
769 "(a size_t) overflows int.");
770 TEUCHOS_TEST_FOR_EXCEPTION(
771 true, std::logic_error, "Should never get here!");
772 }
773 // We've already validated the static casts above.
774 const int numColsA = static_cast<int> (A.getNumVectors ());
775 const int numColsMv = static_cast<int> (mv.getNumVectors ());
776 if (numColsA > numColsMv) {
777 TEUCHOS_TEST_FOR_EXCEPTION(
778 numColsA > numColsMv, std::invalid_argument,
779 errPrefix << "Input multivector 'A' has " << numColsA << " columns, "
780 "but output multivector 'mv' has only " << numColsMv << " columns.");
781 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
782 }
783 if (numColsA == numColsMv) {
784 ::Tpetra::deep_copy (mv, A);
785 } else {
786 Teuchos::RCP<MV> mv_view =
787 CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
788 ::Tpetra::deep_copy (*mv_view, A);
789 }
790 }
791
792
793 static void MvRandom (MV& mv) {
794 mv.randomize ();
795 }
796
797 static void
798 MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
799 {
800 mv.putScalar (alpha);
801 }
802
803 static void MvPrint (const MV& mv, std::ostream& os) {
804 Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
805 mv.describe (fos, Teuchos::VERB_EXTREME);
806 }
807
808#ifdef HAVE_BELOS_TSQR
812 typedef Tpetra::TsqrAdaptor< Tpetra::MultiVector< Scalar, LO, GO, Node > > tsqr_adaptor_type;
813#endif // HAVE_BELOS_TSQR
814 };
815
817 //
818 // Implementation of the Belos::OperatorTraits for Tpetra::Operator.
819 //
821
823 template <class Storage, class LO, class GO, class Node>
824 class OperatorTraits <typename Storage::value_type,
825 Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
826 LO,GO,Node>,
827 Tpetra::Operator<Sacado::MP::Vector<Storage>,
828 LO,GO,Node> >
829 {
830 public:
832 static void
833 Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
834 const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
835 Tpetra::MultiVector<Scalar,LO,GO,Node>& Y,
836 ETrans trans=NOTRANS)
837 {
838 switch (trans) {
839 case NOTRANS:
840 Op.apply(X,Y,Teuchos::NO_TRANS);
841 break;
842 case TRANS:
843 Op.apply(X,Y,Teuchos::TRANS);
844 break;
845 case CONJTRANS:
846 Op.apply(X,Y,Teuchos::CONJ_TRANS);
847 break;
848 default:
849 const std::string scalarName = Teuchos::TypeNameTraits<Scalar>::name();
850 const std::string loName = Teuchos::TypeNameTraits<LO>::name();
851 const std::string goName = Teuchos::TypeNameTraits<GO>::name();
852 const std::string nodeName = Teuchos::TypeNameTraits<Node>::name();
853 const std::string otName = "Belos::OperatorTraits<" + scalarName
854 + "," + loName + "," + goName + "," + nodeName + ">";
855 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, otName << ": Should never "
856 "get here; fell through a switch statement. "
857 "Please report this bug to the Belos developers.");
858 }
859 }
860
861 static bool
862 HasApplyTranspose (const Tpetra::Operator<Scalar,LO,GO,Node>& Op)
863 {
864 return Op.hasTransposeApply ();
865 }
866 };
867
868} // end of Belos namespace
869
870#endif
871// end of file BELOS_TPETRA_ADAPTER_HPP
Kokkos::DefaultExecutionSpace execution_space
static void MvTimesMatAddMv(const dot_type &alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Teuchos::SerialDenseMatrix< int, dot_type > &B, const dot_type &beta, Tpetra::MultiVector< Scalar, LO, GO, Node > &C)
static void MvNorm(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, std::vector< mag_type > &normvec, NormType type=TwoNorm)
For all columns j of mv, set normvec[j] = norm(mv[j]).
static void MvAddMv(Scalar alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, Scalar beta, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
mv := alpha*A + beta*B
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const std::vector< int > &index)
Create and return a deep copy of the given columns of mv.
static Teuchos::RCP< MV > Clone(const MV &X, const int numVecs)
Create a new MultiVector with numVecs columns.
static void MvDot(const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, std::vector< dot_type > &dots)
For all columns j of A, set dots[j] := A[j]^T * B[j].
static void MvTransMv(dot_type alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, Teuchos::SerialDenseMatrix< int, dot_type > &C)
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< Scalar > &alphas)
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< BaseScalar > &alphas)
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const Teuchos::Range1D &index)
Create and return a deep copy of the given columns of mv.
static void Apply(const Tpetra::Operator< Scalar, LO, GO, Node > &Op, const Tpetra::MultiVector< Scalar, LO, GO, Node > &X, Tpetra::MultiVector< Scalar, LO, GO, Node > &Y, ETrans trans=NOTRANS)
std::enable_if< Kokkos::is_view_mp_vector< Kokkos::View< DA, PA... > >::value &&Kokkos::is_view_mp_vector< Kokkos::View< DB, PB... > >::value &&Kokkos::is_view_mp_vector< Kokkos::View< DC, PC... > >::value >::type gemm(const char transA[], const char transB[], typename Kokkos::View< DA, PA... >::const_value_type &alpha, const Kokkos::View< DA, PA... > &A, const Kokkos::View< DB, PB... > &B, typename Kokkos::View< DC, PC... >::const_value_type &beta, const Kokkos::View< DC, PC... > &C)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)