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