Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_TpetraCrsMatrixUQPCEUnitTest.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Teuchos_UnitTestHelpers.hpp"
11 #include "Stokhos_UnitTestHelpers.hpp"
12 
13 // Teuchos
14 #include "Teuchos_XMLParameterListCoreHelpers.hpp"
15 
16 // Tpetra
17 #include "Stokhos_Tpetra_UQ_PCE.hpp"
18 #include "Stokhos_Tpetra_Utilities.hpp"
19 #include "Stokhos_Tpetra_Utilities_UQ_PCE.hpp"
20 #include "Tpetra_Core.hpp"
21 #include "Tpetra_Map.hpp"
22 #include "Tpetra_MultiVector.hpp"
23 #include "Tpetra_Vector.hpp"
24 #include "Tpetra_CrsGraph.hpp"
25 #include "Tpetra_CrsMatrix.hpp"
26 #include "Stokhos_Tpetra_CG.hpp"
27 
28 // Belos solver
29 #ifdef HAVE_STOKHOS_BELOS
30 #include "Belos_TpetraAdapter_UQ_PCE.hpp"
31 #include "BelosLinearProblem.hpp"
32 #include "BelosPseudoBlockGmresSolMgr.hpp"
33 #include "BelosPseudoBlockCGSolMgr.hpp"
34 #endif
35 
36 // Ifpack2 preconditioner
37 #ifdef HAVE_STOKHOS_IFPACK2
38 #include "Stokhos_Ifpack2_UQ_PCE.hpp"
39 #include "Ifpack2_Factory.hpp"
40 #endif
41 
42 // MueLu preconditioner
43 #ifdef HAVE_STOKHOS_MUELU
44 #include "Stokhos_MueLu_UQ_PCE.hpp"
45 #include "MueLu_CreateTpetraPreconditioner.hpp"
46 #endif
47 
48 // Amesos2 solver
49 #ifdef HAVE_STOKHOS_AMESOS2
50 #include "Stokhos_Amesos2_UQ_PCE.hpp"
51 #include "Amesos2_Factory.hpp"
52 #endif
53 
54 // Stokhos
55 #include "Stokhos_LegendreBasis.hpp"
56 #include "Stokhos_CompletePolynomialBasis.hpp"
57 #include "Stokhos_Sparse3Tensor.hpp"
58 
59 // Use "scalar" version of mean-based preconditioner (i.e., a preconditioner
60 // with double as the scalar type). This is currently necessary to get the
61 // MueLu tests to pass on OpenMP and Cuda due to various kernels that don't
62 // work with the PCE scalar type. Also required for RILUK test on Cuda
63 // without UVM (since factorization occurs on the host).
64 #define USE_SCALAR_MEAN_BASED_PREC 1
65 
66 template <typename scalar, typename ordinal>
67 inline
68 scalar generate_vector_coefficient( const ordinal nFEM,
69  const ordinal nStoch,
70  const ordinal iColFEM,
71  const ordinal iStoch )
72 {
73  const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
74  const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
75  return X_fem + X_stoch;
76  //return 1.0;
77 }
78 
79 template <typename scalar, typename ordinal>
80 inline
81 scalar generate_multi_vector_coefficient( const ordinal nFEM,
82  const ordinal nVec,
83  const ordinal nStoch,
84  const ordinal iColFEM,
85  const ordinal iVec,
86  const ordinal iStoch)
87 {
88  const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
89  const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
90  const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
91  return X_fem + X_stoch + X_col;
92  //return 1.0;
93 }
94 
95 template <typename scalar, typename ordinal>
96 inline
97 scalar generate_matrix_coefficient( const ordinal nFEM,
98  const ordinal nStoch,
99  const ordinal iRowFEM,
100  const ordinal iColFEM,
101  const ordinal iStoch )
102 {
103  const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
104  ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
105 
106  const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
107 
108  return A_fem + A_stoch;
109  //return 1.0;
110 }
111 
112 template <typename kokkos_cijk_type, typename ordinal_type>
113 kokkos_cijk_type build_cijk(ordinal_type stoch_dim,
114  ordinal_type poly_ord)
115 {
116  using Teuchos::RCP;
117  using Teuchos::rcp;
118  using Teuchos::Array;
119 
120  typedef typename kokkos_cijk_type::value_type value_type;
121  typedef typename kokkos_cijk_type::execution_space execution_space;
122  typedef Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> one_d_basis;
123  typedef Stokhos::LegendreBasis<ordinal_type,value_type> legendre_basis;
124  typedef Stokhos::CompletePolynomialBasis<ordinal_type,value_type> product_basis;
125  typedef Stokhos::Sparse3Tensor<ordinal_type,value_type> Cijk;
126 
127  // Create product basis
128  Array< RCP<const one_d_basis> > bases(stoch_dim);
129  for (ordinal_type i=0; i<stoch_dim; i++)
130  bases[i] = rcp(new legendre_basis(poly_ord, true));
131  RCP<const product_basis> basis = rcp(new product_basis(bases));
132 
133  // Triple product tensor
134  RCP<Cijk> cijk = basis->computeTripleProductTensor();
135 
136  // Kokkos triple product tensor
137  kokkos_cijk_type kokkos_cijk =
138  Stokhos::create_product_tensor<execution_space>(*basis, *cijk);
139 
140  return kokkos_cijk;
141 }
142 
143 //
144 // Tests
145 //
146 
147 // Stochastic discretizaiton used in the tests
148 const int stoch_dim = 2;
149 const int poly_ord = 3;
150 
151 //
152 // Test vector addition
153 //
154 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
155  Tpetra_CrsMatrix_PCE, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node )
156 {
157  using Teuchos::RCP;
158  using Teuchos::rcp;
159  using Teuchos::ArrayView;
160  using Teuchos::Array;
161  using Teuchos::ArrayRCP;
162 
163  typedef typename Storage::value_type BaseScalar;
164  typedef Sacado::UQ::PCE<Storage> Scalar;
165  typedef typename Scalar::cijk_type Cijk;
166 
167  typedef Teuchos::Comm<int> Tpetra_Comm;
168  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
169  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
170 
171  // Ensure device is initialized
172  if ( !Kokkos::is_initialized() )
173  Kokkos::initialize();
174 
175  // Cijk
176  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
177  Kokkos::setGlobalCijkTensor(cijk);
178  LocalOrdinal pce_size = cijk.dimension();
179 
180  // Comm
181  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
182 
183  // Map
184  GlobalOrdinal nrow = 10;
185  RCP<const Tpetra_Map> map =
186  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
187  nrow, comm);
188  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
189  const size_t num_my_row = myGIDs.size();
190 
191  // Fill vectors
192  RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
193  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
194  ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
195  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
196  Scalar val1(cijk), val2(cijk);
197  for (size_t i=0; i<num_my_row; ++i) {
198  const GlobalOrdinal row = myGIDs[i];
199  for (LocalOrdinal j=0; j<pce_size; ++j) {
200  val1.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row, j);
201  val2.fastAccessCoeff(j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row, j);
202  }
203  x1_view[i] = val1;
204  x2_view[i] = val2;
205  }
206  x1_view = Teuchos::null;
207  x2_view = Teuchos::null;
208 
209  // Add
210  Scalar alpha = 2.1;
211  Scalar beta = 3.7;
212  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
213  y->update(alpha, *x1, beta, *x2, Scalar(0.0));
214 
215  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
216  // Teuchos::VERB_EXTREME);
217 
218  // Check
219  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
220  Scalar val(cijk);
221  BaseScalar tol = 1.0e-14;
222  for (size_t i=0; i<num_my_row; ++i) {
223  const GlobalOrdinal row = myGIDs[i];
224  for (LocalOrdinal j=0; j<pce_size; ++j) {
225  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
226  nrow, pce_size, row, j);
227  val.fastAccessCoeff(j) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
228  }
229  TEST_EQUALITY( y_view[i].size(), pce_size );
230  for (LocalOrdinal j=0; j<pce_size; ++j)
231  TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(j), val.fastAccessCoeff(j), tol );
232  }
233 
234  // Clear global tensor
235  Kokkos::setGlobalCijkTensor(Cijk());
236 }
237 
238 //
239 // Test vector dot product
240 //
241 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
242  Tpetra_CrsMatrix_PCE, VectorDot, Storage, LocalOrdinal, GlobalOrdinal, Node )
243 {
244  using Teuchos::RCP;
245  using Teuchos::rcp;
246  using Teuchos::ArrayView;
247  using Teuchos::Array;
248  using Teuchos::ArrayRCP;
249 
250  typedef typename Storage::value_type BaseScalar;
251  typedef Sacado::UQ::PCE<Storage> Scalar;
252  typedef typename Scalar::cijk_type Cijk;
253 
254  typedef Teuchos::Comm<int> Tpetra_Comm;
255  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
256  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
257  typedef typename Tpetra_Vector::dot_type dot_type;
258 
259  // Ensure device is initialized
260  if ( !Kokkos::is_initialized() )
261  Kokkos::initialize();
262 
263  // Cijk
264  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
265  setGlobalCijkTensor(cijk);
266  LocalOrdinal pce_size = cijk.dimension();
267 
268  // Comm
269  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
270 
271  // Map
272  GlobalOrdinal nrow = 10;
273  RCP<const Tpetra_Map> map =
274  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
275  nrow, comm);
276  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
277  const size_t num_my_row = myGIDs.size();
278 
279  // Fill vectors
280  RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
281  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
282  ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
283  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
284  Scalar val1(cijk), val2(cijk);
285  for (size_t i=0; i<num_my_row; ++i) {
286  const GlobalOrdinal row = myGIDs[i];
287  for (LocalOrdinal j=0; j<pce_size; ++j) {
288  val1.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row, j);
289  val2.fastAccessCoeff(j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row, j);
290  }
291  x1_view[i] = val1;
292  x2_view[i] = val2;
293  }
294  x1_view = Teuchos::null;
295  x2_view = Teuchos::null;
296 
297  // Dot product
298  dot_type dot = x1->dot(*x2);
299 
300  // Check
301 
302  // Local contribution
303  dot_type local_val(0);
304  for (size_t i=0; i<num_my_row; ++i) {
305  const GlobalOrdinal row = myGIDs[i];
306  for (LocalOrdinal j=0; j<pce_size; ++j) {
307  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
308  nrow, pce_size, row, j);
309  local_val += 0.12345 * v * v;
310  }
311  }
312 
313  // Global reduction
314  dot_type val(0);
315  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
316  Teuchos::outArg(val));
317 
318  out << "dot = " << dot << " expected = " << val << std::endl;
319 
320  BaseScalar tol = 1.0e-14;
321  TEST_FLOATING_EQUALITY( dot, val, tol );
322 
323  // Clear global tensor
324  Kokkos::setGlobalCijkTensor(Cijk());
325 }
326 
327 //
328 // Test multi-vector addition
329 //
330 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
331  Tpetra_CrsMatrix_PCE, MultiVectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node )
332 {
333  using Teuchos::RCP;
334  using Teuchos::rcp;
335  using Teuchos::ArrayView;
336  using Teuchos::Array;
337  using Teuchos::ArrayRCP;
338 
339  typedef typename Storage::value_type BaseScalar;
340  typedef Sacado::UQ::PCE<Storage> Scalar;
341  typedef typename Scalar::cijk_type Cijk;
342 
343  typedef Teuchos::Comm<int> Tpetra_Comm;
344  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
345  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
346 
347  // Ensure device is initialized
348  if ( !Kokkos::is_initialized() )
349  Kokkos::initialize();
350 
351  // Cijk
352  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
353  setGlobalCijkTensor(cijk);
354  LocalOrdinal pce_size = cijk.dimension();
355 
356  // Comm
357  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
358 
359  // Map
360  GlobalOrdinal nrow = 10;
361  RCP<const Tpetra_Map> map =
362  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
363  nrow, comm);
364  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
365  const size_t num_my_row = myGIDs.size();
366 
367  // Fill vectors
368  size_t ncol = 5;
369  RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
370  RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
371  ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
372  ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
373  Scalar val1(cijk), val2(cijk);
374  for (size_t i=0; i<num_my_row; ++i) {
375  const GlobalOrdinal row = myGIDs[i];
376  for (size_t j=0; j<ncol; ++j) {
377  for (LocalOrdinal k=0; k<pce_size; ++k) {
378  BaseScalar v =
379  generate_multi_vector_coefficient<BaseScalar,size_t>(
380  nrow, ncol, pce_size, row, j, k);
381  val1.fastAccessCoeff(k) = v;
382  val2.fastAccessCoeff(k) = 0.12345 * v;
383  }
384  x1_view[j][i] = val1;
385  x2_view[j][i] = val2;
386  }
387  }
388  x1_view = Teuchos::null;
389  x2_view = Teuchos::null;
390 
391  // Add
392  Scalar alpha = 2.1;
393  Scalar beta = 3.7;
394  RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
395  y->update(alpha, *x1, beta, *x2, Scalar(0.0));
396 
397  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
398  // Teuchos::VERB_EXTREME);
399 
400  // Check
401  ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
402  Scalar val(cijk);
403  BaseScalar tol = 1.0e-14;
404  for (size_t i=0; i<num_my_row; ++i) {
405  const GlobalOrdinal row = myGIDs[i];
406  for (size_t j=0; j<ncol; ++j) {
407  for (LocalOrdinal k=0; k<pce_size; ++k) {
408  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
409  nrow, ncol, pce_size, row, j, k);
410  val.fastAccessCoeff(k) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
411  }
412  TEST_EQUALITY( y_view[j][i].size(), pce_size );
413  for (LocalOrdinal k=0; k<pce_size; ++k)
414  TEST_FLOATING_EQUALITY( y_view[j][i].fastAccessCoeff(k),
415  val.fastAccessCoeff(k), tol );
416  }
417  }
418 
419  // Clear global tensor
420  Kokkos::setGlobalCijkTensor(Cijk());
421 }
422 
423 //
424 // Test multi-vector dot product
425 //
426 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
427  Tpetra_CrsMatrix_PCE, MultiVectorDot, Storage, LocalOrdinal, GlobalOrdinal, Node )
428 {
429  using Teuchos::RCP;
430  using Teuchos::rcp;
431  using Teuchos::ArrayView;
432  using Teuchos::Array;
433  using Teuchos::ArrayRCP;
434 
435  typedef typename Storage::value_type BaseScalar;
436  typedef Sacado::UQ::PCE<Storage> Scalar;
437  typedef typename Scalar::cijk_type Cijk;
438 
439  typedef Teuchos::Comm<int> Tpetra_Comm;
440  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
441  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
442  typedef typename Tpetra_MultiVector::dot_type dot_type;
443 
444  // Ensure device is initialized
445  if ( !Kokkos::is_initialized() )
446  Kokkos::initialize();
447 
448  // Cijk
449  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
450  setGlobalCijkTensor(cijk);
451  LocalOrdinal pce_size = cijk.dimension();
452 
453  // Comm
454  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
455 
456  // Map
457  GlobalOrdinal nrow = 10;
458  RCP<const Tpetra_Map> map =
459  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
460  nrow, comm);
461  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
462  const size_t num_my_row = myGIDs.size();
463 
464  // Fill vectors
465  size_t ncol = 5;
466  RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
467  RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
468  ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
469  ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
470  Scalar val1(cijk), val2(cijk);
471  for (size_t i=0; i<num_my_row; ++i) {
472  const GlobalOrdinal row = myGIDs[i];
473  for (size_t j=0; j<ncol; ++j) {
474  for (LocalOrdinal k=0; k<pce_size; ++k) {
475  BaseScalar v =
476  generate_multi_vector_coefficient<BaseScalar,size_t>(
477  nrow, ncol, pce_size, row, j, k);
478  val1.fastAccessCoeff(k) = v;
479  val2.fastAccessCoeff(k) = 0.12345 * v;
480  }
481  x1_view[j][i] = val1;
482  x2_view[j][i] = val2;
483  }
484  }
485  x1_view = Teuchos::null;
486  x2_view = Teuchos::null;
487 
488  // Dot product
489  Array<dot_type> dots(ncol);
490  x1->dot(*x2, dots());
491 
492  // Check
493 
494  // Local contribution
495  Array<dot_type> local_vals(ncol, dot_type(0));
496  for (size_t i=0; i<num_my_row; ++i) {
497  const GlobalOrdinal row = myGIDs[i];
498  for (size_t j=0; j<ncol; ++j) {
499  for (LocalOrdinal k=0; k<pce_size; ++k) {
500  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
501  nrow, ncol, pce_size, row, j, k);
502  local_vals[j] += 0.12345 * v * v;
503  }
504  }
505  }
506 
507  // Global reduction
508  Array<dot_type> vals(ncol, dot_type(0));
509  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
510  local_vals.getRawPtr(), vals.getRawPtr());
511 
512  BaseScalar tol = 1.0e-14;
513  for (size_t j=0; j<ncol; ++j) {
514  out << "dots(" << j << ") = " << dots[j]
515  << " expected(" << j << ") = " << vals[j] << std::endl;
516  TEST_FLOATING_EQUALITY( dots[j], vals[j], tol );
517  }
518 
519  // Clear global tensor
520  Kokkos::setGlobalCijkTensor(Cijk());
521 }
522 
523 //
524 // Test multi-vector dot product using subviews
525 //
526 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
527  Tpetra_CrsMatrix_PCE, MultiVectorDotSub, Storage, LocalOrdinal, GlobalOrdinal, Node )
528 {
529  using Teuchos::RCP;
530  using Teuchos::rcp;
531  using Teuchos::ArrayView;
532  using Teuchos::Array;
533  using Teuchos::ArrayRCP;
534 
535  typedef typename Storage::value_type BaseScalar;
536  typedef Sacado::UQ::PCE<Storage> Scalar;
537  typedef typename Scalar::cijk_type Cijk;
538 
539  typedef Teuchos::Comm<int> Tpetra_Comm;
540  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
541  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
542  typedef typename Tpetra_MultiVector::dot_type dot_type;
543 
544  // Ensure device is initialized
545  if ( !Kokkos::is_initialized() )
546  Kokkos::initialize();
547 
548  // Cijk
549  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
550  setGlobalCijkTensor(cijk);
551  LocalOrdinal pce_size = cijk.dimension();
552 
553  // Comm
554  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
555 
556  // Map
557  GlobalOrdinal nrow = 10;
558  RCP<const Tpetra_Map> map =
559  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
560  nrow, comm);
561  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
562  const size_t num_my_row = myGIDs.size();
563 
564  // Fill vectors
565  size_t ncol = 5;
566  RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
567  RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
568  ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
569  ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
570  Scalar val1(cijk), val2(cijk);
571  for (size_t i=0; i<num_my_row; ++i) {
572  const GlobalOrdinal row = myGIDs[i];
573  for (size_t j=0; j<ncol; ++j) {
574  for (LocalOrdinal k=0; k<pce_size; ++k) {
575  BaseScalar v =
576  generate_multi_vector_coefficient<BaseScalar,size_t>(
577  nrow, ncol, pce_size, row, j, k);
578  val1.fastAccessCoeff(k) = v;
579  val2.fastAccessCoeff(k) = 0.12345 * v;
580  }
581  x1_view[j][i] = val1;
582  x2_view[j][i] = val2;
583  }
584  }
585  x1_view = Teuchos::null;
586  x2_view = Teuchos::null;
587 
588  // Get subviews
589  size_t ncol_sub = 2;
590  Teuchos::Array<size_t> cols(ncol_sub);
591  cols[0] = 4; cols[1] = 2;
592  RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
593  RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
594 
595  // Dot product
596  Array<dot_type> dots(ncol_sub);
597  x1_sub->dot(*x2_sub, dots());
598 
599  // Check
600 
601  // Local contribution
602  Array<dot_type> local_vals(ncol_sub, dot_type(0));
603  for (size_t i=0; i<num_my_row; ++i) {
604  const GlobalOrdinal row = myGIDs[i];
605  for (size_t j=0; j<ncol_sub; ++j) {
606  for (LocalOrdinal k=0; k<pce_size; ++k) {
607  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
608  nrow, ncol, pce_size, row, cols[j], k);
609  local_vals[j] += 0.12345 * v * v;
610  }
611  }
612  }
613 
614  // Global reduction
615  Array<dot_type> vals(ncol_sub, dot_type(0));
616  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
617  Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
618  vals.getRawPtr());
619 
620  BaseScalar tol = 1.0e-14;
621  for (size_t j=0; j<ncol_sub; ++j) {
622  out << "dots(" << j << ") = " << dots[j]
623  << " expected(" << j << ") = " << vals[j] << std::endl;
624  TEST_FLOATING_EQUALITY( dots[j], vals[j], tol );
625  }
626 
627  // Clear global tensor
628  Kokkos::setGlobalCijkTensor(Cijk());
629 }
630 
631 //
632 // Test matrix-vector multiplication for a simple banded upper-triangular matrix
633 //
634 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
635  Tpetra_CrsMatrix_PCE, MatrixVectorMultiply, Storage, LocalOrdinal, GlobalOrdinal, Node )
636 {
637  using Teuchos::RCP;
638  using Teuchos::rcp;
639  using Teuchos::ArrayView;
640  using Teuchos::Array;
641  using Teuchos::ArrayRCP;
642 
643  typedef typename Storage::value_type BaseScalar;
644  typedef Sacado::UQ::PCE<Storage> Scalar;
645  typedef typename Scalar::cijk_type Cijk;
646 
647  typedef Teuchos::Comm<int> Tpetra_Comm;
648  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
649  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
650  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
651  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
652 
653  // Ensure device is initialized
654  if ( !Kokkos::is_initialized() )
655  Kokkos::initialize();
656 
657  // Cijk
658  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
659  setGlobalCijkTensor(cijk);
660  LocalOrdinal pce_size = cijk.dimension();
661 
662  // Build banded matrix
663  GlobalOrdinal nrow = 13;
664  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
665  RCP<const Tpetra_Map> map =
666  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
667  nrow, comm);
668  RCP<Tpetra_CrsGraph> graph =
669  rcp(new Tpetra_CrsGraph(map, size_t(2)));
670  Array<GlobalOrdinal> columnIndices(2);
671  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
672  const size_t num_my_row = myGIDs.size();
673  for (size_t i=0; i<num_my_row; ++i) {
674  const GlobalOrdinal row = myGIDs[i];
675  columnIndices[0] = row;
676  size_t ncol = 1;
677  if (row != nrow-1) {
678  columnIndices[1] = row+1;
679  ncol = 2;
680  }
681  graph->insertGlobalIndices(row, columnIndices(0,ncol));
682  }
683  graph->fillComplete();
684  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
685 
686  // Set values in matrix
687  Array<Scalar> vals(2);
688  Scalar val(cijk);
689  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
690  const GlobalOrdinal row = myGIDs[local_row];
691  const size_t num_col = row == nrow - 1 ? 1 : 2;
692  for (size_t local_col=0; local_col<num_col; ++local_col) {
693  const GlobalOrdinal col = row + local_col;
694  columnIndices[local_col] = col;
695  for (LocalOrdinal k=0; k<pce_size; ++k)
696  val.fastAccessCoeff(k) =
697  generate_matrix_coefficient<BaseScalar,size_t>(
698  nrow, pce_size, row, col, k);
699  vals[local_col] = val;
700  }
701  matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
702  }
703  matrix->fillComplete();
704 
705  // Fill vector
706  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
707  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
708  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
709  const GlobalOrdinal row = myGIDs[local_row];
710  for (LocalOrdinal j=0; j<pce_size; ++j)
711  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
712  nrow, pce_size, row, j);
713  x_view[local_row] = val;
714  }
715  x_view = Teuchos::null;
716 
717  // matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
718  // Teuchos::VERB_EXTREME);
719 
720  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
721  // Teuchos::VERB_EXTREME);
722 
723  // Multiply
724  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
725  matrix->apply(*x, *y);
726 
727  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
728  // Teuchos::VERB_EXTREME);
729 
730  // Check
731  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
732  BaseScalar tol = 1.0e-14;
733  typename Cijk::host_mirror_type host_cijk =
734  Kokkos::create_mirror_view(cijk);
735  Kokkos::deep_copy(host_cijk, cijk);
736  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
737  const GlobalOrdinal row = myGIDs[local_row];
738  const size_t num_col = row == nrow - 1 ? 1 : 2;
739  val = 0.0;
740  for (size_t local_col=0; local_col<num_col; ++local_col) {
741  const GlobalOrdinal col = row + local_col;
742  for (LocalOrdinal i=0; i<pce_size; ++i) {
743  const LocalOrdinal num_entry = host_cijk.num_entry(i);
744  const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
745  const LocalOrdinal entry_end = entry_beg + num_entry;
746  BaseScalar tmp = 0;
747  for (LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
748  const LocalOrdinal j = host_cijk.coord(entry,0);
749  const LocalOrdinal k = host_cijk.coord(entry,1);
750  const BaseScalar a_j =
751  generate_matrix_coefficient<BaseScalar,size_t>(
752  nrow, pce_size, row, col, j);
753  const BaseScalar a_k =
754  generate_matrix_coefficient<BaseScalar,size_t>(
755  nrow, pce_size, row, col, k);
756  const BaseScalar x_j =
757  generate_vector_coefficient<BaseScalar,size_t>(
758  nrow, pce_size, col, j);
759  const BaseScalar x_k =
760  generate_vector_coefficient<BaseScalar,size_t>(
761  nrow, pce_size, col, k);
762  tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
763  }
764  val.fastAccessCoeff(i) += tmp;
765  }
766  }
767  TEST_EQUALITY( y_view[local_row].size(), pce_size );
768  for (LocalOrdinal i=0; i<pce_size; ++i)
769  TEST_FLOATING_EQUALITY( y_view[local_row].fastAccessCoeff(i),
770  val.fastAccessCoeff(i), tol );
771  }
772 
773  // Clear global tensor
774  Kokkos::setGlobalCijkTensor(Cijk());
775 }
776 
777 //
778 // Test matrix-multi-vector multiplication for a simple banded upper-triangular matrix
779 //
780 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
781  Tpetra_CrsMatrix_PCE, MatrixMultiVectorMultiply, Storage, LocalOrdinal, GlobalOrdinal, Node )
782 {
783  using Teuchos::RCP;
784  using Teuchos::rcp;
785  using Teuchos::ArrayView;
786  using Teuchos::Array;
787  using Teuchos::ArrayRCP;
788 
789  typedef typename Storage::value_type BaseScalar;
790  typedef Sacado::UQ::PCE<Storage> Scalar;
791  typedef typename Scalar::cijk_type Cijk;
792 
793  typedef Teuchos::Comm<int> Tpetra_Comm;
794  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
795  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
796  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
797  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
798 
799  // Ensure device is initialized
800  if ( !Kokkos::is_initialized() )
801  Kokkos::initialize();
802 
803  // Cijk
804  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
805  setGlobalCijkTensor(cijk);
806  LocalOrdinal pce_size = cijk.dimension();
807 
808  // Build banded matrix
809  GlobalOrdinal nrow = 10;
810  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
811  RCP<const Tpetra_Map> map =
812  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
813  nrow, comm);
814  RCP<Tpetra_CrsGraph> graph =
815  rcp(new Tpetra_CrsGraph(map, size_t(2)));
816  Array<GlobalOrdinal> columnIndices(2);
817  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
818  const size_t num_my_row = myGIDs.size();
819  for (size_t i=0; i<num_my_row; ++i) {
820  const GlobalOrdinal row = myGIDs[i];
821  columnIndices[0] = row;
822  size_t ncol = 1;
823  if (row != nrow-1) {
824  columnIndices[1] = row+1;
825  ncol = 2;
826  }
827  graph->insertGlobalIndices(row, columnIndices(0,ncol));
828  }
829  graph->fillComplete();
830  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
831 
832  // Set values in matrix
833  Array<Scalar> vals(2);
834  Scalar val(cijk);
835  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
836  const GlobalOrdinal row = myGIDs[local_row];
837  const size_t num_col = row == nrow - 1 ? 1 : 2;
838  for (size_t local_col=0; local_col<num_col; ++local_col) {
839  const GlobalOrdinal col = row + local_col;
840  columnIndices[local_col] = col;
841  for (LocalOrdinal k=0; k<pce_size; ++k)
842  val.fastAccessCoeff(k) =
843  generate_matrix_coefficient<BaseScalar,size_t>(
844  nrow, pce_size, row, col, k);
845  vals[local_col] = val;
846  }
847  matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
848  }
849  matrix->fillComplete();
850 
851  // Fill multi-vector
852  size_t ncol = 5;
853  RCP<Tpetra_MultiVector> x = Tpetra::createMultiVector<Scalar>(map, ncol);
854  ArrayRCP< ArrayRCP<Scalar> > x_view = x->get2dViewNonConst();
855  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
856  const GlobalOrdinal row = myGIDs[local_row];
857  for (size_t col=0; col<ncol; ++col) {
858  for (LocalOrdinal k=0; k<pce_size; ++k) {
859  BaseScalar v =
860  generate_multi_vector_coefficient<BaseScalar,size_t>(
861  nrow, ncol, pce_size, row, col, k);
862  val.fastAccessCoeff(k) = v;
863  }
864  x_view[col][local_row] = val;
865  }
866  }
867  x_view = Teuchos::null;
868 
869  // matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
870  // Teuchos::VERB_EXTREME);
871 
872  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
873  // Teuchos::VERB_EXTREME);
874 
875  // Multiply
876  RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
877  matrix->apply(*x, *y);
878 
879  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
880  // Teuchos::VERB_EXTREME);
881 
882  // Check
883  ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
884  BaseScalar tol = 1.0e-14;
885  typename Cijk::host_mirror_type host_cijk =
886  Kokkos::create_mirror_view(cijk);
887  Kokkos::deep_copy(host_cijk, cijk);
888  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
889  const GlobalOrdinal row = myGIDs[local_row];
890  for (size_t xcol=0; xcol<ncol; ++xcol) {
891  const size_t num_col = row == nrow - 1 ? 1 : 2;
892  val = 0.0;
893  for (size_t local_col=0; local_col<num_col; ++local_col) {
894  const GlobalOrdinal col = row + local_col;
895  for (LocalOrdinal i=0; i<pce_size; ++i) {
896  const LocalOrdinal num_entry = host_cijk.num_entry(i);
897  const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
898  const LocalOrdinal entry_end = entry_beg + num_entry;
899  BaseScalar tmp = 0;
900  for (LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
901  const LocalOrdinal j = host_cijk.coord(entry,0);
902  const LocalOrdinal k = host_cijk.coord(entry,1);
903  const BaseScalar a_j =
904  generate_matrix_coefficient<BaseScalar,size_t>(
905  nrow, pce_size, row, col, j);
906  const BaseScalar a_k =
907  generate_matrix_coefficient<BaseScalar,size_t>(
908  nrow, pce_size, row, col, k);
909  const BaseScalar x_j =
910  generate_multi_vector_coefficient<BaseScalar,size_t>(
911  nrow, ncol, pce_size, col, xcol, j);
912  const BaseScalar x_k =
913  generate_multi_vector_coefficient<BaseScalar,size_t>(
914  nrow, ncol, pce_size, col, xcol, k);
915  tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
916  }
917  val.fastAccessCoeff(i) += tmp;
918  }
919  }
920  TEST_EQUALITY( y_view[xcol][local_row].size(), pce_size );
921  for (LocalOrdinal i=0; i<pce_size; ++i)
922  TEST_FLOATING_EQUALITY( y_view[xcol][local_row].fastAccessCoeff(i),
923  val.fastAccessCoeff(i), tol );
924  }
925  }
926 
927  // Clear global tensor
928  Kokkos::setGlobalCijkTensor(Cijk());
929 }
930 
931 //
932 // Test flattening UQ::PCE matrix
933 //
934 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
935  Tpetra_CrsMatrix_PCE, Flatten, Storage, LocalOrdinal, GlobalOrdinal, Node )
936 {
937  using Teuchos::RCP;
938  using Teuchos::rcp;
939  using Teuchos::ArrayView;
940  using Teuchos::Array;
941  using Teuchos::ArrayRCP;
942 
943  typedef typename Storage::value_type BaseScalar;
944  typedef Sacado::UQ::PCE<Storage> Scalar;
945  typedef typename Scalar::cijk_type Cijk;
946 
947  typedef Teuchos::Comm<int> Tpetra_Comm;
948  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
949  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
950  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
951  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
952 
953  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
954  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
955 
956  // Ensure device is initialized
957  if ( !Kokkos::is_initialized() )
958  Kokkos::initialize();
959 
960  // Cijk
961  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
962  setGlobalCijkTensor(cijk);
963  LocalOrdinal pce_size = cijk.dimension();
964 
965  // Build banded matrix
966  GlobalOrdinal nrow = 10;
967  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
968  RCP<const Tpetra_Map> map =
969  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
970  nrow, comm);
971  RCP<Tpetra_CrsGraph> graph =
972  rcp(new Tpetra_CrsGraph(map, size_t(2)));
973  Array<GlobalOrdinal> columnIndices(2);
974  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
975  const size_t num_my_row = myGIDs.size();
976  for (size_t i=0; i<num_my_row; ++i) {
977  const GlobalOrdinal row = myGIDs[i];
978  columnIndices[0] = row;
979  size_t ncol = 1;
980  if (row != nrow-1) {
981  columnIndices[1] = row+1;
982  ncol = 2;
983  }
984  graph->insertGlobalIndices(row, columnIndices(0,ncol));
985  }
986  graph->fillComplete();
987  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
988 
989  // Set values in matrix
990  Array<Scalar> vals(2);
991  Scalar val(cijk);
992  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
993  const GlobalOrdinal row = myGIDs[local_row];
994  const size_t num_col = row == nrow - 1 ? 1 : 2;
995  for (size_t local_col=0; local_col<num_col; ++local_col) {
996  const GlobalOrdinal col = row + local_col;
997  columnIndices[local_col] = col;
998  for (LocalOrdinal k=0; k<pce_size; ++k)
999  val.fastAccessCoeff(k) =
1000  generate_matrix_coefficient<BaseScalar,size_t>(
1001  nrow, pce_size, row, col, k);
1002  vals[local_col] = val;
1003  }
1004  matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1005  }
1006  matrix->fillComplete();
1007 
1008  // Fill vector
1009  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1010  {
1011  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1012  for (size_t i=0; i<num_my_row; ++i) {
1013  const GlobalOrdinal row = myGIDs[i];
1014  for (LocalOrdinal j=0; j<pce_size; ++j)
1015  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
1016  nrow, pce_size, row, j);
1017  x_view[i] = val;
1018  }
1019  }
1020 
1021  // Multiply
1022  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
1023  matrix->apply(*x, *y);
1024 
1025  /*
1026  graph->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1027  Teuchos::VERB_EXTREME);
1028 
1029  matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1030  Teuchos::VERB_EXTREME);
1031 
1032  x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1033  Teuchos::VERB_EXTREME);
1034 
1035  y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1036  Teuchos::VERB_EXTREME);
1037  */
1038 
1039  // Flatten matrix
1040  RCP<const Tpetra_Map> flat_x_map, flat_y_map;
1041  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1042  flat_graph =
1043  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_y_map,
1044  cijk_graph, pce_size);
1045  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1046  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1047 
1048  // Multiply with flattened matix
1049  RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
1050  {
1051  RCP<Flat_Tpetra_Vector> flat_x =
1052  Stokhos::create_flat_vector_view(*x, flat_x_map);
1053  RCP<Flat_Tpetra_Vector> flat_y =
1054  Stokhos::create_flat_vector_view(*y2, flat_y_map);
1055  flat_matrix->apply(*flat_x, *flat_y);
1056 
1057  /*
1058  cijk_graph->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1059  Teuchos::VERB_EXTREME);
1060 
1061  flat_graph->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1062  Teuchos::VERB_EXTREME);
1063 
1064  flat_matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1065  Teuchos::VERB_EXTREME);
1066 
1067  flat_x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1068  Teuchos::VERB_EXTREME);
1069 
1070  flat_y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1071  Teuchos::VERB_EXTREME);
1072  */
1073  }
1074 
1075  // Check
1076  BaseScalar tol = 1.0e-14;
1077  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
1078  ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
1079  for (size_t i=0; i<num_my_row; ++i) {
1080  TEST_EQUALITY( y_view[i].size(), pce_size );
1081  TEST_EQUALITY( y2_view[i].size(), pce_size );
1082  for (LocalOrdinal j=0; j<pce_size; ++j)
1083  TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(j),
1084  y2_view[i].fastAccessCoeff(j), tol );
1085  }
1086 
1087  // Clear global tensor
1088  Kokkos::setGlobalCijkTensor(Cijk());
1089 }
1090 
1091 //
1092 // Test simple CG solve without preconditioning for a 1-D Laplacian matrix
1093 //
1094 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1095  Tpetra_CrsMatrix_PCE, SimpleCG, Storage, LocalOrdinal, GlobalOrdinal, Node )
1096 {
1097  using Teuchos::RCP;
1098  using Teuchos::rcp;
1099  using Teuchos::ArrayView;
1100  using Teuchos::Array;
1101  using Teuchos::ArrayRCP;
1102  using Teuchos::ParameterList;
1103 
1104  typedef typename Storage::value_type BaseScalar;
1105  typedef Sacado::UQ::PCE<Storage> Scalar;
1106  typedef typename Scalar::cijk_type Cijk;
1107 
1108  typedef Teuchos::Comm<int> Tpetra_Comm;
1109  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1110  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1111  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1112  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1113 
1114  // Ensure device is initialized
1115  if ( !Kokkos::is_initialized() )
1116  Kokkos::initialize();
1117 
1118  // Cijk
1119  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1120  setGlobalCijkTensor(cijk);
1121  LocalOrdinal pce_size = cijk.dimension();
1122 
1123  // 1-D Laplacian matrix
1124  GlobalOrdinal nrow = 10;
1125  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1126  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1127  RCP<const Tpetra_Map> map =
1128  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1129  nrow, comm);
1130  RCP<Tpetra_CrsGraph> graph =
1131  rcp(new Tpetra_CrsGraph(map, size_t(3)));
1132  Array<GlobalOrdinal> columnIndices(3);
1133  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1134  const size_t num_my_row = myGIDs.size();
1135  for (size_t i=0; i<num_my_row; ++i) {
1136  const GlobalOrdinal row = myGIDs[i];
1137  if (row == 0 || row == nrow-1) { // Boundary nodes
1138  columnIndices[0] = row;
1139  graph->insertGlobalIndices(row, columnIndices(0,1));
1140  }
1141  else { // Interior nodes
1142  columnIndices[0] = row-1;
1143  columnIndices[1] = row;
1144  columnIndices[2] = row+1;
1145  graph->insertGlobalIndices(row, columnIndices(0,3));
1146  }
1147  }
1148  graph->fillComplete();
1149  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1150 
1151  // Set values in matrix
1152  Array<Scalar> vals(3);
1153  Scalar a_val(cijk);
1154  for (LocalOrdinal j=0; j<pce_size; ++j) {
1155  a_val.fastAccessCoeff(j) =
1156  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(j+1);
1157  }
1158  for (size_t i=0; i<num_my_row; ++i) {
1159  const GlobalOrdinal row = myGIDs[i];
1160  if (row == 0 || row == nrow-1) { // Boundary nodes
1161  columnIndices[0] = row;
1162  vals[0] = Scalar(1.0);
1163  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1164  }
1165  else {
1166  columnIndices[0] = row-1;
1167  columnIndices[1] = row;
1168  columnIndices[2] = row+1;
1169  vals[0] = Scalar(1.0) * a_val;
1170  vals[1] = Scalar(-2.0) * a_val;
1171  vals[2] = Scalar(1.0) * a_val;
1172  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1173  }
1174  }
1175  matrix->fillComplete();
1176 
1177  // matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1178  // Teuchos::VERB_EXTREME);
1179 
1180  // Fill RHS vector
1181  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1182  {
1183  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1184  Scalar b_val(cijk);
1185  for (LocalOrdinal j=0; j<pce_size; ++j) {
1186  b_val.fastAccessCoeff(j) =
1187  BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(j+1);
1188  }
1189  for (size_t i=0; i<num_my_row; ++i) {
1190  const GlobalOrdinal row = myGIDs[i];
1191  if (row == 0 || row == nrow-1)
1192  b_view[i] = Scalar(0.0);
1193  else
1194  b_view[i] = b_val * (h*h);
1195  }
1196  }
1197 
1198  // b->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1199  // Teuchos::VERB_EXTREME);
1200 
1201  // Solve
1202  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1203 #if KOKKOS_VERSION >= 40799
1204  typedef KokkosKernels::ArithTraits<BaseScalar> BST;
1205 #else
1206  typedef Kokkos::ArithTraits<BaseScalar> BST;
1207 #endif
1208  typedef typename BST::mag_type base_mag_type;
1209  typedef typename Tpetra_Vector::mag_type mag_type;
1210  base_mag_type btol = 1e-9;
1211  mag_type tol = btol;
1212  int max_its = 1000;
1213  bool solved = Stokhos::CG_Solve(*matrix, *x, *b, tol, max_its,
1214  out.getOStream().get());
1215  TEST_EQUALITY_CONST( solved, true );
1216 
1217  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1218  // Teuchos::VERB_EXTREME);
1219 
1220  // Check by solving flattened system
1221  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1222  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1223  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1224  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1225  flat_graph =
1226  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1227  cijk_graph, pce_size);
1228  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1229  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1230  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1231  {
1232  RCP<Flat_Tpetra_Vector> flat_x =
1233  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1234  RCP<Flat_Tpetra_Vector> flat_b =
1235  Stokhos::create_flat_vector_view(*b, flat_b_map);
1236  bool solved_flat = Stokhos::CG_Solve(*flat_matrix, *flat_x, *flat_b,
1237  tol, max_its, out.getOStream().get());
1238  TEST_EQUALITY_CONST( solved_flat, true );
1239  }
1240 
1241  btol = 500*btol;
1242  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1243  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1244  for (size_t i=0; i<num_my_row; ++i) {
1245  TEST_EQUALITY( x_view[i].size(), pce_size );
1246  TEST_EQUALITY( x2_view[i].size(), pce_size );
1247 
1248  // Set small values to zero
1249  Scalar v = x_view[i];
1250  Scalar v2 = x2_view[i];
1251  for (LocalOrdinal j=0; j<pce_size; ++j) {
1252  if (j < v.size() && BST::abs(v.coeff(j)) < btol)
1253  v.fastAccessCoeff(j) = BaseScalar(0.0);
1254  if (j < v2.size() && BST::abs(v2.coeff(j)) < btol)
1255  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1256  }
1257 
1258  for (LocalOrdinal j=0; j<pce_size; ++j)
1259  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1260  }
1261 
1262  // Clear global tensor
1263  Kokkos::setGlobalCijkTensor(Cijk());
1264 
1265 }
1266 
1267 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1268 
1269 //
1270 // Test simple CG solve with MueLu preconditioning for a 1-D Laplacian matrix
1271 //
1272 // Currently requires ETI since the specializations needed for mean-based
1273 // are only brought in with ETI
1274 //
1275 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1276  Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
1277 {
1278  using Teuchos::RCP;
1279  using Teuchos::rcp;
1280  using Teuchos::ArrayView;
1281  using Teuchos::Array;
1282  using Teuchos::ArrayRCP;
1283  using Teuchos::ParameterList;
1284  using Teuchos::getParametersFromXmlFile;
1285 
1286  typedef typename Storage::value_type BaseScalar;
1287  typedef Sacado::UQ::PCE<Storage> Scalar;
1288  typedef typename Scalar::cijk_type Cijk;
1289 
1290  typedef Teuchos::Comm<int> Tpetra_Comm;
1291  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1292  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1293  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1294  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1295 
1296  // Ensure device is initialized
1297  if ( !Kokkos::is_initialized() )
1298  Kokkos::initialize();
1299 
1300  // Cijk
1301  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1302  setGlobalCijkTensor(cijk);
1303  LocalOrdinal pce_size = cijk.dimension();
1304 
1305  // 1-D Laplacian matrix
1306  GlobalOrdinal nrow = 10;
1307  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1308  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1309  RCP<const Tpetra_Map> map =
1310  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1311  nrow, comm);
1312  RCP<Tpetra_CrsGraph> graph =
1313  rcp(new Tpetra_CrsGraph(map, size_t(3)));
1314  Array<GlobalOrdinal> columnIndices(3);
1315  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1316  const size_t num_my_row = myGIDs.size();
1317  for (size_t i=0; i<num_my_row; ++i) {
1318  const GlobalOrdinal row = myGIDs[i];
1319  if (row == 0 || row == nrow-1) { // Boundary nodes
1320  columnIndices[0] = row;
1321  graph->insertGlobalIndices(row, columnIndices(0,1));
1322  }
1323  else { // Interior nodes
1324  columnIndices[0] = row-1;
1325  columnIndices[1] = row;
1326  columnIndices[2] = row+1;
1327  graph->insertGlobalIndices(row, columnIndices(0,3));
1328  }
1329  }
1330  graph->fillComplete();
1331  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1332 
1333  // Set values in matrix
1334  Array<Scalar> vals(3);
1335  Scalar a_val(cijk);
1336  for (LocalOrdinal j=0; j<pce_size; ++j) {
1337  a_val.fastAccessCoeff(j) =
1338  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(j+1);
1339  }
1340  for (size_t i=0; i<num_my_row; ++i) {
1341  const GlobalOrdinal row = myGIDs[i];
1342  if (row == 0 || row == nrow-1) { // Boundary nodes
1343  columnIndices[0] = row;
1344  vals[0] = Scalar(1.0);
1345  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1346  }
1347  else {
1348  columnIndices[0] = row-1;
1349  columnIndices[1] = row;
1350  columnIndices[2] = row+1;
1351  vals[0] = Scalar(1.0) * a_val;
1352  vals[1] = Scalar(-2.0) * a_val;
1353  vals[2] = Scalar(1.0) * a_val;
1354  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1355  }
1356  }
1357  matrix->fillComplete();
1358 
1359  // Fill RHS vector
1360  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1361  {
1362  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1363  Scalar b_val(cijk);
1364  for (LocalOrdinal j=0; j<pce_size; ++j) {
1365  b_val.fastAccessCoeff(j) =
1366  BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(j+1);
1367  }
1368  for (size_t i=0; i<num_my_row; ++i) {
1369  const GlobalOrdinal row = myGIDs[i];
1370  if (row == 0 || row == nrow-1)
1371  b_view[i] = Scalar(0.0);
1372  else
1373  b_view[i] = b_val * (h*h);
1374  }
1375  }
1376 
1377  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1378 #if KOKKOS_VERSION >= 40799
1379  typedef KokkosKernels::ArithTraits<BaseScalar> BST;
1380 #else
1381  typedef Kokkos::ArithTraits<BaseScalar> BST;
1382 #endif
1383  typedef typename BST::mag_type base_mag_type;
1384  typedef typename Tpetra_Vector::mag_type mag_type;
1385  base_mag_type btol = 1e-9;
1386  mag_type tol = btol;
1387  int max_its = 1000;
1388  {
1389  // Create preconditioner
1390  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1391  RCP<ParameterList> muelu_params =
1392  getParametersFromXmlFile("muelu_cheby.xml");
1393 #if USE_SCALAR_MEAN_BASED_PREC
1394  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1395  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1396  RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1397  Stokhos::build_mean_scalar_matrix(*matrix);
1398  RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1399  RCP<Scalar_OP> M_s =
1400  MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1401  RCP<OP> M = rcp(new Stokhos::MeanBasedTpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(M_s));
1402 #else
1403  Cijk mean_cijk =
1404  Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
1405  Kokkos::setGlobalCijkTensor(mean_cijk);
1406  RCP<Tpetra_CrsMatrix> mean_matrix = Stokhos::build_mean_matrix(*matrix);
1407  RCP<OP> mean_matrix_op = mean_matrix;
1408  RCP<OP> M =
1409  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1410  Kokkos::setGlobalCijkTensor(cijk);
1411 #endif
1412 
1413  // Solve
1414  bool solved = Stokhos::PCG_Solve(*matrix, *x, *b, *M, tol, max_its,
1415  out.getOStream().get());
1416  TEST_EQUALITY_CONST( solved, true );
1417  }
1418 
1419  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1420  // Teuchos::VERB_EXTREME);
1421 
1422  // Check by solving flattened system
1423  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1424  {
1425  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1426  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1427  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1428  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1429  flat_graph =
1430  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1431  cijk_graph, pce_size);
1432  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1433  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1434  RCP<Flat_Tpetra_Vector> flat_x =
1435  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1436  RCP<Flat_Tpetra_Vector> flat_b =
1437  Stokhos::create_flat_vector_view(*b, flat_b_map);
1438  // typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatPrec;
1439  // RCP<FlatPrec> flat_M =
1440  // MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(flat_matrix, *muelu_params);
1441  // bool solved_flat = Stokhos::PCG_Solve(*flat_matrix, *flat_x, *flat_b, *flat_M,
1442  // tol, max_its, out.getOStream().get());
1443  bool solved_flat = Stokhos::CG_Solve(*flat_matrix, *flat_x, *flat_b,
1444  tol, max_its, out.getOStream().get());
1445  TEST_EQUALITY_CONST( solved_flat, true );
1446  }
1447 
1448  btol = 500*btol;
1449  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1450  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1451  for (size_t i=0; i<num_my_row; ++i) {
1452  TEST_EQUALITY( x_view[i].size(), pce_size );
1453  TEST_EQUALITY( x2_view[i].size(), pce_size );
1454 
1455  // Set small values to zero
1456  Scalar v = x_view[i];
1457  Scalar v2 = x2_view[i];
1458  for (LocalOrdinal j=0; j<pce_size; ++j) {
1459  if (j < v.size() && BST::abs(v.coeff(j)) < btol)
1460  v.fastAccessCoeff(j) = BaseScalar(0.0);
1461  if (j < v2.size() && BST::abs(v2.coeff(j)) < btol)
1462  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1463  }
1464 
1465  for (LocalOrdinal j=0; j<pce_size; ++j)
1466  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1467  }
1468 
1469  // Clear global tensor
1470  Kokkos::setGlobalCijkTensor(Cijk());
1471 
1472 }
1473 
1474 #else
1475 
1476 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1477  Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node ) {}
1478 
1479 #endif
1480 
1481 #if defined(HAVE_STOKHOS_BELOS)
1482 
1483 //
1484 // Test Belos GMRES solve for a simple banded upper-triangular matrix
1485 //
1486 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1487  Tpetra_CrsMatrix_PCE, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1488 {
1489  using Teuchos::RCP;
1490  using Teuchos::rcp;
1491  using Teuchos::ArrayView;
1492  using Teuchos::Array;
1493  using Teuchos::ArrayRCP;
1494  using Teuchos::ParameterList;
1495 
1496  typedef typename Storage::value_type BaseScalar;
1497  typedef Sacado::UQ::PCE<Storage> Scalar;
1498  typedef typename Scalar::cijk_type Cijk;
1499 
1500  typedef Teuchos::Comm<int> Tpetra_Comm;
1501  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1502  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1503  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1504  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1505 
1506  // Ensure device is initialized
1507  if ( !Kokkos::is_initialized() )
1508  Kokkos::initialize();
1509 
1510  // Cijk
1511  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1512  setGlobalCijkTensor(cijk);
1513  LocalOrdinal pce_size = cijk.dimension();
1514 
1515  // Build banded matrix
1516  GlobalOrdinal nrow = 10;
1517  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1518  RCP<const Tpetra_Map> map =
1519  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1520  nrow, comm);
1521  RCP<Tpetra_CrsGraph> graph =
1522  rcp(new Tpetra_CrsGraph(map, size_t(2)));
1523  Array<GlobalOrdinal> columnIndices(2);
1524  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1525  const size_t num_my_row = myGIDs.size();
1526  for (size_t i=0; i<num_my_row; ++i) {
1527  const GlobalOrdinal row = myGIDs[i];
1528  columnIndices[0] = row;
1529  size_t ncol = 1;
1530  if (row != nrow-1) {
1531  columnIndices[1] = row+1;
1532  ncol = 2;
1533  }
1534  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1535  }
1536  graph->fillComplete();
1537  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1538 
1539  // Set values in matrix
1540  Array<Scalar> vals(2);
1541  Scalar val(cijk);
1542  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
1543  const GlobalOrdinal row = myGIDs[local_row];
1544  const size_t num_col = row == nrow - 1 ? 1 : 2;
1545  for (size_t local_col=0; local_col<num_col; ++local_col) {
1546  const GlobalOrdinal col = row + local_col;
1547  columnIndices[local_col] = col;
1548  for (LocalOrdinal k=0; k<pce_size; ++k)
1549  val.fastAccessCoeff(k) =
1550  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1551  vals[local_col] = val;
1552  }
1553  matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1554  }
1555  matrix->fillComplete();
1556 
1557  // Fill RHS vector
1558  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1559  {
1560  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1561  for (size_t i=0; i<num_my_row; ++i) {
1562  b_view[i] = Scalar(1.0);
1563  }
1564  }
1565 
1566  // Solve
1567  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1568  typedef BaseScalar BelosScalar;
1569  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1570  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1571  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1572  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1573  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1574  RCP<ParameterList> belosParams = rcp(new ParameterList);
1575  typename ST::magnitudeType tol = 1e-9;
1576  belosParams->set("Flexible Gmres", false);
1577  belosParams->set("Num Blocks", 100);
1578  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1579  belosParams->set("Maximum Iterations", 100);
1580  belosParams->set("Verbosity", 33);
1581  belosParams->set("Output Style", 1);
1582  belosParams->set("Output Frequency", 1);
1583  belosParams->set("Output Stream", out.getOStream());
1584  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1585  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1586  problem->setProblem();
1587  Belos::ReturnType ret = solver->solve();
1588  TEST_EQUALITY_CONST( ret, Belos::Converged );
1589 
1590  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1591  // Teuchos::VERB_EXTREME);
1592 
1593  // Check by solving flattened system
1594  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1595  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1596  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1597  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1598  flat_graph =
1599  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1600  cijk_graph, pce_size);
1601  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1602  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1603  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1604  {
1605  RCP<Flat_Tpetra_Vector> flat_x =
1606  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1607  RCP<Flat_Tpetra_Vector> flat_b =
1608  Stokhos::create_flat_vector_view(*b, flat_b_map);
1609  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1610  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1611  typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1612  RCP< FBLinProb > flat_problem =
1613  rcp(new FBLinProb(flat_matrix, flat_x, flat_b));
1614  RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1615  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1616  belosParams));
1617  flat_problem->setProblem();
1618  Belos::ReturnType flat_ret = flat_solver->solve();
1619  TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1620  }
1621 
1622  typename ST::magnitudeType btol = 100*tol;
1623  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1624  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1625  for (size_t i=0; i<num_my_row; ++i) {
1626  TEST_EQUALITY( x_view[i].size(), pce_size );
1627  TEST_EQUALITY( x2_view[i].size(), pce_size );
1628 
1629  // Set small values to zero
1630  Scalar v = x_view[i];
1631  Scalar v2 = x2_view[i];
1632  for (LocalOrdinal j=0; j<pce_size; ++j) {
1633  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
1634  v.fastAccessCoeff(j) = BaseScalar(0.0);
1635  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
1636  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1637  }
1638 
1639  for (LocalOrdinal j=0; j<pce_size; ++j)
1640  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1641  }
1642 
1643  // Clear global tensor
1644  Kokkos::setGlobalCijkTensor(Cijk());
1645 }
1646 
1647 #else
1648 
1649 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1650  Tpetra_CrsMatrix_PCE, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1651 {}
1652 
1653 #endif
1654 
1655 // Test currently doesn't work (in serial) because of our bad division strategy
1656 
1657 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
1658 
1659 //
1660 // Test Belos GMRES solve with Ifpack2 RILUK preconditioning for a
1661 // simple banded upper-triangular matrix
1662 //
1663 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1664  Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
1665 {
1666  using Teuchos::RCP;
1667  using Teuchos::rcp;
1668  using Teuchos::ArrayView;
1669  using Teuchos::Array;
1670  using Teuchos::ArrayRCP;
1671  using Teuchos::ParameterList;
1672 
1673  typedef typename Storage::value_type BaseScalar;
1674  typedef Sacado::UQ::PCE<Storage> Scalar;
1675  typedef typename Scalar::cijk_type Cijk;
1676 
1677  typedef Teuchos::Comm<int> Tpetra_Comm;
1678  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1679  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1680  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1681  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1682 
1683  // Ensure device is initialized
1684  if ( !Kokkos::is_initialized() )
1685  Kokkos::initialize();
1686 
1687  // Cijk
1688  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1689  setGlobalCijkTensor(cijk);
1690  LocalOrdinal pce_size = cijk.dimension();
1691 
1692  // Build banded matrix
1693  GlobalOrdinal nrow = 10;
1694  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1695  RCP<const Tpetra_Map> map =
1696  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1697  nrow, comm);
1698  RCP<Tpetra_CrsGraph> graph =
1699  rcp(new Tpetra_CrsGraph(map, size_t(2)));
1700  Array<GlobalOrdinal> columnIndices(2);
1701  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1702  const size_t num_my_row = myGIDs.size();
1703  for (size_t i=0; i<num_my_row; ++i) {
1704  const GlobalOrdinal row = myGIDs[i];
1705  columnIndices[0] = row;
1706  size_t ncol = 1;
1707  if (row != nrow-1) {
1708  columnIndices[1] = row+1;
1709  ncol = 2;
1710  }
1711  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1712  }
1713  graph->fillComplete();
1714  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1715 
1716  // Set values in matrix
1717  Array<Scalar> vals(2);
1718  Scalar val(cijk);
1719  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
1720  const GlobalOrdinal row = myGIDs[local_row];
1721  const size_t num_col = row == nrow - 1 ? 1 : 2;
1722  for (size_t local_col=0; local_col<num_col; ++local_col) {
1723  const GlobalOrdinal col = row + local_col;
1724  columnIndices[local_col] = col;
1725  for (LocalOrdinal k=0; k<pce_size; ++k)
1726  val.fastAccessCoeff(k) =
1727  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1728  vals[local_col] = val;
1729  }
1730  matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1731  }
1732  matrix->fillComplete();
1733 
1734  // Fill RHS vector
1735  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1736  {
1737  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1738  for (size_t i=0; i<num_my_row; ++i) {
1739  b_view[i] = Scalar(1.0);
1740  }
1741  }
1742 
1743  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1744  RCP<ParameterList> belosParams = rcp(new ParameterList);
1745  typedef BaseScalar BelosScalar;
1746  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1747  typename ST::magnitudeType tol = 1e-9;
1748  {
1749  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1750  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1751  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1752 
1753 #if USE_SCALAR_MEAN_BASED_PREC
1754  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1755  RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1756  Stokhos::build_mean_scalar_matrix(*matrix);
1757  typedef Ifpack2::Preconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Prec;
1758  Ifpack2::Factory factory;
1759  RCP<Scalar_Prec> M_s = factory.create<Scalar_Tpetra_CrsMatrix>("RILUK", mean_matrix);
1760  M_s->initialize();
1761  M_s->compute();
1762  RCP<OP> M = rcp(new Stokhos::MeanBasedTpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(M_s));
1763 #else
1764  RCP<Tpetra_CrsMatrix> mean_matrix = Stokhos::build_mean_matrix(*matrix);
1765  typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
1766  Ifpack2::Factory factory;
1767  RCP<Prec> M = factory.create<Tpetra_CrsMatrix>("RILUK", mean_matrix);
1768  M->initialize();
1769  M->compute();
1770 #endif
1771 
1772  // Solve
1773  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1774  problem->setRightPrec(M);
1775  problem->setProblem();
1776  belosParams->set("Flexible Gmres", false);
1777  belosParams->set("Num Blocks", 100);
1778  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1779  belosParams->set("Maximum Iterations", 100);
1780  belosParams->set("Verbosity", 33);
1781  belosParams->set("Output Style", 1);
1782  belosParams->set("Output Frequency", 1);
1783  belosParams->set("Output Stream", out.getOStream());
1784  //belosParams->set("Orthogonalization", "TSQR");
1785  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1786  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1787  Belos::ReturnType ret = solver->solve();
1788  TEST_EQUALITY_CONST( ret, Belos::Converged );
1789  }
1790 
1791  // Check by solving flattened system
1792  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1793  {
1794  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1795  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1796  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1797  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1798  flat_graph =
1799  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1800  cijk_graph, pce_size);
1801  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1802  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1803  RCP<Flat_Tpetra_Vector> flat_x =
1804  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1805  RCP<Flat_Tpetra_Vector> flat_b =
1806  Stokhos::create_flat_vector_view(*b, flat_b_map);
1807  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1808  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1809  typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1810  RCP< FBLinProb > flat_problem =
1811  rcp(new FBLinProb(flat_matrix, flat_x, flat_b));
1812  RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1813  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1814  belosParams));
1815  flat_problem->setProblem();
1816  Belos::ReturnType flat_ret = flat_solver->solve();
1817  TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1818  }
1819 
1820  typename ST::magnitudeType btol = 100*tol;
1821  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1822  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1823  for (size_t i=0; i<num_my_row; ++i) {
1824  TEST_EQUALITY( x_view[i].size(), pce_size );
1825  TEST_EQUALITY( x2_view[i].size(), pce_size );
1826 
1827  // Set small values to zero
1828  Scalar v = x_view[i];
1829  Scalar v2 = x2_view[i];
1830  for (LocalOrdinal j=0; j<pce_size; ++j) {
1831  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
1832  v.fastAccessCoeff(j) = BaseScalar(0.0);
1833  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
1834  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1835  }
1836 
1837  for (LocalOrdinal j=0; j<pce_size; ++j)
1838  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1839  }
1840 
1841  // Clear global tensor
1842  Kokkos::setGlobalCijkTensor(Cijk());
1843 }
1844 
1845 #else
1846 
1847 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1848  Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
1849 {}
1850 
1851 #endif
1852 
1853 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1854 
1855 //
1856 // Test Belos CG solve with MueLu preconditioning for a 1-D Laplacian matrix
1857 //
1858 // Currently requires ETI since the specializations needed for mean-based
1859 // are only brought in with ETI
1860 //
1861 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1862  Tpetra_CrsMatrix_PCE, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
1863 {
1864  using Teuchos::RCP;
1865  using Teuchos::rcp;
1866  using Teuchos::ArrayView;
1867  using Teuchos::Array;
1868  using Teuchos::ArrayRCP;
1869  using Teuchos::ParameterList;
1870  using Teuchos::getParametersFromXmlFile;
1871 
1872  typedef typename Storage::value_type BaseScalar;
1873  typedef Sacado::UQ::PCE<Storage> Scalar;
1874  typedef typename Scalar::cijk_type Cijk;
1875 
1876  typedef Teuchos::Comm<int> Tpetra_Comm;
1877  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1878  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1879  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1880  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1881 
1882  // Ensure device is initialized
1883  if ( !Kokkos::is_initialized() )
1884  Kokkos::initialize();
1885 
1886  // Cijk
1887  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1888  setGlobalCijkTensor(cijk);
1889  LocalOrdinal pce_size = cijk.dimension();
1890 
1891  // 1-D Laplacian matrix
1892  GlobalOrdinal nrow = 10;
1893  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1894  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1895  RCP<const Tpetra_Map> map =
1896  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1897  nrow, comm);
1898  RCP<Tpetra_CrsGraph> graph =
1899  rcp(new Tpetra_CrsGraph(map, size_t(3)));
1900  Array<GlobalOrdinal> columnIndices(3);
1901  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1902  const size_t num_my_row = myGIDs.size();
1903  for (size_t i=0; i<num_my_row; ++i) {
1904  const GlobalOrdinal row = myGIDs[i];
1905  if (row == 0 || row == nrow-1) { // Boundary nodes
1906  columnIndices[0] = row;
1907  graph->insertGlobalIndices(row, columnIndices(0,1));
1908  }
1909  else { // Interior nodes
1910  columnIndices[0] = row-1;
1911  columnIndices[1] = row;
1912  columnIndices[2] = row+1;
1913  graph->insertGlobalIndices(row, columnIndices(0,3));
1914  }
1915  }
1916  graph->fillComplete();
1917  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1918 
1919  // Set values in matrix
1920  Array<Scalar> vals(3);
1921  Scalar a_val(cijk);
1922  for (LocalOrdinal j=0; j<pce_size; ++j) {
1923  a_val.fastAccessCoeff(j) =
1924  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(j+1);
1925  }
1926  for (size_t i=0; i<num_my_row; ++i) {
1927  const GlobalOrdinal row = myGIDs[i];
1928  if (row == 0 || row == nrow-1) { // Boundary nodes
1929  columnIndices[0] = row;
1930  vals[0] = Scalar(1.0);
1931  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1932  }
1933  else {
1934  columnIndices[0] = row-1;
1935  columnIndices[1] = row;
1936  columnIndices[2] = row+1;
1937  vals[0] = Scalar(1.0) * a_val;
1938  vals[1] = Scalar(-2.0) * a_val;
1939  vals[2] = Scalar(1.0) * a_val;
1940  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1941  }
1942  }
1943  matrix->fillComplete();
1944 
1945  // Fill RHS vector
1946  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1947  {
1948  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1949  Scalar b_val(cijk);
1950  for (LocalOrdinal j=0; j<pce_size; ++j) {
1951  b_val.fastAccessCoeff(j) =
1952  BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(j+1);
1953  }
1954  for (size_t i=0; i<num_my_row; ++i) {
1955  const GlobalOrdinal row = myGIDs[i];
1956  if (row == 0 || row == nrow-1)
1957  b_view[i] = Scalar(0.0);
1958  else
1959  b_view[i] = b_val * (h*h);
1960  }
1961  }
1962 
1963  // Create preconditioner
1964  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1965  RCP<ParameterList> muelu_params =
1966  getParametersFromXmlFile("muelu_cheby.xml");
1967 #if USE_SCALAR_MEAN_BASED_PREC
1968  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1969  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1970  RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1971  Stokhos::build_mean_scalar_matrix(*matrix);
1972  RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1973  RCP<Scalar_OP> M_s =
1974  MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1975  RCP<OP> M = rcp(new Stokhos::MeanBasedTpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(M_s));
1976 #else
1977  Cijk mean_cijk =
1978  Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
1979  Kokkos::setGlobalCijkTensor(mean_cijk);
1980  RCP<Tpetra_CrsMatrix> mean_matrix = Stokhos::build_mean_matrix(*matrix);
1981  RCP<OP> mean_matrix_op = mean_matrix;
1982  RCP<OP> M =
1983  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1984  Kokkos::setGlobalCijkTensor(cijk);
1985 #endif
1986 
1987  // Solve
1988  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1989  typedef BaseScalar BelosScalar;
1990  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1991  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1992  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1993  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1994  problem->setRightPrec(M);
1995  problem->setProblem();
1996  RCP<ParameterList> belosParams = rcp(new ParameterList);
1997  typename ST::magnitudeType tol = 1e-9;
1998  belosParams->set("Num Blocks", 100);
1999  belosParams->set("Convergence Tolerance", BelosScalar(tol));
2000  belosParams->set("Maximum Iterations", 100);
2001  belosParams->set("Verbosity", 33);
2002  belosParams->set("Output Style", 1);
2003  belosParams->set("Output Frequency", 1);
2004  belosParams->set("Output Stream", out.getOStream());
2005  //belosParams->set("Orthogonalization", "TSQR");
2006  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2007  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2008  // RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2009  // rcp(new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2010  Belos::ReturnType ret = solver->solve();
2011  TEST_EQUALITY_CONST( ret, Belos::Converged );
2012 
2013  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2014  // Teuchos::VERB_EXTREME);
2015 
2016  // Check by solving flattened system
2017  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2018  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2019  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2020  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2021  flat_graph =
2022  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
2023  cijk_graph, pce_size);
2024  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2025  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
2026  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2027  {
2028  RCP<Flat_Tpetra_Vector> flat_x =
2029  Stokhos::create_flat_vector_view(*x2, flat_x_map);
2030  RCP<Flat_Tpetra_Vector> flat_b =
2031  Stokhos::create_flat_vector_view(*b, flat_b_map);
2032  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
2033  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
2034  typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
2035  RCP< FBLinProb > flat_problem =
2036  rcp(new FBLinProb(flat_matrix, flat_x, flat_b));
2037  RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
2038  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
2039  belosParams));
2040  // RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
2041  // rcp(new Belos::PseudoBlockCGSolMgr<BelosScalar,FMV,FOP>(flat_problem,
2042  // belosParams));
2043  flat_problem->setProblem();
2044  Belos::ReturnType flat_ret = flat_solver->solve();
2045  TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
2046  }
2047 
2048  typename ST::magnitudeType btol = 100*tol;
2049  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2050  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2051  for (size_t i=0; i<num_my_row; ++i) {
2052  TEST_EQUALITY( x_view[i].size(), pce_size );
2053  TEST_EQUALITY( x2_view[i].size(), pce_size );
2054 
2055  // Set small values to zero
2056  Scalar v = x_view[i];
2057  Scalar v2 = x2_view[i];
2058  for (LocalOrdinal j=0; j<pce_size; ++j) {
2059  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
2060  v.fastAccessCoeff(j) = BaseScalar(0.0);
2061  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
2062  v2.fastAccessCoeff(j) = BaseScalar(0.0);
2063  }
2064 
2065  for (LocalOrdinal j=0; j<pce_size; ++j)
2066  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
2067  }
2068 
2069  // Clear global tensor
2070  Kokkos::setGlobalCijkTensor(Cijk());
2071 
2072 }
2073 
2074 #else
2075 
2076 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
2077  Tpetra_CrsMatrix_PCE, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
2078 {}
2079 
2080 #endif
2081 
2082 #if defined(HAVE_STOKHOS_AMESOS2)
2083 
2084 //
2085 // Test Amesos2 solve for a 1-D Laplacian matrix
2086 //
2087 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
2088  Tpetra_CrsMatrix_PCE, Amesos2, Storage, LocalOrdinal, GlobalOrdinal, Node )
2089 {
2090  using Teuchos::RCP;
2091  using Teuchos::rcp;
2092  using Teuchos::ArrayView;
2093  using Teuchos::Array;
2094  using Teuchos::ArrayRCP;
2095  using Teuchos::ParameterList;
2096 
2097  typedef typename Storage::value_type BaseScalar;
2098  typedef Sacado::UQ::PCE<Storage> Scalar;
2099  typedef typename Scalar::cijk_type Cijk;
2100 
2101  typedef Teuchos::Comm<int> Tpetra_Comm;
2102  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2103  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2104  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
2105  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2106  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2107 
2108  // Ensure device is initialized
2109  if ( !Kokkos::is_initialized() )
2110  Kokkos::initialize();
2111 
2112  // Cijk
2113  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
2114  setGlobalCijkTensor(cijk);
2115  LocalOrdinal pce_size = cijk.dimension();
2116 
2117  // 1-D Laplacian matrix
2118  GlobalOrdinal nrow = 10;
2119  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
2120  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2121  RCP<const Tpetra_Map> map =
2122  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2123  nrow, comm);
2124  RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map, size_t(3));
2125  Array<GlobalOrdinal> columnIndices(3);
2126  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2127  const size_t num_my_row = myGIDs.size();
2128  for (size_t i=0; i<num_my_row; ++i) {
2129  const GlobalOrdinal row = myGIDs[i];
2130  if (row == 0 || row == nrow-1) { // Boundary nodes
2131  columnIndices[0] = row;
2132  graph->insertGlobalIndices(row, columnIndices(0,1));
2133  }
2134  else { // Interior nodes
2135  columnIndices[0] = row-1;
2136  columnIndices[1] = row;
2137  columnIndices[2] = row+1;
2138  graph->insertGlobalIndices(row, columnIndices(0,3));
2139  }
2140  }
2141  graph->fillComplete();
2142  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
2143 
2144  // Set values in matrix
2145  Array<Scalar> vals(3);
2146  Scalar a_val(cijk);
2147  for (LocalOrdinal j=0; j<pce_size; ++j) {
2148  a_val.fastAccessCoeff(j) =
2149  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(j+1);
2150  }
2151  for (size_t i=0; i<num_my_row; ++i) {
2152  const GlobalOrdinal row = myGIDs[i];
2153  if (row == 0 || row == nrow-1) { // Boundary nodes
2154  columnIndices[0] = row;
2155  vals[0] = Scalar(1.0);
2156  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2157  }
2158  else {
2159  columnIndices[0] = row-1;
2160  columnIndices[1] = row;
2161  columnIndices[2] = row+1;
2162  vals[0] = Scalar(1.0) * a_val;
2163  vals[1] = Scalar(-2.0) * a_val;
2164  vals[2] = Scalar(1.0) * a_val;
2165  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2166  }
2167  }
2168  matrix->fillComplete();
2169 
2170  // Fill RHS vector
2171  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2172  {
2173  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2174  Scalar b_val(cijk);
2175  for (LocalOrdinal j=0; j<pce_size; ++j) {
2176  b_val.fastAccessCoeff(j) =
2177  BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(j+1);
2178  }
2179  for (size_t i=0; i<num_my_row; ++i) {
2180  const GlobalOrdinal row = myGIDs[i];
2181  if (row == 0 || row == nrow-1)
2182  b_view[i] = Scalar(0.0);
2183  else
2184  b_view[i] = b_val * (h*h);
2185  }
2186  }
2187 
2188  // Solve
2189  typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2190  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2191  std::string solver_name;
2192 #if defined(HAVE_AMESOS2_BASKER)
2193  solver_name = "basker";
2194 #elif defined(HAVE_AMESOS2_KLU2)
2195  solver_name = "klu2";
2196 #elif defined(HAVE_AMESOS2_SUPERLUDIST)
2197  solver_name = "superlu_dist";
2198 #elif defined(HAVE_AMESOS2_SUPERLUMT)
2199  solver_name = "superlu_mt";
2200 #elif defined(HAVE_AMESOS2_SUPERLU)
2201  solver_name = "superlu";
2202 #elif defined(HAVE_AMESOS2_PARDISO_MKL)
2203  solver_name = "pardisomkl";
2204 #elif defined(HAVE_AMESOS2_LAPACK)
2205  solver_name = "lapack";
2206 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2207  solver_name = "lapack";
2208 #else
2209  // if there are no solvers, we just return as a successful test
2210  success = true;
2211  return;
2212 #endif
2213  out << "Solving linear system with " << solver_name << std::endl;
2214  {
2215  RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2216  solver_name, matrix, x, b);
2217  solver->solve();
2218  }
2219  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2220  // Teuchos::VERB_EXTREME);
2221 
2222  // Check by solving flattened system
2223  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2224  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_MultiVector;
2225  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2226  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2227  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2228  flat_graph =
2229  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
2230  cijk_graph, pce_size);
2231  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2232  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
2233  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2234  {
2235  RCP<Flat_Tpetra_Vector> flat_x =
2236  Stokhos::create_flat_vector_view(*x2, flat_x_map);
2237  RCP<Flat_Tpetra_Vector> flat_b =
2238  Stokhos::create_flat_vector_view(*b, flat_b_map);
2239  typedef Amesos2::Solver<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector> Flat_Solver;
2240  RCP<Flat_Solver> flat_solver =
2241  Amesos2::create<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector>(
2242  solver_name, flat_matrix, flat_x, flat_b);
2243  flat_solver->solve();
2244  }
2245 
2246 #if KOKKOS_VERSION >= 40799
2247  typedef KokkosKernels::ArithTraits<BaseScalar> ST;
2248 #else
2249  typedef Kokkos::ArithTraits<BaseScalar> ST;
2250 #endif
2251  typename ST::mag_type btol = 1e-12;
2252  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2253  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2254  for (size_t i=0; i<num_my_row; ++i) {
2255  TEST_EQUALITY( x_view[i].size(), pce_size );
2256  TEST_EQUALITY( x2_view[i].size(), pce_size );
2257 
2258  // Set small values to zero
2259  Scalar v = x_view[i];
2260  Scalar v2 = x2_view[i];
2261  for (LocalOrdinal j=0; j<pce_size; ++j) {
2262  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
2263  v.fastAccessCoeff(j) = BaseScalar(0.0);
2264  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
2265  v2.fastAccessCoeff(j) = BaseScalar(0.0);
2266  }
2267 
2268  for (LocalOrdinal j=0; j<pce_size; ++j)
2269  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
2270  }
2271 
2272  // Clear global tensor
2273  Kokkos::setGlobalCijkTensor(Cijk());
2274 }
2275 
2276 #else
2277 
2278 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
2279  Tpetra_CrsMatrix_PCE, Amesos2, Storage, LocalOrdinal, GlobalOrdinal, Node )
2280 {}
2281 
2282 #endif
2283 
2284 #define CRSMATRIX_UQ_PCE_TESTS_SLGN(S, LO, GO, N) \
2285  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorAdd, S, LO, GO, N ) \
2286  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorDot, S, LO, GO, N ) \
2287  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorAdd, S, LO, GO, N ) \
2288  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDot, S, LO, GO, N ) \
2289  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDotSub, S, LO, GO, N ) \
2290  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixVectorMultiply, S, LO, GO, N ) \
2291  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2292  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Flatten, S, LO, GO, N ) \
2293  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimpleCG, S, LO, GO, N ) \
2294  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, S, LO, GO, N ) \
2295  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES, S, LO, GO, N ) \
2296  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, S, LO, GO, N ) \
2297  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosCG_Muelu, S, LO, GO, N ) \
2298  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Amesos2, S, LO, GO, N )
2299 
2300 #define CRSMATRIX_UQ_PCE_TESTS_N(N) \
2301  typedef Stokhos::DeviceForNode2<N>::type Device; \
2302  typedef Stokhos::DynamicStorage<int,double,Device::execution_space> DS; \
2303  using default_local_ordinal_type = Tpetra::Map<>::local_ordinal_type; \
2304  using default_global_ordinal_type = Tpetra::Map<>::global_ordinal_type; \
2305  CRSMATRIX_UQ_PCE_TESTS_SLGN(DS, default_local_ordinal_type, default_global_ordinal_type, N)
Stokhos_Sparse3Tensor.hpp
Stokhos::create_flat_matrix
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &cijk_graph, const CijkType &cijk_dev)
Definition: Stokhos_Tpetra_Utilities_UQ_PCE.hpp:374
Stokhos::build_mean_matrix
Teuchos::RCP< Tpetra::CrsMatrix< Scalar, LO, GO, N > > build_mean_matrix(const Tpetra::CrsMatrix< Scalar, LO, GO, N > &A)
Definition: Stokhos_Tpetra_Utilities.hpp:207
Stokhos::CG_Solve
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
Definition: Stokhos_Tpetra_CG.hpp:19
Stokhos_Amesos2_UQ_PCE.hpp
Teuchos_UnitTestHelpers.hpp
generate_multi_vector_coefficient
scalar generate_multi_vector_coefficient(const ordinal nFEM, const ordinal nVec, const ordinal nStoch, const ordinal iColFEM, const ordinal iVec, const ordinal iStoch)
Definition: Stokhos_TpetraCrsMatrixMPVectorUnitTest.hpp:70
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
Definition: Stokhos_TpetraCrsMatrixMPVectorUnitTest.hpp:93
stoch_dim
const int stoch_dim
Definition: Stokhos_TpetraCrsMatrixUQPCEUnitTest.hpp:148
poly_ord
const int poly_ord
Definition: Stokhos_TpetraCrsMatrixUQPCEUnitTest.hpp:149
Stokhos::Sparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Definition: Stokhos_Sparse3Tensor.hpp:24
TEST_EQUALITY_CONST
#define TEST_EQUALITY_CONST(v1, v2)
DynamicVecTest::execution_space
Kokkos::DefaultExecutionSpace execution_space
Definition: Stokhos_SacadoMPVectorCommTests.cpp:651
Stokhos_Tpetra_UQ_PCE.hpp
Stokhos::create_flat_vector_view
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
Definition: Stokhos_Tpetra_Utilities_UQ_PCE.hpp:201
KokkosKernelsUnitTest::generate_vector_coefficient
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
Definition: Stokhos_KokkosArrayKernelsUnitTest.hpp:120
Teuchos_XMLParameterListCoreHelpers.hpp
KokkosBlas::dot
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< XD, XP...> >::value &&Kokkos::is_view_uq_pce< Kokkos::View< YD, YP...> >::value, typename Kokkos::Details::InnerProductSpaceTraits< typename Kokkos::View< XD, XP...>::non_const_value_type >::dot_type >::type dot(const Kokkos::View< XD, XP...> &x, const Kokkos::View< YD, YP...> &y)
Definition: Kokkos_Blas1_UQ_PCE.hpp:32
Kokkos::create_mirror_view
Stokhos::CrsMatrix< ValueType, Device, Layout >::host_mirror_type create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Definition: Stokhos_CrsMatrix.hpp:646
Stokhos::build_mean_scalar_matrix
Teuchos::RCP< Tpetra::CrsMatrix< typename Scalar::value_type, LO, GO, N > > build_mean_scalar_matrix(const Tpetra::CrsMatrix< Scalar, LO, GO, N > &A)
Definition: Stokhos_Tpetra_Utilities.hpp:234
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:495
Scalar
pce_type Scalar
Definition: linear2d_diffusion_scalar_types.hpp:17
Stokhos_MueLu_UQ_PCE.hpp
Stokhos_Tpetra_Utilities_UQ_PCE.hpp
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Belos_TpetraAdapter_UQ_PCE.hpp
Stokhos::PCG_Solve
bool PCG_Solve(const Matrix &A, Vector &x, const Vector &b, const Prec &M, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
Definition: Stokhos_Tpetra_CG.hpp:70
LocalOrdinal
int LocalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:18
Stokhos_Ifpack2_UQ_PCE.hpp
Kokkos::deep_copy
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Definition: Stokhos_CrsMatrix.hpp:657
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1150
fastAccessCoeff
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j)-expr2.val(j)
GlobalOrdinal
int GlobalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:19
Kokkos::setGlobalCijkTensor
void setGlobalCijkTensor(const cijk_type &cijk)
Definition: Kokkos_View_UQ_PCE.hpp:24
Stokhos::CompletePolynomialBasis
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Definition: Stokhos_CompletePolynomialBasis.hpp:40
build_cijk
kokkos_cijk_type build_cijk(ordinal_type stoch_dim, ordinal_type poly_ord)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:101
Stokhos_UnitTestHelpers.hpp
Stokhos_LegendreBasis.hpp
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Kokkos::cijk
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:88
Stokhos::LegendreBasis
Legendre polynomial basis.
Definition: Stokhos_LegendreBasis.hpp:33
KokkosKernelsUnitTest::generate_matrix_coefficient
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
Definition: Stokhos_KokkosArrayKernelsUnitTest.hpp:103
val
expr val()
Stokhos::OneDOrthogPolyBasis
Abstract base class for 1-D orthogonal polynomials.
Definition: Stokhos_OneDOrthogPolyBasis.hpp:49
Stokhos_Tpetra_CG.hpp
Stokhos_Tpetra_Utilities.hpp
Node
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: experimental/linear2d_diffusion_pce.cpp:47
Stokhos::create_flat_pce_graph
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_flat_pce_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, const CijkType &cijk, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_range_map, Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &cijk_graph, const size_t matrix_pce_size)
Definition: Stokhos_Tpetra_Utilities_UQ_PCE.hpp:89
Stokhos_CompletePolynomialBasis.hpp
BelosGMRES
BelosGMRES
Definition: tpetra_mat_vec.cpp:257
TEST_EQUALITY
#define TEST_EQUALITY(v1, v2)
Generated by   doxygen 1.8.5