Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Belos_TpetraAdapter_UQ_PCE.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_UQ_PCE_HPP
43#define BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
44
45#include "BelosTpetraAdapter.hpp"
47#include "KokkosBlas.hpp"
48
49#ifdef HAVE_BELOS_TSQR
51#endif // HAVE_BELOS_TSQR
52
53namespace Belos {
54
56 //
57 // Implementation of Belos::MultiVecTraits for Tpetra::MultiVector.
58 //
60
71 template<class Storage, class LO, class GO, class Node>
72 class MultiVecTraits<typename Storage::value_type,
73 Tpetra::MultiVector< Sacado::UQ::PCE<Storage>,
74 LO, GO, Node > > {
75 public:
76 typedef typename Storage::ordinal_type s_ordinal;
77 typedef typename Storage::value_type BaseScalar;
79 typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
80 public:
81 typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::dot_type dot_type;
82 typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::mag_type mag_type;
83
89 static Teuchos::RCP<MV> Clone (const MV& X, const int numVecs) {
90 Teuchos::RCP<MV> Y (new MV (X.getMap (), numVecs, false));
91 Y->setCopyOrView (Teuchos::View);
92 return Y;
93 }
94
96 static Teuchos::RCP<MV> CloneCopy (const MV& X)
97 {
98 // Make a deep copy of X. The one-argument copy constructor
99 // does a shallow copy by default; the second argument tells it
100 // to do a deep copy.
101 Teuchos::RCP<MV> X_copy (new MV (X, Teuchos::Copy));
102 // Make Tpetra::MultiVector use the new view semantics. This is
103 // a no-op for the Kokkos refactor version of Tpetra; it only
104 // does something for the "classic" version of Tpetra. This
105 // shouldn't matter because Belos only handles MV through RCP
106 // and through this interface anyway, but it doesn't hurt to set
107 // it and make sure that it works.
108 X_copy->setCopyOrView (Teuchos::View);
109 return X_copy;
110 }
111
123 static Teuchos::RCP<MV>
124 CloneCopy (const MV& mv, const std::vector<int>& index)
125 {
126#ifdef HAVE_TPETRA_DEBUG
127 const char fnName[] = "Belos::MultiVecTraits::CloneCopy(mv,index)";
128 const size_t inNumVecs = mv.getNumVectors ();
129 TEUCHOS_TEST_FOR_EXCEPTION(
130 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
131 std::runtime_error, fnName << ": All indices must be nonnegative.");
132 TEUCHOS_TEST_FOR_EXCEPTION(
133 index.size () > 0 &&
134 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
135 std::runtime_error,
136 fnName << ": All indices must be strictly less than the number of "
137 "columns " << inNumVecs << " of the input multivector mv.");
138#endif // HAVE_TPETRA_DEBUG
139
140 // Tpetra wants an array of size_t, not of int.
141 Teuchos::Array<size_t> columns (index.size ());
142 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
143 columns[j] = index[j];
144 }
145 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
146 // continuous column index range in MultiVector::subCopy, so we
147 // don't have to check here.
148 Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
149 X_copy->setCopyOrView (Teuchos::View);
150 return X_copy;
151 }
152
159 static Teuchos::RCP<MV>
160 CloneCopy (const MV& mv, const Teuchos::Range1D& index)
161 {
162 const bool validRange = index.size() > 0 &&
163 index.lbound() >= 0 &&
164 index.ubound() < GetNumberVecs(mv);
165 if (! validRange) { // invalid range; generate error message
166 std::ostringstream os;
167 os << "Belos::MultiVecTraits::CloneCopy(mv,index=["
168 << index.lbound() << "," << index.ubound() << "]): ";
169 TEUCHOS_TEST_FOR_EXCEPTION(
170 index.size() == 0, std::invalid_argument,
171 os.str() << "Empty index range is not allowed.");
172 TEUCHOS_TEST_FOR_EXCEPTION(
173 index.lbound() < 0, std::invalid_argument,
174 os.str() << "Index range includes negative index/ices, which is not "
175 "allowed.");
176 TEUCHOS_TEST_FOR_EXCEPTION(
177 index.ubound() >= GetNumberVecs(mv), std::invalid_argument,
178 os.str() << "Index range exceeds number of vectors "
179 << mv.getNumVectors() << " in the input multivector.");
180 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
181 os.str() << "Should never get here!");
182 }
183 Teuchos::RCP<MV> X_copy = mv.subCopy (index);
184 X_copy->setCopyOrView (Teuchos::View);
185 return X_copy;
186 }
187
188
189 static Teuchos::RCP<MV>
190 CloneViewNonConst (MV& mv, const std::vector<int>& index)
191 {
192#ifdef HAVE_TPETRA_DEBUG
193 const char fnName[] = "Belos::MultiVecTraits::CloneViewNonConst(mv,index)";
194 const size_t numVecs = mv.getNumVectors ();
195 TEUCHOS_TEST_FOR_EXCEPTION(
196 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
197 std::invalid_argument,
198 fnName << ": All indices must be nonnegative.");
199 TEUCHOS_TEST_FOR_EXCEPTION(
200 index.size () > 0 &&
201 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
202 std::invalid_argument,
203 fnName << ": All indices must be strictly less than the number of "
204 "columns " << numVecs << " in the input MultiVector mv.");
205#endif // HAVE_TPETRA_DEBUG
206
207 // Tpetra wants an array of size_t, not of int.
208 Teuchos::Array<size_t> columns (index.size ());
209 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
210 columns[j] = index[j];
211 }
212 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
213 // continuous column index range in
214 // MultiVector::subViewNonConst, so we don't have to check here.
215 Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
216 X_view->setCopyOrView (Teuchos::View);
217 return X_view;
218 }
219
220 static Teuchos::RCP<MV>
221 CloneViewNonConst (MV& mv, const Teuchos::Range1D& index)
222 {
223 // NOTE (mfh 11 Jan 2011) We really should check for possible
224 // overflow of int here. However, the number of columns in a
225 // multivector typically fits in an int.
226 const int numCols = static_cast<int> (mv.getNumVectors());
227 const bool validRange = index.size() > 0 &&
228 index.lbound() >= 0 && index.ubound() < numCols;
229 if (! validRange) {
230 std::ostringstream os;
231 os << "Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
232 << index.lbound() << ", " << index.ubound() << "]): ";
233 TEUCHOS_TEST_FOR_EXCEPTION(
234 index.size() == 0, std::invalid_argument,
235 os.str() << "Empty index range is not allowed.");
236 TEUCHOS_TEST_FOR_EXCEPTION(
237 index.lbound() < 0, std::invalid_argument,
238 os.str() << "Index range includes negative inde{x,ices}, which is "
239 "not allowed.");
240 TEUCHOS_TEST_FOR_EXCEPTION(
241 index.ubound() >= numCols, std::invalid_argument,
242 os.str() << "Index range exceeds number of vectors " << numCols
243 << " in the input multivector.");
244 TEUCHOS_TEST_FOR_EXCEPTION(
245 true, std::logic_error,
246 os.str() << "Should never get here!");
247 }
248 Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
249 X_view->setCopyOrView (Teuchos::View);
250 return X_view;
251 }
252
253 static Teuchos::RCP<const MV>
254 CloneView (const MV& mv, const std::vector<int>& index)
255 {
256#ifdef HAVE_TPETRA_DEBUG
257 const char fnName[] = "Belos::MultiVecTraits<Scalar, "
258 "Tpetra::MultiVector<...> >::CloneView(mv,index)";
259 const size_t numVecs = mv.getNumVectors ();
260 TEUCHOS_TEST_FOR_EXCEPTION(
261 *std::min_element (index.begin (), index.end ()) < 0,
262 std::invalid_argument,
263 fnName << ": All indices must be nonnegative.");
264 TEUCHOS_TEST_FOR_EXCEPTION(
265 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
266 std::invalid_argument,
267 fnName << ": All indices must be strictly less than the number of "
268 "columns " << numVecs << " in the input MultiVector mv.");
269#endif // HAVE_TPETRA_DEBUG
270
271 // Tpetra wants an array of size_t, not of int.
272 Teuchos::Array<size_t> columns (index.size ());
273 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
274 columns[j] = index[j];
275 }
276 // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
277 // continuous column index range in MultiVector::subView, so we
278 // don't have to check here.
279 Teuchos::RCP<const MV> X_view = mv.subView (columns);
280 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
281 return X_view;
282 }
283
284 static Teuchos::RCP<const MV>
285 CloneView (const MV& mv, const Teuchos::Range1D& index)
286 {
287 // NOTE (mfh 11 Jan 2011) We really should check for possible
288 // overflow of int here. However, the number of columns in a
289 // multivector typically fits in an int.
290 const int numCols = static_cast<int> (mv.getNumVectors());
291 const bool validRange = index.size() > 0 &&
292 index.lbound() >= 0 && index.ubound() < numCols;
293 if (! validRange) {
294 std::ostringstream os;
295 os << "Belos::MultiVecTraits::CloneView(mv, index=["
296 << index.lbound () << ", " << index.ubound() << "]): ";
297 TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
298 os.str() << "Empty index range is not allowed.");
299 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
300 os.str() << "Index range includes negative index/ices, which is not "
301 "allowed.");
302 TEUCHOS_TEST_FOR_EXCEPTION(
303 index.ubound() >= numCols, std::invalid_argument,
304 os.str() << "Index range exceeds number of vectors " << numCols
305 << " in the input multivector.");
306 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
307 os.str() << "Should never get here!");
308 }
309 Teuchos::RCP<const MV> X_view = mv.subView (index);
310 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
311 return X_view;
312 }
313
314 static ptrdiff_t GetGlobalLength (const MV& mv) {
315 return static_cast<ptrdiff_t> (mv.getGlobalLength ());
316 }
317
318 static int GetNumberVecs (const MV& mv) {
319 return static_cast<int> (mv.getNumVectors ());
320 }
321
322 static bool HasConstantStride (const MV& mv) {
323 return mv.isConstantStride ();
324 }
325
326 static void
328 const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
329 const Teuchos::SerialDenseMatrix<int,dot_type>& B,
330 const dot_type& beta,
331 Tpetra::MultiVector<Scalar,LO,GO,Node>& C)
332 {
333 using Teuchos::RCP;
334 using Teuchos::rcp;
335
336 // Check if numRowsB == numColsB == 1, in which case we can call update()
337 const int numRowsB = B.numRows ();
338 const int numColsB = B.numCols ();
339 const int strideB = B.stride ();
340 if (numRowsB == 1 && numColsB == 1) {
341 C.update (alpha*B(0,0), A, beta);
342 return;
343 }
344
345 // Ensure A and C have constant stride
346 RCP<const MV> Atmp;
347 RCP<MV> Ctmp;
348 if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
349 else Atmp = rcp(&A,false);
350
351 if (C.isConstantStride() == false) Ctmp = rcp (new MV (C, Teuchos::Copy));
352 else Ctmp = rcp(&C,false);
353
354 // Create flattened view's
355 typedef Tpetra::MultiVector<dot_type,LO,GO,Node> FMV;
356 typedef Tpetra::MultiVector<const dot_type,LO,GO,Node> CFMV;
357 typedef typename FMV::dual_view_type::t_dev flat_view_type;
358 typedef typename CFMV::dual_view_type::t_dev const_flat_view_type;
359 typedef typename flat_view_type::execution_space execution_space;
360 const_flat_view_type flat_A_view = Atmp->template getLocalView<execution_space>(Tpetra::Access::ReadOnly);
361 flat_view_type flat_C_view = Ctmp->template getLocalView<execution_space>(Tpetra::Access::ReadWrite);
362
363 // Create a view for B on the host
364 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> b_host_view_type;
365 b_host_view_type B_view_host_input( B.values(), strideB, numColsB);
366 auto B_view_host = Kokkos::subview( B_view_host_input,
367 Kokkos::pair<int,int>(0, numRowsB),
368 Kokkos::pair<int,int>(0, numColsB));
369
370 // Create view for B on the device -- need to be careful to get the
371 // right stride to match B
372 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> b_view_type;
373 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> b_1d_view_type;
374 b_1d_view_type B_1d_view_dev(Kokkos::ViewAllocateWithoutInitializing("B"), numRowsB*numColsB);
375 b_view_type B_view_dev( B_1d_view_dev.data(), numRowsB, numColsB);
376 //Kokkos::deep_copy(B_view_dev, B_view_host);
377 // Device-to-host copies requires dest to be contiguous, which
378 // C_view_host may not be. So do 1 column at a time.
379 for (int j=0; j<numColsB; ++j)
380 Kokkos::deep_copy(Kokkos::subview(B_view_dev,Kokkos::ALL,j),
381 Kokkos::subview(B_view_host,Kokkos::ALL,j));
382
383 // Do local multiply
384 {
385 const char ctransA = 'N', ctransB = 'N';
387 &ctransA, &ctransB,
388 alpha, flat_A_view, B_view_dev, beta, flat_C_view);
389 }
390 // Copy back to C if we made a copy
391 if (C.isConstantStride() == false)
392 C.assign(*Ctmp);
393 }
394
402 static void
404 const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
405 Scalar beta,
406 const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
407 Tpetra::MultiVector<Scalar,LO,GO,Node>& mv)
408 {
409 mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
410 }
411
412 static void
413 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
414 const Scalar& alpha)
415 {
416 mv.scale (alpha);
417 }
418
419 static void
420 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
421 const std::vector<BaseScalar>& alphas)
422 {
423 std::vector<Scalar> alphas_mp(alphas.size());
424 const size_t sz = alphas.size();
425 for (size_t i=0; i<sz; ++i)
426 alphas_mp[i] = alphas[i];
427 mv.scale (alphas_mp);
428 }
429
430 static void
431 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
432 const std::vector<Scalar>& alphas)
433 {
434 mv.scale (alphas);
435 }
436
437 static void
439 const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
440 const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
441 Teuchos::SerialDenseMatrix<int,dot_type>& C)
442 {
443 using Teuchos::Comm;
444 using Teuchos::RCP;
445 using Teuchos::rcp;
446 using Teuchos::REDUCE_SUM;
447 using Teuchos::reduceAll;
448
449 // Check if numRowsC == numColsC == 1, in which case we can call dot()
450 const int numRowsC = C.numRows ();
451 const int numColsC = C.numCols ();
452 const int strideC = C.stride ();
453 if (numRowsC == 1 && numColsC == 1) {
454 if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
455 // Short-circuit, as required by BLAS semantics.
456 C(0,0) = alpha;
457 return;
458 }
459 A.dot (B, Teuchos::ArrayView<dot_type> (C.values (), 1));
460 if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
461 C(0,0) *= alpha;
462 }
463 return;
464 }
465
466 // Ensure A and B have constant stride
467 RCP<const MV> Atmp, Btmp;
468 if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
469 else Atmp = rcp(&A,false);
470
471 if (B.isConstantStride() == false) Btmp = rcp (new MV (B, Teuchos::Copy));
472 else Btmp = rcp(&B,false);
473
474 // Create flattened Kokkos::MultiVector's
475 typedef Tpetra::MultiVector<const dot_type,LO,GO,Node> FMV;
476 typedef typename FMV::dual_view_type::t_dev flat_view_type;
477 typedef typename flat_view_type::execution_space execution_space;
478 flat_view_type flat_A_view = Atmp->template getLocalView<execution_space>(Tpetra::Access::ReadOnly);
479 flat_view_type flat_B_view = Btmp->template getLocalView<execution_space>(Tpetra::Access::ReadOnly);
480
481 // Create a view for C on the host
482 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> c_host_view_type;
483 c_host_view_type C_view_host_input( C.values(), strideC, numColsC);
484 auto C_view_host = Kokkos::subview(C_view_host_input,
485 Kokkos::pair<int,int>(0, numRowsC),
486 Kokkos::pair<int,int>(0, numColsC));
487
488 // Create view for C on the device -- need to be careful to get the
489 // right stride to match C (allow setting to 0 for first-touch)
490 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> c_view_type;
491 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> c_1d_view_type;
492 c_1d_view_type C_1d_view_dev("C", numRowsC*numColsC);
493 c_view_type C_view_dev( C_1d_view_dev.data(), numRowsC, numColsC);
494
495 // Do local multiply
496 {
497 const char ctransA = 'C', ctransB = 'N';
499 &ctransA, &ctransB,
500 alpha, flat_A_view, flat_B_view,
501 Kokkos::Details::ArithTraits<dot_type>::zero(),
502 C_view_dev);
503 }
504 // reduce across processors -- could check for RDMA
505 RCP<const Comm<int> > pcomm = A.getMap()->getComm ();
506 if (pcomm->getSize () == 1) {
507 //Kokkos::deep_copy(C_view_host, C_view_dev);
508 // Device-to-host copies requires dest to be contiguous, which
509 // C_view_host may not be. So do 1 column at a time.
510 for (int j=0; j<numColsC; ++j)
511 Kokkos::deep_copy(Kokkos::subview(C_view_host,Kokkos::ALL,j),
512 Kokkos::subview(C_view_dev,Kokkos::ALL,j));
513 }
514 else {
515 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, Kokkos::HostSpace> c_1d_host_view_type;
516 c_1d_host_view_type C_1d_view_tmp(Kokkos::ViewAllocateWithoutInitializing("C_tmp"), strideC*numColsC);
517 c_host_view_type C_view_tmp( C_1d_view_tmp.data(),
518 strideC, numColsC);
519 for (int j=0; j<numColsC; ++j)
520 Kokkos::deep_copy(Kokkos::subview(C_view_tmp,
521 Kokkos::pair<int,int>(0, numRowsC),
522 j),
523 Kokkos::subview(C_view_dev, Kokkos::ALL, j));
524 reduceAll<int> (*pcomm, REDUCE_SUM, strideC*numColsC,
525 C_view_tmp.data(),
526 C_view_host.data());
527 }
528 }
529
531 static void
532 MvDot (const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
533 const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
534 std::vector<dot_type>& dots)
535 {
536 const size_t numVecs = A.getNumVectors ();
537
538 TEUCHOS_TEST_FOR_EXCEPTION(
539 numVecs != B.getNumVectors (), std::invalid_argument,
540 "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
541 "A and B must have the same number of columns. "
542 "A has " << numVecs << " column(s), "
543 "but B has " << B.getNumVectors () << " column(s).");
544#ifdef HAVE_TPETRA_DEBUG
545 TEUCHOS_TEST_FOR_EXCEPTION(
546 dots.size() < numVecs, std::invalid_argument,
547 "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
548 "The output array 'dots' must have room for all dot products. "
549 "A and B each have " << numVecs << " column(s), "
550 "but 'dots' only has " << dots.size() << " entry(/ies).");
551#endif // HAVE_TPETRA_DEBUG
552
553 Teuchos::ArrayView<dot_type> av (dots);
554 A.dot (B, av (0, numVecs));
555 }
556
558 static void
559 MvNorm (const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
560 std::vector<mag_type> &normvec,
561 NormType type=TwoNorm)
562 {
563
564#ifdef HAVE_TPETRA_DEBUG
565 typedef std::vector<int>::size_type size_type;
566 TEUCHOS_TEST_FOR_EXCEPTION(
567 normvec.size () < static_cast<size_type> (mv.getNumVectors ()),
568 std::invalid_argument,
569 "Belos::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
570 "argument must have at least as many entries as the number of vectors "
571 "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
572 << " < mv.getNumVectors() = " << mv.getNumVectors () << ".");
573#endif
574
575 Teuchos::ArrayView<mag_type> av(normvec);
576 switch (type) {
577 case OneNorm:
578 mv.norm1(av(0,mv.getNumVectors()));
579 break;
580 case TwoNorm:
581 mv.norm2(av(0,mv.getNumVectors()));
582 break;
583 case InfNorm:
584 mv.normInf(av(0,mv.getNumVectors()));
585 break;
586 default:
587 // Throw logic_error rather than invalid_argument, because if
588 // we get here, it's probably the fault of a Belos solver,
589 // rather than a user giving Belos an invalid input.
590 TEUCHOS_TEST_FOR_EXCEPTION(
591 true, std::logic_error,
592 "Belos::MultiVecTraits::MvNorm: Invalid NormType value " << type
593 << ". Valid values are OneNorm=" << OneNorm << ", TwoNorm="
594 << TwoNorm <<", and InfNorm=" << InfNorm << ". If you are a Belos "
595 "user and have not modified Belos in any way, and you get this "
596 "message, then this is probably a bug in the Belos solver you were "
597 "using. Please report this to the Belos developers.");
598 }
599 }
600
601 static void
602 SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
603 {
604 using Teuchos::Range1D;
605 using Teuchos::RCP;
606 const size_t inNumVecs = A.getNumVectors ();
607#ifdef HAVE_TPETRA_DEBUG
608 TEUCHOS_TEST_FOR_EXCEPTION(
609 inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
610 "Belos::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
611 "have no more entries as the number of columns in the input MultiVector"
612 " A. A.getNumVectors() = " << inNumVecs << " < index.size () = "
613 << index.size () << ".");
614#endif // HAVE_TPETRA_DEBUG
615 RCP<MV> mvsub = CloneViewNonConst (mv, index);
616 if (inNumVecs > static_cast<size_t> (index.size ())) {
617 RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
618 ::Tpetra::deep_copy (*mvsub, *Asub);
619 } else {
620 ::Tpetra::deep_copy (*mvsub, A);
621 }
622 }
623
624 static void
625 SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
626 {
627 // Range1D bounds are signed; size_t is unsigned.
628 // Assignment of Tpetra::MultiVector is a deep copy.
629
630 // Tpetra::MultiVector::getNumVectors() returns size_t. It's
631 // fair to assume that the number of vectors won't overflow int,
632 // since the typical use case of multivectors involves few
633 // columns, but it's friendly to check just in case.
634 const size_t maxInt =
635 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
636 const bool overflow =
637 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
638 if (overflow) {
639 std::ostringstream os;
640 os << "Belos::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
641 << ", " << index.ubound () << "], mv): ";
642 TEUCHOS_TEST_FOR_EXCEPTION(
643 maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
644 "of columns (size_t) in the input MultiVector 'A' overflows int.");
645 TEUCHOS_TEST_FOR_EXCEPTION(
646 maxInt < mv.getNumVectors (), std::range_error, os.str () << "Number "
647 "of columns (size_t) in the output MultiVector 'mv' overflows int.");
648 }
649 // We've already validated the static casts above.
650 const int numColsA = static_cast<int> (A.getNumVectors ());
651 const int numColsMv = static_cast<int> (mv.getNumVectors ());
652 // 'index' indexes into mv; it's the index set of the target.
653 const bool validIndex =
654 index.lbound () >= 0 && index.ubound () < numColsMv;
655 // We can't take more columns out of A than A has.
656 const bool validSource = index.size () <= numColsA;
657
658 if (! validIndex || ! validSource) {
659 std::ostringstream os;
660 os << "Belos::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
661 << ", " << index.ubound () << "], mv): ";
662 TEUCHOS_TEST_FOR_EXCEPTION(
663 index.lbound() < 0, std::invalid_argument,
664 os.str() << "Range lower bound must be nonnegative.");
665 TEUCHOS_TEST_FOR_EXCEPTION(
666 index.ubound() >= numColsMv, std::invalid_argument,
667 os.str() << "Range upper bound must be less than the number of "
668 "columns " << numColsA << " in the 'mv' output argument.");
669 TEUCHOS_TEST_FOR_EXCEPTION(
670 index.size() > numColsA, std::invalid_argument,
671 os.str() << "Range must have no more elements than the number of "
672 "columns " << numColsA << " in the 'A' input argument.");
673 TEUCHOS_TEST_FOR_EXCEPTION(
674 true, std::logic_error, "Should never get here!");
675 }
676
677 // View of the relevant column(s) of the target multivector mv.
678 // We avoid view creation overhead by only creating a view if
679 // the index range is different than [0, (# columns in mv) - 1].
680 Teuchos::RCP<MV> mv_view;
681 if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
682 mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
683 } else {
684 mv_view = CloneViewNonConst (mv, index);
685 }
686
687 // View of the relevant column(s) of the source multivector A.
688 // If A has fewer columns than mv_view, then create a view of
689 // the first index.size() columns of A.
690 Teuchos::RCP<const MV> A_view;
691 if (index.size () == numColsA) {
692 A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
693 } else {
694 A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
695 }
696
697 ::Tpetra::deep_copy (*mv_view, *A_view);
698 }
699
700 static void Assign (const MV& A, MV& mv)
701 {
702 const char errPrefix[] = "Belos::MultiVecTraits::Assign(A, mv): ";
703
704 // Range1D bounds are signed; size_t is unsigned.
705 // Assignment of Tpetra::MultiVector is a deep copy.
706
707 // Tpetra::MultiVector::getNumVectors() returns size_t. It's
708 // fair to assume that the number of vectors won't overflow int,
709 // since the typical use case of multivectors involves few
710 // columns, but it's friendly to check just in case.
711 const size_t maxInt =
712 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
713 const bool overflow =
714 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
715 if (overflow) {
716 TEUCHOS_TEST_FOR_EXCEPTION(
717 maxInt < A.getNumVectors(), std::range_error,
718 errPrefix << "Number of columns in the input multivector 'A' "
719 "(a size_t) overflows int.");
720 TEUCHOS_TEST_FOR_EXCEPTION(
721 maxInt < mv.getNumVectors(), std::range_error,
722 errPrefix << "Number of columns in the output multivector 'mv' "
723 "(a size_t) overflows int.");
724 TEUCHOS_TEST_FOR_EXCEPTION(
725 true, std::logic_error, "Should never get here!");
726 }
727 // We've already validated the static casts above.
728 const int numColsA = static_cast<int> (A.getNumVectors ());
729 const int numColsMv = static_cast<int> (mv.getNumVectors ());
730 if (numColsA > numColsMv) {
731 TEUCHOS_TEST_FOR_EXCEPTION(
732 numColsA > numColsMv, std::invalid_argument,
733 errPrefix << "Input multivector 'A' has " << numColsA << " columns, "
734 "but output multivector 'mv' has only " << numColsMv << " columns.");
735 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
736 }
737 if (numColsA == numColsMv) {
738 ::Tpetra::deep_copy (mv, A);
739 } else {
740 Teuchos::RCP<MV> mv_view =
741 CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
742 ::Tpetra::deep_copy (*mv_view, A);
743 }
744 }
745
746
747 static void MvRandom (MV& mv) {
748 mv.randomize ();
749 }
750
751 static void
752 MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
753 {
754 mv.putScalar (alpha);
755 }
756
757 static void MvPrint (const MV& mv, std::ostream& os) {
758 Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
759 mv.describe (fos, Teuchos::VERB_EXTREME);
760 }
761
762#ifdef HAVE_BELOS_TSQR
766 typedef Tpetra::TsqrAdaptor< Tpetra::MultiVector< Scalar, LO, GO, Node > > tsqr_adaptor_type;
767#endif // HAVE_BELOS_TSQR
768 };
769
771 //
772 // Implementation of the Belos::OperatorTraits for Tpetra::Operator.
773 //
775
777 template <class Storage, class LO, class GO, class Node>
778 class OperatorTraits <typename Storage::value_type,
779 Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
780 LO,GO,Node>,
781 Tpetra::Operator<Sacado::UQ::PCE<Storage>,
782 LO,GO,Node> >
783 {
784 public:
786 static void
787 Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
788 const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
789 Tpetra::MultiVector<Scalar,LO,GO,Node>& Y,
790 ETrans trans=NOTRANS)
791 {
792 switch (trans) {
793 case NOTRANS:
794 Op.apply(X,Y,Teuchos::NO_TRANS);
795 break;
796 case TRANS:
797 Op.apply(X,Y,Teuchos::TRANS);
798 break;
799 case CONJTRANS:
800 Op.apply(X,Y,Teuchos::CONJ_TRANS);
801 break;
802 default:
803 const std::string scalarName = Teuchos::TypeNameTraits<Scalar>::name();
804 const std::string loName = Teuchos::TypeNameTraits<LO>::name();
805 const std::string goName = Teuchos::TypeNameTraits<GO>::name();
806 const std::string nodeName = Teuchos::TypeNameTraits<Node>::name();
807 const std::string otName = "Belos::OperatorTraits<" + scalarName
808 + "," + loName + "," + goName + "," + nodeName + ">";
809 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, otName << ": Should never "
810 "get here; fell through a switch statement. "
811 "Please report this bug to the Belos developers.");
812 }
813 }
814
815 static bool
816 HasApplyTranspose (const Tpetra::Operator<Scalar,LO,GO,Node>& Op)
817 {
818 return Op.hasTransposeApply ();
819 }
820 };
821
822} // end of Belos namespace
823
824#endif
825// end of file BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
Kokkos::DefaultExecutionSpace execution_space
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 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 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 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 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 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 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 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 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)