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::HostMirror 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::HostMirror 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  typedef Kokkos::ArithTraits<BaseScalar> BST;
1204  typedef typename BST::mag_type base_mag_type;
1205  typedef typename Tpetra_Vector::mag_type mag_type;
1206  base_mag_type btol = 1e-9;
1207  mag_type tol = btol;
1208  int max_its = 1000;
1209  bool solved = Stokhos::CG_Solve(*matrix, *x, *b, tol, max_its,
1210  out.getOStream().get());
1211  TEST_EQUALITY_CONST( solved, true );
1212 
1213  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1214  // Teuchos::VERB_EXTREME);
1215 
1216  // Check by solving flattened system
1217  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1218  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1219  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1220  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1221  flat_graph =
1222  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1223  cijk_graph, pce_size);
1224  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1225  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1226  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1227  {
1228  RCP<Flat_Tpetra_Vector> flat_x =
1229  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1230  RCP<Flat_Tpetra_Vector> flat_b =
1231  Stokhos::create_flat_vector_view(*b, flat_b_map);
1232  bool solved_flat = Stokhos::CG_Solve(*flat_matrix, *flat_x, *flat_b,
1233  tol, max_its, out.getOStream().get());
1234  TEST_EQUALITY_CONST( solved_flat, true );
1235  }
1236 
1237  btol = 500*btol;
1238  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1239  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1240  for (size_t i=0; i<num_my_row; ++i) {
1241  TEST_EQUALITY( x_view[i].size(), pce_size );
1242  TEST_EQUALITY( x2_view[i].size(), pce_size );
1243 
1244  // Set small values to zero
1245  Scalar v = x_view[i];
1246  Scalar v2 = x2_view[i];
1247  for (LocalOrdinal j=0; j<pce_size; ++j) {
1248  if (j < v.size() && BST::abs(v.coeff(j)) < btol)
1249  v.fastAccessCoeff(j) = BaseScalar(0.0);
1250  if (j < v2.size() && BST::abs(v2.coeff(j)) < btol)
1251  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1252  }
1253 
1254  for (LocalOrdinal j=0; j<pce_size; ++j)
1255  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1256  }
1257 
1258  // Clear global tensor
1259  Kokkos::setGlobalCijkTensor(Cijk());
1260 
1261 }
1262 
1263 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1264 
1265 //
1266 // Test simple CG solve with MueLu preconditioning for a 1-D Laplacian matrix
1267 //
1268 // Currently requires ETI since the specializations needed for mean-based
1269 // are only brought in with ETI
1270 //
1271 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1272  Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
1273 {
1274  using Teuchos::RCP;
1275  using Teuchos::rcp;
1276  using Teuchos::ArrayView;
1277  using Teuchos::Array;
1278  using Teuchos::ArrayRCP;
1279  using Teuchos::ParameterList;
1280  using Teuchos::getParametersFromXmlFile;
1281 
1282  typedef typename Storage::value_type BaseScalar;
1283  typedef Sacado::UQ::PCE<Storage> Scalar;
1284  typedef typename Scalar::cijk_type Cijk;
1285 
1286  typedef Teuchos::Comm<int> Tpetra_Comm;
1287  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1288  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1289  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1290  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1291 
1292  // Ensure device is initialized
1293  if ( !Kokkos::is_initialized() )
1294  Kokkos::initialize();
1295 
1296  // Cijk
1297  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1298  setGlobalCijkTensor(cijk);
1299  LocalOrdinal pce_size = cijk.dimension();
1300 
1301  // 1-D Laplacian matrix
1302  GlobalOrdinal nrow = 10;
1303  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1304  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1305  RCP<const Tpetra_Map> map =
1306  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1307  nrow, comm);
1308  RCP<Tpetra_CrsGraph> graph =
1309  rcp(new Tpetra_CrsGraph(map, size_t(3)));
1310  Array<GlobalOrdinal> columnIndices(3);
1311  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1312  const size_t num_my_row = myGIDs.size();
1313  for (size_t i=0; i<num_my_row; ++i) {
1314  const GlobalOrdinal row = myGIDs[i];
1315  if (row == 0 || row == nrow-1) { // Boundary nodes
1316  columnIndices[0] = row;
1317  graph->insertGlobalIndices(row, columnIndices(0,1));
1318  }
1319  else { // Interior nodes
1320  columnIndices[0] = row-1;
1321  columnIndices[1] = row;
1322  columnIndices[2] = row+1;
1323  graph->insertGlobalIndices(row, columnIndices(0,3));
1324  }
1325  }
1326  graph->fillComplete();
1327  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1328 
1329  // Set values in matrix
1330  Array<Scalar> vals(3);
1331  Scalar a_val(cijk);
1332  for (LocalOrdinal j=0; j<pce_size; ++j) {
1333  a_val.fastAccessCoeff(j) =
1334  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(j+1);
1335  }
1336  for (size_t i=0; i<num_my_row; ++i) {
1337  const GlobalOrdinal row = myGIDs[i];
1338  if (row == 0 || row == nrow-1) { // Boundary nodes
1339  columnIndices[0] = row;
1340  vals[0] = Scalar(1.0);
1341  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1342  }
1343  else {
1344  columnIndices[0] = row-1;
1345  columnIndices[1] = row;
1346  columnIndices[2] = row+1;
1347  vals[0] = Scalar(1.0) * a_val;
1348  vals[1] = Scalar(-2.0) * a_val;
1349  vals[2] = Scalar(1.0) * a_val;
1350  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1351  }
1352  }
1353  matrix->fillComplete();
1354 
1355  // Fill RHS vector
1356  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1357  {
1358  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1359  Scalar b_val(cijk);
1360  for (LocalOrdinal j=0; j<pce_size; ++j) {
1361  b_val.fastAccessCoeff(j) =
1362  BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(j+1);
1363  }
1364  for (size_t i=0; i<num_my_row; ++i) {
1365  const GlobalOrdinal row = myGIDs[i];
1366  if (row == 0 || row == nrow-1)
1367  b_view[i] = Scalar(0.0);
1368  else
1369  b_view[i] = b_val * (h*h);
1370  }
1371  }
1372 
1373  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1374  typedef Kokkos::ArithTraits<BaseScalar> BST;
1375  typedef typename BST::mag_type base_mag_type;
1376  typedef typename Tpetra_Vector::mag_type mag_type;
1377  base_mag_type btol = 1e-9;
1378  mag_type tol = btol;
1379  int max_its = 1000;
1380  {
1381  // Create preconditioner
1382  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1383  RCP<ParameterList> muelu_params =
1384  getParametersFromXmlFile("muelu_cheby.xml");
1385 #if USE_SCALAR_MEAN_BASED_PREC
1386  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1387  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1388  RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1389  Stokhos::build_mean_scalar_matrix(*matrix);
1390  RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1391  RCP<Scalar_OP> M_s =
1392  MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1393  RCP<OP> M = rcp(new Stokhos::MeanBasedTpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(M_s));
1394 #else
1395  Cijk mean_cijk =
1396  Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
1397  Kokkos::setGlobalCijkTensor(mean_cijk);
1398  RCP<Tpetra_CrsMatrix> mean_matrix = Stokhos::build_mean_matrix(*matrix);
1399  RCP<OP> mean_matrix_op = mean_matrix;
1400  RCP<OP> M =
1401  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1402  Kokkos::setGlobalCijkTensor(cijk);
1403 #endif
1404 
1405  // Solve
1406  bool solved = Stokhos::PCG_Solve(*matrix, *x, *b, *M, tol, max_its,
1407  out.getOStream().get());
1408  TEST_EQUALITY_CONST( solved, true );
1409  }
1410 
1411  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1412  // Teuchos::VERB_EXTREME);
1413 
1414  // Check by solving flattened system
1415  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1416  {
1417  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1418  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1419  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1420  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1421  flat_graph =
1422  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1423  cijk_graph, pce_size);
1424  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1425  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1426  RCP<Flat_Tpetra_Vector> flat_x =
1427  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1428  RCP<Flat_Tpetra_Vector> flat_b =
1429  Stokhos::create_flat_vector_view(*b, flat_b_map);
1430  // typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatPrec;
1431  // RCP<FlatPrec> flat_M =
1432  // MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(flat_matrix, *muelu_params);
1433  // bool solved_flat = Stokhos::PCG_Solve(*flat_matrix, *flat_x, *flat_b, *flat_M,
1434  // tol, max_its, out.getOStream().get());
1435  bool solved_flat = Stokhos::CG_Solve(*flat_matrix, *flat_x, *flat_b,
1436  tol, max_its, out.getOStream().get());
1437  TEST_EQUALITY_CONST( solved_flat, true );
1438  }
1439 
1440  btol = 500*btol;
1441  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1442  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1443  for (size_t i=0; i<num_my_row; ++i) {
1444  TEST_EQUALITY( x_view[i].size(), pce_size );
1445  TEST_EQUALITY( x2_view[i].size(), pce_size );
1446 
1447  // Set small values to zero
1448  Scalar v = x_view[i];
1449  Scalar v2 = x2_view[i];
1450  for (LocalOrdinal j=0; j<pce_size; ++j) {
1451  if (j < v.size() && BST::abs(v.coeff(j)) < btol)
1452  v.fastAccessCoeff(j) = BaseScalar(0.0);
1453  if (j < v2.size() && BST::abs(v2.coeff(j)) < btol)
1454  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1455  }
1456 
1457  for (LocalOrdinal j=0; j<pce_size; ++j)
1458  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1459  }
1460 
1461  // Clear global tensor
1462  Kokkos::setGlobalCijkTensor(Cijk());
1463 
1464 }
1465 
1466 #else
1467 
1468 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1469  Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node ) {}
1470 
1471 #endif
1472 
1473 #if defined(HAVE_STOKHOS_BELOS)
1474 
1475 //
1476 // Test Belos GMRES solve for a simple banded upper-triangular matrix
1477 //
1478 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1479  Tpetra_CrsMatrix_PCE, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1480 {
1481  using Teuchos::RCP;
1482  using Teuchos::rcp;
1483  using Teuchos::ArrayView;
1484  using Teuchos::Array;
1485  using Teuchos::ArrayRCP;
1486  using Teuchos::ParameterList;
1487 
1488  typedef typename Storage::value_type BaseScalar;
1489  typedef Sacado::UQ::PCE<Storage> Scalar;
1490  typedef typename Scalar::cijk_type Cijk;
1491 
1492  typedef Teuchos::Comm<int> Tpetra_Comm;
1493  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1494  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1495  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1496  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1497 
1498  // Ensure device is initialized
1499  if ( !Kokkos::is_initialized() )
1500  Kokkos::initialize();
1501 
1502  // Cijk
1503  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1504  setGlobalCijkTensor(cijk);
1505  LocalOrdinal pce_size = cijk.dimension();
1506 
1507  // Build banded matrix
1508  GlobalOrdinal nrow = 10;
1509  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1510  RCP<const Tpetra_Map> map =
1511  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1512  nrow, comm);
1513  RCP<Tpetra_CrsGraph> graph =
1514  rcp(new Tpetra_CrsGraph(map, size_t(2)));
1515  Array<GlobalOrdinal> columnIndices(2);
1516  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1517  const size_t num_my_row = myGIDs.size();
1518  for (size_t i=0; i<num_my_row; ++i) {
1519  const GlobalOrdinal row = myGIDs[i];
1520  columnIndices[0] = row;
1521  size_t ncol = 1;
1522  if (row != nrow-1) {
1523  columnIndices[1] = row+1;
1524  ncol = 2;
1525  }
1526  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1527  }
1528  graph->fillComplete();
1529  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1530 
1531  // Set values in matrix
1532  Array<Scalar> vals(2);
1533  Scalar val(cijk);
1534  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
1535  const GlobalOrdinal row = myGIDs[local_row];
1536  const size_t num_col = row == nrow - 1 ? 1 : 2;
1537  for (size_t local_col=0; local_col<num_col; ++local_col) {
1538  const GlobalOrdinal col = row + local_col;
1539  columnIndices[local_col] = col;
1540  for (LocalOrdinal k=0; k<pce_size; ++k)
1541  val.fastAccessCoeff(k) =
1542  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1543  vals[local_col] = val;
1544  }
1545  matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1546  }
1547  matrix->fillComplete();
1548 
1549  // Fill RHS vector
1550  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1551  {
1552  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1553  for (size_t i=0; i<num_my_row; ++i) {
1554  b_view[i] = Scalar(1.0);
1555  }
1556  }
1557 
1558  // Solve
1559  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1560  typedef BaseScalar BelosScalar;
1561  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1562  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1563  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1564  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1565  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1566  RCP<ParameterList> belosParams = rcp(new ParameterList);
1567  typename ST::magnitudeType tol = 1e-9;
1568  belosParams->set("Flexible Gmres", false);
1569  belosParams->set("Num Blocks", 100);
1570  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1571  belosParams->set("Maximum Iterations", 100);
1572  belosParams->set("Verbosity", 33);
1573  belosParams->set("Output Style", 1);
1574  belosParams->set("Output Frequency", 1);
1575  belosParams->set("Output Stream", out.getOStream());
1576  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1577  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1578  problem->setProblem();
1579  Belos::ReturnType ret = solver->solve();
1580  TEST_EQUALITY_CONST( ret, Belos::Converged );
1581 
1582  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1583  // Teuchos::VERB_EXTREME);
1584 
1585  // Check by solving flattened system
1586  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1587  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1588  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1589  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1590  flat_graph =
1591  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1592  cijk_graph, pce_size);
1593  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1594  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1595  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1596  {
1597  RCP<Flat_Tpetra_Vector> flat_x =
1598  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1599  RCP<Flat_Tpetra_Vector> flat_b =
1600  Stokhos::create_flat_vector_view(*b, flat_b_map);
1601  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1602  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1603  typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1604  RCP< FBLinProb > flat_problem =
1605  rcp(new FBLinProb(flat_matrix, flat_x, flat_b));
1606  RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1607  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1608  belosParams));
1609  flat_problem->setProblem();
1610  Belos::ReturnType flat_ret = flat_solver->solve();
1611  TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1612  }
1613 
1614  typename ST::magnitudeType btol = 100*tol;
1615  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1616  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1617  for (size_t i=0; i<num_my_row; ++i) {
1618  TEST_EQUALITY( x_view[i].size(), pce_size );
1619  TEST_EQUALITY( x2_view[i].size(), pce_size );
1620 
1621  // Set small values to zero
1622  Scalar v = x_view[i];
1623  Scalar v2 = x2_view[i];
1624  for (LocalOrdinal j=0; j<pce_size; ++j) {
1625  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
1626  v.fastAccessCoeff(j) = BaseScalar(0.0);
1627  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
1628  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1629  }
1630 
1631  for (LocalOrdinal j=0; j<pce_size; ++j)
1632  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1633  }
1634 
1635  // Clear global tensor
1636  Kokkos::setGlobalCijkTensor(Cijk());
1637 }
1638 
1639 #else
1640 
1641 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1642  Tpetra_CrsMatrix_PCE, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1643 {}
1644 
1645 #endif
1646 
1647 // Test currently doesn't work (in serial) because of our bad division strategy
1648 
1649 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
1650 
1651 //
1652 // Test Belos GMRES solve with Ifpack2 RILUK preconditioning for a
1653 // simple banded upper-triangular matrix
1654 //
1655 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1656  Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
1657 {
1658  using Teuchos::RCP;
1659  using Teuchos::rcp;
1660  using Teuchos::ArrayView;
1661  using Teuchos::Array;
1662  using Teuchos::ArrayRCP;
1663  using Teuchos::ParameterList;
1664 
1665  typedef typename Storage::value_type BaseScalar;
1666  typedef Sacado::UQ::PCE<Storage> Scalar;
1667  typedef typename Scalar::cijk_type Cijk;
1668 
1669  typedef Teuchos::Comm<int> Tpetra_Comm;
1670  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1671  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1672  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1673  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1674 
1675  // Ensure device is initialized
1676  if ( !Kokkos::is_initialized() )
1677  Kokkos::initialize();
1678 
1679  // Cijk
1680  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1681  setGlobalCijkTensor(cijk);
1682  LocalOrdinal pce_size = cijk.dimension();
1683 
1684  // Build banded matrix
1685  GlobalOrdinal nrow = 10;
1686  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1687  RCP<const Tpetra_Map> map =
1688  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1689  nrow, comm);
1690  RCP<Tpetra_CrsGraph> graph =
1691  rcp(new Tpetra_CrsGraph(map, size_t(2)));
1692  Array<GlobalOrdinal> columnIndices(2);
1693  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1694  const size_t num_my_row = myGIDs.size();
1695  for (size_t i=0; i<num_my_row; ++i) {
1696  const GlobalOrdinal row = myGIDs[i];
1697  columnIndices[0] = row;
1698  size_t ncol = 1;
1699  if (row != nrow-1) {
1700  columnIndices[1] = row+1;
1701  ncol = 2;
1702  }
1703  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1704  }
1705  graph->fillComplete();
1706  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1707 
1708  // Set values in matrix
1709  Array<Scalar> vals(2);
1710  Scalar val(cijk);
1711  for (size_t local_row=0; local_row<num_my_row; ++local_row) {
1712  const GlobalOrdinal row = myGIDs[local_row];
1713  const size_t num_col = row == nrow - 1 ? 1 : 2;
1714  for (size_t local_col=0; local_col<num_col; ++local_col) {
1715  const GlobalOrdinal col = row + local_col;
1716  columnIndices[local_col] = col;
1717  for (LocalOrdinal k=0; k<pce_size; ++k)
1718  val.fastAccessCoeff(k) =
1719  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1720  vals[local_col] = val;
1721  }
1722  matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1723  }
1724  matrix->fillComplete();
1725 
1726  // Fill RHS vector
1727  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1728  {
1729  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1730  for (size_t i=0; i<num_my_row; ++i) {
1731  b_view[i] = Scalar(1.0);
1732  }
1733  }
1734 
1735  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1736  RCP<ParameterList> belosParams = rcp(new ParameterList);
1737  typedef BaseScalar BelosScalar;
1738  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1739  typename ST::magnitudeType tol = 1e-9;
1740  {
1741  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1742  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1743  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1744 
1745 #if USE_SCALAR_MEAN_BASED_PREC
1746  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1747  RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1748  Stokhos::build_mean_scalar_matrix(*matrix);
1749  typedef Ifpack2::Preconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Prec;
1750  Ifpack2::Factory factory;
1751  RCP<Scalar_Prec> M_s = factory.create<Scalar_Tpetra_CrsMatrix>("RILUK", mean_matrix);
1752  M_s->initialize();
1753  M_s->compute();
1754  RCP<OP> M = rcp(new Stokhos::MeanBasedTpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(M_s));
1755 #else
1756  RCP<Tpetra_CrsMatrix> mean_matrix = Stokhos::build_mean_matrix(*matrix);
1757  typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
1758  Ifpack2::Factory factory;
1759  RCP<Prec> M = factory.create<Tpetra_CrsMatrix>("RILUK", mean_matrix);
1760  M->initialize();
1761  M->compute();
1762 #endif
1763 
1764  // Solve
1765  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1766  problem->setRightPrec(M);
1767  problem->setProblem();
1768  belosParams->set("Flexible Gmres", false);
1769  belosParams->set("Num Blocks", 100);
1770  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1771  belosParams->set("Maximum Iterations", 100);
1772  belosParams->set("Verbosity", 33);
1773  belosParams->set("Output Style", 1);
1774  belosParams->set("Output Frequency", 1);
1775  belosParams->set("Output Stream", out.getOStream());
1776  //belosParams->set("Orthogonalization", "TSQR");
1777  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1778  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1779  Belos::ReturnType ret = solver->solve();
1780  TEST_EQUALITY_CONST( ret, Belos::Converged );
1781  }
1782 
1783  // Check by solving flattened system
1784  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1785  {
1786  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1787  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1788  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1789  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1790  flat_graph =
1791  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1792  cijk_graph, pce_size);
1793  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1794  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
1795  RCP<Flat_Tpetra_Vector> flat_x =
1796  Stokhos::create_flat_vector_view(*x2, flat_x_map);
1797  RCP<Flat_Tpetra_Vector> flat_b =
1798  Stokhos::create_flat_vector_view(*b, flat_b_map);
1799  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1800  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1801  typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1802  RCP< FBLinProb > flat_problem =
1803  rcp(new FBLinProb(flat_matrix, flat_x, flat_b));
1804  RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1805  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1806  belosParams));
1807  flat_problem->setProblem();
1808  Belos::ReturnType flat_ret = flat_solver->solve();
1809  TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1810  }
1811 
1812  typename ST::magnitudeType btol = 100*tol;
1813  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1814  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1815  for (size_t i=0; i<num_my_row; ++i) {
1816  TEST_EQUALITY( x_view[i].size(), pce_size );
1817  TEST_EQUALITY( x2_view[i].size(), pce_size );
1818 
1819  // Set small values to zero
1820  Scalar v = x_view[i];
1821  Scalar v2 = x2_view[i];
1822  for (LocalOrdinal j=0; j<pce_size; ++j) {
1823  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
1824  v.fastAccessCoeff(j) = BaseScalar(0.0);
1825  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
1826  v2.fastAccessCoeff(j) = BaseScalar(0.0);
1827  }
1828 
1829  for (LocalOrdinal j=0; j<pce_size; ++j)
1830  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
1831  }
1832 
1833  // Clear global tensor
1834  Kokkos::setGlobalCijkTensor(Cijk());
1835 }
1836 
1837 #else
1838 
1839 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1840  Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
1841 {}
1842 
1843 #endif
1844 
1845 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1846 
1847 //
1848 // Test Belos CG solve with MueLu preconditioning for a 1-D Laplacian matrix
1849 //
1850 // Currently requires ETI since the specializations needed for mean-based
1851 // are only brought in with ETI
1852 //
1853 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
1854  Tpetra_CrsMatrix_PCE, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
1855 {
1856  using Teuchos::RCP;
1857  using Teuchos::rcp;
1858  using Teuchos::ArrayView;
1859  using Teuchos::Array;
1860  using Teuchos::ArrayRCP;
1861  using Teuchos::ParameterList;
1862  using Teuchos::getParametersFromXmlFile;
1863 
1864  typedef typename Storage::value_type BaseScalar;
1865  typedef Sacado::UQ::PCE<Storage> Scalar;
1866  typedef typename Scalar::cijk_type Cijk;
1867 
1868  typedef Teuchos::Comm<int> Tpetra_Comm;
1869  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1870  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1871  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1872  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1873 
1874  // Ensure device is initialized
1875  if ( !Kokkos::is_initialized() )
1876  Kokkos::initialize();
1877 
1878  // Cijk
1879  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
1880  setGlobalCijkTensor(cijk);
1881  LocalOrdinal pce_size = cijk.dimension();
1882 
1883  // 1-D Laplacian matrix
1884  GlobalOrdinal nrow = 10;
1885  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1886  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1887  RCP<const Tpetra_Map> map =
1888  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1889  nrow, comm);
1890  RCP<Tpetra_CrsGraph> graph =
1891  rcp(new Tpetra_CrsGraph(map, size_t(3)));
1892  Array<GlobalOrdinal> columnIndices(3);
1893  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1894  const size_t num_my_row = myGIDs.size();
1895  for (size_t i=0; i<num_my_row; ++i) {
1896  const GlobalOrdinal row = myGIDs[i];
1897  if (row == 0 || row == nrow-1) { // Boundary nodes
1898  columnIndices[0] = row;
1899  graph->insertGlobalIndices(row, columnIndices(0,1));
1900  }
1901  else { // Interior nodes
1902  columnIndices[0] = row-1;
1903  columnIndices[1] = row;
1904  columnIndices[2] = row+1;
1905  graph->insertGlobalIndices(row, columnIndices(0,3));
1906  }
1907  }
1908  graph->fillComplete();
1909  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1910 
1911  // Set values in matrix
1912  Array<Scalar> vals(3);
1913  Scalar a_val(cijk);
1914  for (LocalOrdinal j=0; j<pce_size; ++j) {
1915  a_val.fastAccessCoeff(j) =
1916  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(j+1);
1917  }
1918  for (size_t i=0; i<num_my_row; ++i) {
1919  const GlobalOrdinal row = myGIDs[i];
1920  if (row == 0 || row == nrow-1) { // Boundary nodes
1921  columnIndices[0] = row;
1922  vals[0] = Scalar(1.0);
1923  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1924  }
1925  else {
1926  columnIndices[0] = row-1;
1927  columnIndices[1] = row;
1928  columnIndices[2] = row+1;
1929  vals[0] = Scalar(1.0) * a_val;
1930  vals[1] = Scalar(-2.0) * a_val;
1931  vals[2] = Scalar(1.0) * a_val;
1932  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1933  }
1934  }
1935  matrix->fillComplete();
1936 
1937  // Fill RHS vector
1938  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1939  {
1940  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1941  Scalar b_val(cijk);
1942  for (LocalOrdinal j=0; j<pce_size; ++j) {
1943  b_val.fastAccessCoeff(j) =
1944  BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(j+1);
1945  }
1946  for (size_t i=0; i<num_my_row; ++i) {
1947  const GlobalOrdinal row = myGIDs[i];
1948  if (row == 0 || row == nrow-1)
1949  b_view[i] = Scalar(0.0);
1950  else
1951  b_view[i] = b_val * (h*h);
1952  }
1953  }
1954 
1955  // Create preconditioner
1956  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1957  RCP<ParameterList> muelu_params =
1958  getParametersFromXmlFile("muelu_cheby.xml");
1959 #if USE_SCALAR_MEAN_BASED_PREC
1960  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1961  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1962  RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1963  Stokhos::build_mean_scalar_matrix(*matrix);
1964  RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1965  RCP<Scalar_OP> M_s =
1966  MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1967  RCP<OP> M = rcp(new Stokhos::MeanBasedTpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(M_s));
1968 #else
1969  Cijk mean_cijk =
1970  Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
1971  Kokkos::setGlobalCijkTensor(mean_cijk);
1972  RCP<Tpetra_CrsMatrix> mean_matrix = Stokhos::build_mean_matrix(*matrix);
1973  RCP<OP> mean_matrix_op = mean_matrix;
1974  RCP<OP> M =
1975  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1976  Kokkos::setGlobalCijkTensor(cijk);
1977 #endif
1978 
1979  // Solve
1980  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1981  typedef BaseScalar BelosScalar;
1982  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1983  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1984  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1985  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1986  problem->setRightPrec(M);
1987  problem->setProblem();
1988  RCP<ParameterList> belosParams = rcp(new ParameterList);
1989  typename ST::magnitudeType tol = 1e-9;
1990  belosParams->set("Num Blocks", 100);
1991  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1992  belosParams->set("Maximum Iterations", 100);
1993  belosParams->set("Verbosity", 33);
1994  belosParams->set("Output Style", 1);
1995  belosParams->set("Output Frequency", 1);
1996  belosParams->set("Output Stream", out.getOStream());
1997  //belosParams->set("Orthogonalization", "TSQR");
1998  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1999  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2000  // RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2001  // rcp(new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2002  Belos::ReturnType ret = solver->solve();
2003  TEST_EQUALITY_CONST( ret, Belos::Converged );
2004 
2005  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2006  // Teuchos::VERB_EXTREME);
2007 
2008  // Check by solving flattened system
2009  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2010  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2011  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2012  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2013  flat_graph =
2014  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
2015  cijk_graph, pce_size);
2016  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2017  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
2018  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2019  {
2020  RCP<Flat_Tpetra_Vector> flat_x =
2021  Stokhos::create_flat_vector_view(*x2, flat_x_map);
2022  RCP<Flat_Tpetra_Vector> flat_b =
2023  Stokhos::create_flat_vector_view(*b, flat_b_map);
2024  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
2025  typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
2026  typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
2027  RCP< FBLinProb > flat_problem =
2028  rcp(new FBLinProb(flat_matrix, flat_x, flat_b));
2029  RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
2030  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
2031  belosParams));
2032  // RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
2033  // rcp(new Belos::PseudoBlockCGSolMgr<BelosScalar,FMV,FOP>(flat_problem,
2034  // belosParams));
2035  flat_problem->setProblem();
2036  Belos::ReturnType flat_ret = flat_solver->solve();
2037  TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
2038  }
2039 
2040  typename ST::magnitudeType btol = 100*tol;
2041  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2042  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2043  for (size_t i=0; i<num_my_row; ++i) {
2044  TEST_EQUALITY( x_view[i].size(), pce_size );
2045  TEST_EQUALITY( x2_view[i].size(), pce_size );
2046 
2047  // Set small values to zero
2048  Scalar v = x_view[i];
2049  Scalar v2 = x2_view[i];
2050  for (LocalOrdinal j=0; j<pce_size; ++j) {
2051  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
2052  v.fastAccessCoeff(j) = BaseScalar(0.0);
2053  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
2054  v2.fastAccessCoeff(j) = BaseScalar(0.0);
2055  }
2056 
2057  for (LocalOrdinal j=0; j<pce_size; ++j)
2058  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
2059  }
2060 
2061  // Clear global tensor
2062  Kokkos::setGlobalCijkTensor(Cijk());
2063 
2064 }
2065 
2066 #else
2067 
2068 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
2069  Tpetra_CrsMatrix_PCE, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
2070 {}
2071 
2072 #endif
2073 
2074 #if defined(HAVE_STOKHOS_AMESOS2)
2075 
2076 //
2077 // Test Amesos2 solve for a 1-D Laplacian matrix
2078 //
2079 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
2080  Tpetra_CrsMatrix_PCE, Amesos2, Storage, LocalOrdinal, GlobalOrdinal, Node )
2081 {
2082  using Teuchos::RCP;
2083  using Teuchos::rcp;
2084  using Teuchos::ArrayView;
2085  using Teuchos::Array;
2086  using Teuchos::ArrayRCP;
2087  using Teuchos::ParameterList;
2088 
2089  typedef typename Storage::value_type BaseScalar;
2090  typedef Sacado::UQ::PCE<Storage> Scalar;
2091  typedef typename Scalar::cijk_type Cijk;
2092 
2093  typedef Teuchos::Comm<int> Tpetra_Comm;
2094  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2095  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2096  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
2097  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2098  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2099 
2100  // Ensure device is initialized
2101  if ( !Kokkos::is_initialized() )
2102  Kokkos::initialize();
2103 
2104  // Cijk
2105  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
2106  setGlobalCijkTensor(cijk);
2107  LocalOrdinal pce_size = cijk.dimension();
2108 
2109  // 1-D Laplacian matrix
2110  GlobalOrdinal nrow = 10;
2111  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
2112  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2113  RCP<const Tpetra_Map> map =
2114  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2115  nrow, comm);
2116  RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map, size_t(3));
2117  Array<GlobalOrdinal> columnIndices(3);
2118  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2119  const size_t num_my_row = myGIDs.size();
2120  for (size_t i=0; i<num_my_row; ++i) {
2121  const GlobalOrdinal row = myGIDs[i];
2122  if (row == 0 || row == nrow-1) { // Boundary nodes
2123  columnIndices[0] = row;
2124  graph->insertGlobalIndices(row, columnIndices(0,1));
2125  }
2126  else { // Interior nodes
2127  columnIndices[0] = row-1;
2128  columnIndices[1] = row;
2129  columnIndices[2] = row+1;
2130  graph->insertGlobalIndices(row, columnIndices(0,3));
2131  }
2132  }
2133  graph->fillComplete();
2134  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
2135 
2136  // Set values in matrix
2137  Array<Scalar> vals(3);
2138  Scalar a_val(cijk);
2139  for (LocalOrdinal j=0; j<pce_size; ++j) {
2140  a_val.fastAccessCoeff(j) =
2141  BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(j+1);
2142  }
2143  for (size_t i=0; i<num_my_row; ++i) {
2144  const GlobalOrdinal row = myGIDs[i];
2145  if (row == 0 || row == nrow-1) { // Boundary nodes
2146  columnIndices[0] = row;
2147  vals[0] = Scalar(1.0);
2148  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2149  }
2150  else {
2151  columnIndices[0] = row-1;
2152  columnIndices[1] = row;
2153  columnIndices[2] = row+1;
2154  vals[0] = Scalar(1.0) * a_val;
2155  vals[1] = Scalar(-2.0) * a_val;
2156  vals[2] = Scalar(1.0) * a_val;
2157  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2158  }
2159  }
2160  matrix->fillComplete();
2161 
2162  // Fill RHS vector
2163  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2164  {
2165  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2166  Scalar b_val(cijk);
2167  for (LocalOrdinal j=0; j<pce_size; ++j) {
2168  b_val.fastAccessCoeff(j) =
2169  BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(j+1);
2170  }
2171  for (size_t i=0; i<num_my_row; ++i) {
2172  const GlobalOrdinal row = myGIDs[i];
2173  if (row == 0 || row == nrow-1)
2174  b_view[i] = Scalar(0.0);
2175  else
2176  b_view[i] = b_val * (h*h);
2177  }
2178  }
2179 
2180  // Solve
2181  typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2182  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2183  std::string solver_name;
2184 #if defined(HAVE_AMESOS2_BASKER)
2185  solver_name = "basker";
2186 #elif defined(HAVE_AMESOS2_KLU2)
2187  solver_name = "klu2";
2188 #elif defined(HAVE_AMESOS2_SUPERLUDIST)
2189  solver_name = "superlu_dist";
2190 #elif defined(HAVE_AMESOS2_SUPERLUMT)
2191  solver_name = "superlu_mt";
2192 #elif defined(HAVE_AMESOS2_SUPERLU)
2193  solver_name = "superlu";
2194 #elif defined(HAVE_AMESOS2_PARDISO_MKL)
2195  solver_name = "pardisomkl";
2196 #elif defined(HAVE_AMESOS2_LAPACK)
2197  solver_name = "lapack";
2198 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2199  solver_name = "lapack";
2200 #else
2201  // if there are no solvers, we just return as a successful test
2202  success = true;
2203  return;
2204 #endif
2205  out << "Solving linear system with " << solver_name << std::endl;
2206  {
2207  RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2208  solver_name, matrix, x, b);
2209  solver->solve();
2210  }
2211  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2212  // Teuchos::VERB_EXTREME);
2213 
2214  // Check by solving flattened system
2215  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2216  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_MultiVector;
2217  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2218  RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2219  RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2220  flat_graph =
2221  Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
2222  cijk_graph, pce_size);
2223  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2224  Stokhos::create_flat_matrix(*matrix, flat_graph, cijk_graph, cijk);
2225  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2226  {
2227  RCP<Flat_Tpetra_Vector> flat_x =
2228  Stokhos::create_flat_vector_view(*x2, flat_x_map);
2229  RCP<Flat_Tpetra_Vector> flat_b =
2230  Stokhos::create_flat_vector_view(*b, flat_b_map);
2231  typedef Amesos2::Solver<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector> Flat_Solver;
2232  RCP<Flat_Solver> flat_solver =
2233  Amesos2::create<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector>(
2234  solver_name, flat_matrix, flat_x, flat_b);
2235  flat_solver->solve();
2236  }
2237 
2238  typedef Kokkos::ArithTraits<BaseScalar> ST;
2239  typename ST::mag_type btol = 1e-12;
2240  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2241  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2242  for (size_t i=0; i<num_my_row; ++i) {
2243  TEST_EQUALITY( x_view[i].size(), pce_size );
2244  TEST_EQUALITY( x2_view[i].size(), pce_size );
2245 
2246  // Set small values to zero
2247  Scalar v = x_view[i];
2248  Scalar v2 = x2_view[i];
2249  for (LocalOrdinal j=0; j<pce_size; ++j) {
2250  if (j < v.size() && ST::magnitude(v.coeff(j)) < btol)
2251  v.fastAccessCoeff(j) = BaseScalar(0.0);
2252  if (j < v2.size() && ST::magnitude(v2.coeff(j)) < btol)
2253  v2.fastAccessCoeff(j) = BaseScalar(0.0);
2254  }
2255 
2256  for (LocalOrdinal j=0; j<pce_size; ++j)
2257  TEST_FLOATING_EQUALITY(v.coeff(j), v2.coeff(j), btol);
2258  }
2259 
2260  // Clear global tensor
2261  Kokkos::setGlobalCijkTensor(Cijk());
2262 }
2263 
2264 #else
2265 
2266 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
2267  Tpetra_CrsMatrix_PCE, Amesos2, Storage, LocalOrdinal, GlobalOrdinal, Node )
2268 {}
2269 
2270 #endif
2271 
2272 #define CRSMATRIX_UQ_PCE_TESTS_SLGN(S, LO, GO, N) \
2273  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorAdd, S, LO, GO, N ) \
2274  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorDot, S, LO, GO, N ) \
2275  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorAdd, S, LO, GO, N ) \
2276  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDot, S, LO, GO, N ) \
2277  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDotSub, S, LO, GO, N ) \
2278  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixVectorMultiply, S, LO, GO, N ) \
2279  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2280  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Flatten, S, LO, GO, N ) \
2281  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimpleCG, S, LO, GO, N ) \
2282  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, S, LO, GO, N ) \
2283  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES, S, LO, GO, N ) \
2284  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, S, LO, GO, N ) \
2285  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosCG_Muelu, S, LO, GO, N ) \
2286  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Amesos2, S, LO, GO, N )
2287 
2288 #define CRSMATRIX_UQ_PCE_TESTS_N(N) \
2289  typedef Stokhos::DeviceForNode2<N>::type Device; \
2290  typedef Stokhos::DynamicStorage<int,double,Device::execution_space> DS; \
2291  using default_local_ordinal_type = Tpetra::Map<>::local_ordinal_type; \
2292  using default_global_ordinal_type = Tpetra::Map<>::global_ordinal_type; \
2293  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
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:656
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)
Kokkos::create_mirror_view
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Definition: Stokhos_CrsMatrix.hpp:645
Generated on Fri Nov 22 2024 09:18:15 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5