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