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