Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_TpetraCrsMatrixMPVectorUnitTest.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
12 #include "Stokhos_Ensemble_Sizes.hpp"
13 
14 // Teuchos
16 
17 // Tpetra
20 #include "Tpetra_Core.hpp"
21 #include "Tpetra_Map.hpp"
22 #include "Tpetra_MultiVector.hpp"
23 #include "Tpetra_Vector.hpp"
24 #include "Tpetra_CrsGraph.hpp"
25 #include "Tpetra_CrsMatrix.hpp"
26 #include "Tpetra_Details_WrappedDualView.hpp"
27 #include "Stokhos_Tpetra_CG.hpp"
28 
29 // Belos solver
30 #ifdef HAVE_STOKHOS_BELOS
32 #include "BelosLinearProblem.hpp"
33 #include "BelosPseudoBlockGmresSolMgr.hpp"
34 #include "BelosPseudoBlockCGSolMgr.hpp"
35 #endif
36 
37 // Ifpack2 preconditioner
38 #ifdef HAVE_STOKHOS_IFPACK2
40 #include "Ifpack2_Factory.hpp"
41 #endif
42 
43 // MueLu preconditioner
44 #ifdef HAVE_STOKHOS_MUELU
46 #include "MueLu_CreateTpetraPreconditioner.hpp"
47 #endif
48 
49 // Amesos2 solver
50 #ifdef HAVE_STOKHOS_AMESOS2
52 #include "Amesos2_Factory.hpp"
53 #endif
54 
55 template <typename scalar, typename ordinal>
56 inline
57 scalar generate_vector_coefficient( const ordinal nFEM,
58  const ordinal nStoch,
59  const ordinal iColFEM,
60  const ordinal iStoch )
61 {
62  const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
63  const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
64  return X_fem + X_stoch;
65  //return 1.0;
66 }
67 
68 template <typename scalar, typename ordinal>
69 inline
70 scalar generate_multi_vector_coefficient( const ordinal nFEM,
71  const ordinal nVec,
72  const ordinal nStoch,
73  const ordinal iColFEM,
74  const ordinal iVec,
75  const ordinal iStoch)
76 {
77  const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
78  const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
79  const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
80  return X_fem + X_stoch + X_col;
81  //return 1.0;
82 }
83 
84 //
85 // Tests
86 //
87 
88 const int VectorSize = STOKHOS_DEFAULT_ENSEMBLE_SIZE;
89 
90 //
91 // Test vector addition
92 //
94  Tpetra_CrsMatrix_MP, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node )
95 {
96  using Teuchos::RCP;
97  using Teuchos::rcp;
98  using Teuchos::ArrayView;
99  using Teuchos::Array;
100  using Teuchos::ArrayRCP;
101 
102  typedef typename Storage::value_type BaseScalar;
104 
105  typedef Teuchos::Comm<int> Tpetra_Comm;
106  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
107  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
108 
109  // Ensure device is initialized
110  if ( !Kokkos::is_initialized() )
111  Kokkos::initialize();
112 
113  // Comm
114  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
115 
116  // Map
117  GlobalOrdinal nrow = 10;
118  RCP<const Tpetra_Map> map =
119  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
120  nrow, comm);
121  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
122  const size_t num_my_row = myGIDs.size();
123 
124  // Fill vectors
125  RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
126  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
127  {
128  auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
129  auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
130  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
131  for (size_t i=0; i<num_my_row; ++i) {
132  const GlobalOrdinal row = myGIDs[i];
133  for (LocalOrdinal j=0; j<VectorSize; ++j) {
134  val1.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
135  val2.fastAccessCoeff(j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
136  }
137  x1_view(i,0) = val1;
138  x2_view(i,0) = val2;
139  }
140  }
141 
142  // Add
143  Scalar alpha = 2.1;
144  Scalar beta = 3.7;
145  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
146  y->update(alpha, *x1, beta, *x2, Scalar(0.0));
147 
148  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
149  // Teuchos::VERB_EXTREME);
150 
151  // Check
152  auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
153  Scalar val(VectorSize, BaseScalar(0.0));
154  BaseScalar tol = 1.0e-14;
155  for (size_t i=0; i<num_my_row; ++i) {
156  const GlobalOrdinal row = myGIDs[i];
157  for (LocalOrdinal j=0; j<VectorSize; ++j) {
158  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
159  nrow, VectorSize, row, j);
160  val.fastAccessCoeff(j) = alpha.coeff(j)*v + 0.12345*beta.coeff(j)*v;
161  }
162  TEST_EQUALITY( y_view(i,0).size(), VectorSize );
163  for (LocalOrdinal j=0; j<VectorSize; ++j)
164  TEST_FLOATING_EQUALITY( y_view(i,0).fastAccessCoeff(j), val.fastAccessCoeff(j), tol );
165  }
166 }
167 
168 //
169 // Test vector dot product
170 //
172  Tpetra_CrsMatrix_MP, VectorDot, Storage, LocalOrdinal, GlobalOrdinal, Node )
173 {
174  using Teuchos::RCP;
175  using Teuchos::rcp;
176  using Teuchos::ArrayView;
177  using Teuchos::Array;
178  using Teuchos::ArrayRCP;
179 
180  typedef typename Storage::value_type BaseScalar;
182 
183  typedef Teuchos::Comm<int> Tpetra_Comm;
184  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
185  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
186  typedef typename Tpetra_Vector::dot_type dot_type;
187 
188  // Ensure device is initialized
189  if ( !Kokkos::is_initialized() )
190  Kokkos::initialize();
191 
192  // Comm
193  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
194 
195  // Map
196  GlobalOrdinal nrow = 10;
197  RCP<const Tpetra_Map> map =
198  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
199  nrow, comm);
200  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
201  const size_t num_my_row = myGIDs.size();
202 
203  // Fill vectors
204  RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
205  RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
206  {
207  auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
208  auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
209  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
210  for (size_t i=0; i<num_my_row; ++i) {
211  const GlobalOrdinal row = myGIDs[i];
212  for (LocalOrdinal j=0; j<VectorSize; ++j) {
213  val1.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
214  val2.fastAccessCoeff(j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
215  }
216  x1_view(i,0) = val1;
217  x2_view(i,0) = val2;
218  }
219  }
220 
221  // Dot product
222  dot_type dot = x1->dot(*x2);
223 
224  // Check
225 
226 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
227 
228  // Local contribution
229  dot_type local_val(0);
230  for (size_t i=0; i<num_my_row; ++i) {
231  const GlobalOrdinal row = myGIDs[i];
232  for (LocalOrdinal j=0; j<VectorSize; ++j) {
233  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
234  nrow, VectorSize, row, j);
235  local_val += 0.12345 * v * v;
236  }
237  }
238 
239  // Global reduction
240  dot_type val(0);
241  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
242  Teuchos::outArg(val));
243 
244 #else
245 
246  // Local contribution
247  dot_type local_val(VectorSize, 0.0);
248  for (size_t i=0; i<num_my_row; ++i) {
249  const GlobalOrdinal row = myGIDs[i];
250  for (LocalOrdinal j=0; j<VectorSize; ++j) {
251  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
252  nrow, VectorSize, row, j);
253  local_val.fastAccessCoeff(j) += 0.12345 * v * v;
254  }
255  }
256 
257  // Global reduction
258  dot_type val(VectorSize, 0.0);
259  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
260  Teuchos::outArg(val));
261 
262 #endif
263 
264  out << "dot = " << dot << " expected = " << val << std::endl;
265 
266  BaseScalar tol = 1.0e-14;
267  TEST_FLOATING_EQUALITY( dot, val, tol );
268 }
269 
270 //
271 // Test multi-vector addition
272 //
274  Tpetra_CrsMatrix_MP, MultiVectorAdd, 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;
284 
285  typedef Teuchos::Comm<int> Tpetra_Comm;
286  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
287  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
288 
289  // Ensure device is initialized
290  if ( !Kokkos::is_initialized() )
291  Kokkos::initialize();
292 
293  // Comm
294  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
295 
296  // Map
297  GlobalOrdinal nrow = 10;
298  RCP<const Tpetra_Map> map =
299  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
300  nrow, comm);
301  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
302  const size_t num_my_row = myGIDs.size();
303 
304  // Fill vectors
305  size_t ncol = 5;
306  RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
307  RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
308  {
309  auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
310  auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
311  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
312  for (size_t i=0; i<num_my_row; ++i) {
313  const GlobalOrdinal row = myGIDs[i];
314  for (size_t j=0; j<ncol; ++j) {
315  for (LocalOrdinal k=0; k<VectorSize; ++k) {
316  BaseScalar v =
317  generate_multi_vector_coefficient<BaseScalar,size_t>(
318  nrow, ncol, VectorSize, row, j, k);
319  val1.fastAccessCoeff(k) = v;
320  val2.fastAccessCoeff(k) = 0.12345 * v;
321  }
322  x1_view(i,j) = val1;
323  x2_view(i,j) = val2;
324  }
325  }
326  }
327 
328  // Add
329  Scalar alpha = 2.1;
330  Scalar beta = 3.7;
331  RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
332  y->update(alpha, *x1, beta, *x2, Scalar(0.0));
333 
334  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
335  // Teuchos::VERB_EXTREME);
336 
337  // Check
338  auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
339  Scalar val(VectorSize, BaseScalar(0.0));
340  BaseScalar tol = 1.0e-14;
341  for (size_t i=0; i<num_my_row; ++i) {
342  const GlobalOrdinal row = myGIDs[i];
343  for (size_t j=0; j<ncol; ++j) {
344  for (LocalOrdinal k=0; k<VectorSize; ++k) {
345  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
346  nrow, ncol, VectorSize, row, j, k);
347  val.fastAccessCoeff(k) = alpha.coeff(k)*v + 0.12345*beta.coeff(k)*v;
348  }
349  TEST_EQUALITY( y_view(i,j).size(), VectorSize );
350  for (LocalOrdinal k=0; k<VectorSize; ++k)
352  val.fastAccessCoeff(k), tol );
353  }
354  }
355 }
356 
357 //
358 // Test multi-vector dot product
359 //
361  Tpetra_CrsMatrix_MP, MultiVectorDot, Storage, LocalOrdinal, GlobalOrdinal, Node )
362 {
363  using Teuchos::RCP;
364  using Teuchos::rcp;
365  using Teuchos::ArrayView;
366  using Teuchos::Array;
367  using Teuchos::ArrayRCP;
368 
369  typedef typename Storage::value_type BaseScalar;
371 
372  typedef Teuchos::Comm<int> Tpetra_Comm;
373  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
374  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
375  typedef typename Tpetra_MultiVector::dot_type dot_type;
376 
377  // Ensure device is initialized
378  if ( !Kokkos::is_initialized() )
379  Kokkos::initialize();
380 
381  // Comm
382  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
383 
384  // Map
385  GlobalOrdinal nrow = 10;
386  RCP<const Tpetra_Map> map =
387  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
388  nrow, comm);
389  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
390  const size_t num_my_row = myGIDs.size();
391 
392  // Fill vectors
393  size_t ncol = 5;
394  RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
395  RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
396  {
397  auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
398  auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
399  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
400  for (size_t i=0; i<num_my_row; ++i) {
401  const GlobalOrdinal row = myGIDs[i];
402  for (size_t j=0; j<ncol; ++j) {
403  for (LocalOrdinal k=0; k<VectorSize; ++k) {
404  BaseScalar v =
405  generate_multi_vector_coefficient<BaseScalar,size_t>(
406  nrow, ncol, VectorSize, row, j, k);
407  val1.fastAccessCoeff(k) = v;
408  val2.fastAccessCoeff(k) = 0.12345 * v;
409  }
410  x1_view(i,j) = val1;
411  x2_view(i,j) = val2;
412  }
413  }
414  }
415 
416  // Dot product
417  Array<dot_type> dots(ncol);
418  x1->dot(*x2, dots());
419 
420  // Check
421 
422 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
423 
424  // Local contribution
425  Array<dot_type> local_vals(ncol, dot_type(0));
426  for (size_t i=0; i<num_my_row; ++i) {
427  const GlobalOrdinal row = myGIDs[i];
428  for (size_t j=0; j<ncol; ++j) {
429  for (LocalOrdinal k=0; k<VectorSize; ++k) {
430  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
431  nrow, ncol, VectorSize, row, j, k);
432  local_vals[j] += 0.12345 * v * v;
433  }
434  }
435  }
436 
437  // Global reduction
438  Array<dot_type> vals(ncol, dot_type(0));
439  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
440  local_vals.getRawPtr(), vals.getRawPtr());
441 
442 #else
443 
444  // Local contribution
445  Array<dot_type> local_vals(ncol, dot_type(VectorSize, 0.0));
446  for (size_t i=0; i<num_my_row; ++i) {
447  const GlobalOrdinal row = myGIDs[i];
448  for (size_t j=0; j<ncol; ++j) {
449  for (LocalOrdinal k=0; k<VectorSize; ++k) {
450  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
451  nrow, ncol, VectorSize, row, j, k);
452  local_vals[j].fastAccessCoeff(k) += 0.12345 * v * v;
453  }
454  }
455  }
456 
457  // Global reduction
458  Array<dot_type> vals(ncol, dot_type(VectorSize, 0.0));
459  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
460  local_vals.getRawPtr(), vals.getRawPtr());
461 
462 #endif
463 
464  BaseScalar tol = 1.0e-14;
465  for (size_t j=0; j<ncol; ++j) {
466  out << "dots(" << j << ") = " << dots[j]
467  << " expected(" << j << ") = " << vals[j] << std::endl;
468  TEST_FLOATING_EQUALITY( dots[j], vals[j], tol );
469  }
470 }
471 
472 //
473 // Test multi-vector dot product using subviews
474 //
476  Tpetra_CrsMatrix_MP, MultiVectorDotSub, Storage, LocalOrdinal, GlobalOrdinal, Node )
477 {
478  using Teuchos::RCP;
479  using Teuchos::rcp;
480  using Teuchos::ArrayView;
481  using Teuchos::Array;
482  using Teuchos::ArrayRCP;
483 
484  typedef typename Storage::value_type BaseScalar;
486 
487  typedef Teuchos::Comm<int> Tpetra_Comm;
488  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
489  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
490  typedef typename Tpetra_MultiVector::dot_type dot_type;
491 
492  // Ensure device is initialized
493  if ( !Kokkos::is_initialized() )
494  Kokkos::initialize();
495 
496  // Comm
497  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
498 
499  // Map
500  GlobalOrdinal nrow = 10;
501  RCP<const Tpetra_Map> map =
502  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
503  nrow, comm);
504  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
505  const size_t num_my_row = myGIDs.size();
506 
507  // Fill vectors
508  size_t ncol = 5;
509  RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
510  RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
511  {
512  auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
513  auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
514  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
515  for (size_t i=0; i<num_my_row; ++i) {
516  const GlobalOrdinal row = myGIDs[i];
517  for (size_t j=0; j<ncol; ++j) {
518  for (LocalOrdinal k=0; k<VectorSize; ++k) {
519  BaseScalar v =
520  generate_multi_vector_coefficient<BaseScalar,size_t>(
521  nrow, ncol, VectorSize, row, j, k);
522  val1.fastAccessCoeff(k) = v;
523  val2.fastAccessCoeff(k) = 0.12345 * v;
524  }
525  x1_view(i,j) = val1;
526  x2_view(i,j) = val2;
527  }
528  }
529  }
530 
531  // Get subviews
532  size_t ncol_sub = 2;
533  Teuchos::Array<size_t> cols(ncol_sub);
534  cols[0] = 4; cols[1] = 2;
535  RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
536  RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
537 
538  // Dot product
539  Array<dot_type> dots(ncol_sub);
540  x1_sub->dot(*x2_sub, dots());
541 
542  // Check
543 
544 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
545 
546  // Local contribution
547  Array<dot_type> local_vals(ncol_sub, dot_type(0));
548  for (size_t i=0; i<num_my_row; ++i) {
549  const GlobalOrdinal row = myGIDs[i];
550  for (size_t j=0; j<ncol_sub; ++j) {
551  for (LocalOrdinal k=0; k<VectorSize; ++k) {
552  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
553  nrow, ncol, VectorSize, row, cols[j], k);
554  local_vals[j] += 0.12345 * v * v;
555  }
556  }
557  }
558 
559  // Global reduction
560  Array<dot_type> vals(ncol_sub, dot_type(0));
561  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
562  Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
563  vals.getRawPtr());
564 
565 #else
566 
567  // Local contribution
568  Array<dot_type> local_vals(ncol_sub, dot_type(VectorSize, 0.0));
569  for (size_t i=0; i<num_my_row; ++i) {
570  const GlobalOrdinal row = myGIDs[i];
571  for (size_t j=0; j<ncol_sub; ++j) {
572  for (LocalOrdinal k=0; k<VectorSize; ++k) {
573  BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
574  nrow, ncol, VectorSize, row, cols[j], k);
575  local_vals[j].fastAccessCoeff(k) += 0.12345 * v * v;
576  }
577  }
578  }
579 
580  // Global reduction
581  Array<dot_type> vals(ncol_sub, dot_type(VectorSize, 0.0));
582  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
583  Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
584  vals.getRawPtr());
585 
586 #endif
587 
588  BaseScalar tol = 1.0e-14;
589  for (size_t j=0; j<ncol_sub; ++j) {
590  out << "dots(" << j << ") = " << dots[j]
591  << " expected(" << j << ") = " << vals[j] << std::endl;
592  TEST_FLOATING_EQUALITY( dots[j], vals[j], tol );
593  }
594 }
595 
596 //
597 // Test matrix-vector multiplication for a simple banded upper-triangular matrix
598 //
600  Tpetra_CrsMatrix_MP, MatrixVectorMultiply, Storage, LocalOrdinal, GlobalOrdinal, Node )
601 {
602  using Teuchos::RCP;
603  using Teuchos::rcp;
604  using Teuchos::ArrayView;
605  using Teuchos::Array;
606  using Teuchos::ArrayRCP;
607 
608  typedef typename Storage::value_type BaseScalar;
610 
611  typedef Teuchos::Comm<int> Tpetra_Comm;
612  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
613  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
614  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
615  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
616 
617  // Ensure device is initialized
618  if ( !Kokkos::is_initialized() )
619  Kokkos::initialize();
620 
621  // Build banded matrix
622  GlobalOrdinal nrow = 10;
623  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
624  RCP<const Tpetra_Map> map =
625  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
626  nrow, comm);
627  RCP<Tpetra_CrsGraph> graph =
628  rcp(new Tpetra_CrsGraph(map, size_t(2)));
629  Array<GlobalOrdinal> columnIndices(2);
630  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
631  const size_t num_my_row = myGIDs.size();
632  for (size_t i=0; i<num_my_row; ++i) {
633  const GlobalOrdinal row = myGIDs[i];
634  columnIndices[0] = row;
635  size_t ncol = 1;
636  if (row != nrow-1) {
637  columnIndices[1] = row+1;
638  ncol = 2;
639  }
640  graph->insertGlobalIndices(row, columnIndices(0,ncol));
641  }
642  graph->fillComplete();
643  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
644 
645  // Set values in matrix
646  Array<Scalar> vals(2);
647  Scalar val(VectorSize, BaseScalar(0.0));
648  for (size_t i=0; i<num_my_row; ++i) {
649  const GlobalOrdinal row = myGIDs[i];
650  columnIndices[0] = row;
651  for (LocalOrdinal j=0; j<VectorSize; ++j)
652  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
653  nrow, VectorSize, row, j);
654  vals[0] = val;
655  size_t ncol = 1;
656 
657  if (row != nrow-1) {
658  columnIndices[1] = row+1;
659  for (LocalOrdinal j=0; j<VectorSize; ++j)
660  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
661  nrow, VectorSize, row+1, j);
662  vals[1] = val;
663  ncol = 2;
664  }
665  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
666  }
667  matrix->fillComplete();
668 
669  // Fill vector
670  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
671  {
672  auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
673  for (size_t i=0; i<num_my_row; ++i) {
674  const GlobalOrdinal row = myGIDs[i];
675  for (LocalOrdinal j=0; j<VectorSize; ++j)
676  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
677  nrow, VectorSize, row, j);
678  x_view(i,0) = val;
679  }
680  }
681 
682  // matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
683  // Teuchos::VERB_EXTREME);
684 
685  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
686  // Teuchos::VERB_EXTREME);
687 
688  // Multiply
689  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
690  matrix->apply(*x, *y);
691 
692  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
693  // Teuchos::VERB_EXTREME);
694 
695  // Check
696  auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
697  BaseScalar tol = 1.0e-14;
698  for (size_t i=0; i<num_my_row; ++i) {
699  const GlobalOrdinal row = myGIDs[i];
700  for (LocalOrdinal j=0; j<VectorSize; ++j) {
701  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
702  nrow, VectorSize, row, j);
703  val.fastAccessCoeff(j) = v*v;
704  }
705  if (row != nrow-1) {
706  for (LocalOrdinal j=0; j<VectorSize; ++j) {
707  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
708  nrow, VectorSize, row+1, j);
709  val.fastAccessCoeff(j) += v*v;
710  }
711  }
712  TEST_EQUALITY( y_view(i,0).size(), VectorSize );
713  for (LocalOrdinal j=0; j<VectorSize; ++j)
714  TEST_FLOATING_EQUALITY( y_view(i,0).fastAccessCoeff(j), val.fastAccessCoeff(j), tol );
715  }
716 }
717 
718 //
719 // Test matrix-multi-vector multiplication for a simple banded upper-triangular matrix
720 //
722  Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, Storage, LocalOrdinal, GlobalOrdinal, Node )
723 {
724  using Teuchos::RCP;
725  using Teuchos::rcp;
726  using Teuchos::ArrayView;
727  using Teuchos::Array;
728  using Teuchos::ArrayRCP;
729 
730  typedef typename Storage::value_type BaseScalar;
732 
733  typedef Teuchos::Comm<int> Tpetra_Comm;
734  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
735  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
736  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
737  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
738 
739  // Ensure device is initialized
740  if ( !Kokkos::is_initialized() )
741  Kokkos::initialize();
742 
743  // Build banded matrix
744  GlobalOrdinal nrow = 10;
745  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
746  RCP<const Tpetra_Map> map =
747  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
748  nrow, comm);
749  RCP<Tpetra_CrsGraph> graph =
750  rcp(new Tpetra_CrsGraph(map, size_t(2)));
751  Array<GlobalOrdinal> columnIndices(2);
752  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
753  const size_t num_my_row = myGIDs.size();
754  for (size_t i=0; i<num_my_row; ++i) {
755  const GlobalOrdinal row = myGIDs[i];
756  columnIndices[0] = row;
757  size_t ncol = 1;
758  if (row != nrow-1) {
759  columnIndices[1] = row+1;
760  ncol = 2;
761  }
762  graph->insertGlobalIndices(row, columnIndices(0,ncol));
763  }
764  graph->fillComplete();
765  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
766 
767  // Set values in matrix
768  Array<Scalar> vals(2);
769  Scalar val(VectorSize, BaseScalar(0.0));
770  for (size_t i=0; i<num_my_row; ++i) {
771  const GlobalOrdinal row = myGIDs[i];
772  columnIndices[0] = row;
773  for (LocalOrdinal j=0; j<VectorSize; ++j)
774  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
775  nrow, VectorSize, row, j);
776  vals[0] = val;
777  size_t ncol = 1;
778 
779  if (row != nrow-1) {
780  columnIndices[1] = row+1;
781  for (LocalOrdinal j=0; j<VectorSize; ++j)
782  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
783  nrow, VectorSize, row+1, j);
784  vals[1] = val;
785  ncol = 2;
786  }
787  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
788  }
789  matrix->fillComplete();
790 
791  // Fill multi-vector
792  size_t ncol = 5;
793  RCP<Tpetra_MultiVector> x = Tpetra::createMultiVector<Scalar>(map, ncol);
794  {
795  auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
796  for (size_t i=0; i<num_my_row; ++i) {
797  const GlobalOrdinal row = myGIDs[i];
798  for (size_t j=0; j<ncol; ++j) {
799  for (LocalOrdinal k=0; k<VectorSize; ++k) {
800  BaseScalar v =
801  generate_multi_vector_coefficient<BaseScalar,size_t>(
802  nrow, ncol, VectorSize, row, j, k);
803  val.fastAccessCoeff(k) = v;
804  }
805  x_view(i,j) = val;
806  }
807  }
808  }
809 
810  // matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
811  // Teuchos::VERB_EXTREME);
812 
813  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
814  // Teuchos::VERB_EXTREME);
815 
816  // Multiply
817  RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
818  matrix->apply(*x, *y);
819 
820  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
821  // Teuchos::VERB_EXTREME);
822 
823  // Check
824  auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
825  BaseScalar tol = 1.0e-14;
826  for (size_t i=0; i<num_my_row; ++i) {
827  const GlobalOrdinal row = myGIDs[i];
828  for (size_t j=0; j<ncol; ++j) {
829  for (LocalOrdinal k=0; k<VectorSize; ++k) {
830  BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
831  nrow, VectorSize, row, k);
832  BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
833  nrow, ncol, VectorSize, row, j, k);
834  val.fastAccessCoeff(k) = v1*v2;
835  }
836  if (row != nrow-1) {
837  for (LocalOrdinal k=0; k<VectorSize; ++k) {
838  BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
839  nrow, VectorSize, row+1, k);
840  BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
841  nrow, ncol, VectorSize, row+1, j, k);
842  val.fastAccessCoeff(k) += v1*v2;
843  }
844  }
845  TEST_EQUALITY( y_view(i,j).size(), VectorSize );
846  for (LocalOrdinal k=0; k<VectorSize; ++k)
848  val.fastAccessCoeff(k), tol );
849  }
850  }
851 }
852 
853 //
854 // Test flattening MP::Vector matrix
855 //
857  Tpetra_CrsMatrix_MP, Flatten, Storage, LocalOrdinal, GlobalOrdinal, Node )
858 {
859  using Teuchos::RCP;
860  using Teuchos::rcp;
861  using Teuchos::ArrayView;
862  using Teuchos::Array;
863  using Teuchos::ArrayRCP;
864 
865  typedef typename Storage::value_type BaseScalar;
867 
868  typedef Teuchos::Comm<int> Tpetra_Comm;
869  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
870  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
871  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
872  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
873 
874  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
875  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
876 
877  // Ensure device is initialized
878  if ( !Kokkos::is_initialized() )
879  Kokkos::initialize();
880 
881  // Build banded matrix
882  GlobalOrdinal nrow = 10;
883  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
884  RCP<const Tpetra_Map> map =
885  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
886  nrow, comm);
887  RCP<Tpetra_CrsGraph> graph =
888  rcp(new Tpetra_CrsGraph(map, size_t(2)));
889  Array<GlobalOrdinal> columnIndices(2);
890  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
891  const size_t num_my_row = myGIDs.size();
892  for (size_t i=0; i<num_my_row; ++i) {
893  const GlobalOrdinal row = myGIDs[i];
894  columnIndices[0] = row;
895  size_t ncol = 1;
896  if (row != nrow-1) {
897  columnIndices[1] = row+1;
898  ncol = 2;
899  }
900  graph->insertGlobalIndices(row, columnIndices(0,ncol));
901  }
902  graph->fillComplete();
903  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
904 
905  // Set values in matrix
906  Array<Scalar> vals(2);
907  Scalar val(VectorSize, BaseScalar(0.0));
908  for (size_t i=0; i<num_my_row; ++i) {
909  const GlobalOrdinal row = myGIDs[i];
910  columnIndices[0] = row;
911  for (LocalOrdinal j=0; j<VectorSize; ++j)
912  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
913  nrow, VectorSize, row, j);
914  vals[0] = val;
915  size_t ncol = 1;
916 
917  if (row != nrow-1) {
918  columnIndices[1] = row+1;
919  for (LocalOrdinal j=0; j<VectorSize; ++j)
920  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
921  nrow, VectorSize, row+1, j);
922  vals[1] = val;
923  ncol = 2;
924  }
925  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
926  }
927  matrix->fillComplete();
928 
929  // Fill vector
930  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
931  {
932  auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
933  for (size_t i=0; i<num_my_row; ++i) {
934  const GlobalOrdinal row = myGIDs[i];
935  for (LocalOrdinal j=0; j<VectorSize; ++j)
936  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
937  nrow, VectorSize, row, j);
938  x_view(i,0) = val;
939  }
940  }
941 
942  // Multiply
943  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
944  matrix->apply(*x, *y);
945 
946  // Flatten matrix
947  RCP<const Tpetra_Map> flat_x_map, flat_y_map;
948  RCP<const Tpetra_CrsGraph> flat_graph =
949  Stokhos::create_flat_mp_graph(*graph, flat_x_map, flat_y_map, VectorSize);
950  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
951  Stokhos::create_flat_matrix(*matrix, flat_graph, VectorSize);
952 
953  // Multiply with flattened matix
954  RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
955  {
956  RCP<Flat_Tpetra_Vector> flat_x =
957  Stokhos::create_flat_vector_view(*x, flat_x_map);
958  RCP<Flat_Tpetra_Vector> flat_y =
959  Stokhos::create_flat_vector_view(*y2, flat_y_map);
960  flat_matrix->apply(*flat_x, *flat_y);
961  }
962 
963  // flat_y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
964  // Teuchos::VERB_EXTREME);
965 
966  // Check
967  BaseScalar tol = 1.0e-14;
968  auto y_view = y-> getLocalViewHost(Tpetra::Access::ReadOnly);
969  auto y2_view = y2->getLocalViewHost(Tpetra::Access::ReadOnly);
970  for (size_t i=0; i<num_my_row; ++i) {
971  TEST_EQUALITY( y_view( i,0).size(), VectorSize );
972  TEST_EQUALITY( y2_view(i,0).size(), VectorSize );
973  for (LocalOrdinal j=0; j<VectorSize; ++j)
975  y2_view(i,0).fastAccessCoeff(j), tol );
976  }
977 }
978 
979 //
980 // Test interaction between Tpetra WrappedDualView and MP::Vector
981 //
983  Tpetra_CrsMatrix_MP, WrappedDualView, Storage, LocalOrdinal, GlobalOrdinal, Node )
984 {
985  //BMK 6-2021: This test is required because a View of MP::Vector has slightly different behavior than a typical Kokkos::View.
986  //If you construct a Kokkos::View with a label and 0 extent, it gets a non-null allocation.
987  //But for View<MP::Vector>, the same constructor produces a null data pointer but
988  //an active reference counting node (use_count() > 0).
989  //This test makes sure that Tpetra WrappedDualView works correctly with a View where data() == nullptr but use_count() > 0.
990  using Teuchos::RCP;
991  using Teuchos::rcp;
992  using Teuchos::ArrayView;
993  using Teuchos::Array;
994  using Teuchos::ArrayRCP;
995 
996  //typedef typename Storage::value_type BaseScalar;
998 
999  using DualViewType = Kokkos::DualView<Scalar*, typename Node::device_type>;
1000  using WDV = Tpetra::Details::WrappedDualView<DualViewType>;
1001  using values_view = typename DualViewType::t_dev;
1002 
1003  // Ensure device is initialized
1004  if ( !Kokkos::is_initialized() )
1005  Kokkos::initialize();
1006 
1007  WDV wdv;
1008  {
1009  values_view myView("emptyTestView", 0);
1010  wdv = WDV(myView);
1011  }
1012  size_t use_h = wdv.getHostView(Tpetra::Access::ReadOnly).use_count();
1013  size_t use_d = wdv.getDeviceView(Tpetra::Access::ReadOnly).use_count();
1014  //The WrappedDualView is now the only object holding references to the host and device views,
1015  //so they should have identical use counts.
1016  TEST_EQUALITY(use_h, use_d);
1017 }
1018 
1019 //
1020 // Test simple CG solve without preconditioning for a 1-D Laplacian matrix
1021 //
1023  Tpetra_CrsMatrix_MP, SimpleCG, Storage, LocalOrdinal, GlobalOrdinal, Node )
1024 {
1025  using Teuchos::RCP;
1026  using Teuchos::rcp;
1027  using Teuchos::ArrayView;
1028  using Teuchos::Array;
1029  using Teuchos::ArrayRCP;
1030  using Teuchos::ParameterList;
1031 
1032  typedef typename Storage::value_type BaseScalar;
1034 
1035  typedef Teuchos::Comm<int> Tpetra_Comm;
1036  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1037  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1038  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1039  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1040 
1041  // Ensure device is initialized
1042  if ( !Kokkos::is_initialized() )
1043  Kokkos::initialize();
1044 
1045  // 1-D Laplacian matrix
1046  GlobalOrdinal nrow = 50;
1047  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1048  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1049  RCP<const Tpetra_Map> map =
1050  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1051  nrow, comm);
1052  RCP<Tpetra_CrsGraph> graph =
1053  rcp(new Tpetra_CrsGraph(map, size_t(3)));
1054  Array<GlobalOrdinal> columnIndices(3);
1055  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1056  const size_t num_my_row = myGIDs.size();
1057  for (size_t i=0; i<num_my_row; ++i) {
1058  const GlobalOrdinal row = myGIDs[i];
1059  if (row == 0 || row == nrow-1) { // Boundary nodes
1060  columnIndices[0] = row;
1061  graph->insertGlobalIndices(row, columnIndices(0,1));
1062  }
1063  else { // Interior nodes
1064  columnIndices[0] = row-1;
1065  columnIndices[1] = row;
1066  columnIndices[2] = row+1;
1067  graph->insertGlobalIndices(row, columnIndices(0,3));
1068  }
1069  }
1070  graph->fillComplete();
1071  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1072 
1073  // Set values in matrix
1074  Array<Scalar> vals(3);
1075  Scalar a_val(VectorSize, BaseScalar(0.0));
1076  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1077  a_val.fastAccessCoeff(j) =
1078  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1079  }
1080  for (size_t i=0; i<num_my_row; ++i) {
1081  const GlobalOrdinal row = myGIDs[i];
1082  if (row == 0 || row == nrow-1) { // Boundary nodes
1083  columnIndices[0] = row;
1084  vals[0] = Scalar(1.0);
1085  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1086  }
1087  else {
1088  columnIndices[0] = row-1;
1089  columnIndices[1] = row;
1090  columnIndices[2] = row+1;
1091  vals[0] = Scalar(-1.0) * a_val;
1092  vals[1] = Scalar(2.0) * a_val;
1093  vals[2] = Scalar(-1.0) * a_val;
1094  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1095  }
1096  }
1097  matrix->fillComplete();
1098 
1099  matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1101 
1102  // Fill RHS vector
1103  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1104  Scalar b_val;
1105  {
1106  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1107  b_val = Scalar(VectorSize, BaseScalar(0.0));
1108  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1109  b_val.fastAccessCoeff(j) =
1110  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1111  }
1112  for (size_t i=0; i<num_my_row; ++i) {
1113  const GlobalOrdinal row = myGIDs[i];
1114  if (row == 0 || row == nrow-1)
1115  b_view(i,0) = Scalar(0.0);
1116  else
1117  b_view(i,0) = -Scalar(b_val * h * h);
1118  }
1119  }
1120 
1121  // Solve
1122  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1123 #if KOKKOS_VERSION >= 40799
1124  typedef KokkosKernels::ArithTraits<BaseScalar> BST;
1125 #else
1126  typedef Kokkos::ArithTraits<BaseScalar> BST;
1127 #endif
1128  typedef typename BST::mag_type base_mag_type;
1129  typedef typename Tpetra_Vector::mag_type mag_type;
1130  base_mag_type btol = 1e-9;
1131  mag_type tol = btol;
1132  int max_its = 1000;
1133  bool solved = Stokhos::CG_Solve(*matrix, *x, *b, tol, max_its,
1134  out.getOStream().get());
1135  TEST_EQUALITY_CONST( solved, true );
1136 
1137  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1138  // Teuchos::VERB_EXTREME);
1139 
1140  // Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1)
1141  btol = 1000*btol;
1142  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1143  Scalar val(VectorSize, BaseScalar(0.0));
1144  for (size_t i=0; i<num_my_row; ++i) {
1145  const GlobalOrdinal row = myGIDs[i];
1146  BaseScalar xx = row * h;
1147  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1148  val.fastAccessCoeff(j) =
1149  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0));
1150  }
1151  TEST_EQUALITY( x_view(i,0).size(), VectorSize );
1152 
1153  // Set small values to zero
1154  Scalar v = x_view(i,0);
1155  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1156  if (BST::abs(v.coeff(j)) < btol)
1157  v.fastAccessCoeff(j) = BaseScalar(0.0);
1158  if (BST::abs(val.coeff(j)) < btol)
1159  val.fastAccessCoeff(j) = BaseScalar(0.0);
1160  }
1161 
1162  for (LocalOrdinal j=0; j<VectorSize; ++j)
1163  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), btol);
1164  }
1165 
1166 }
1167 
1168 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2)
1169 
1170 //
1171 // Test simple CG solve with MueLu preconditioning for a 1-D Laplacian matrix
1172 //
1174  Tpetra_CrsMatrix_MP, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
1175 {
1176  using Teuchos::RCP;
1177  using Teuchos::rcp;
1178  using Teuchos::ArrayView;
1179  using Teuchos::Array;
1180  using Teuchos::ArrayRCP;
1181  using Teuchos::ParameterList;
1182  using Teuchos::getParametersFromXmlFile;
1183 
1184  typedef typename Storage::value_type BaseScalar;
1186 
1187  typedef Teuchos::Comm<int> Tpetra_Comm;
1188  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1189  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1190  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1191  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1192 
1193  // Ensure device is initialized
1194  if ( !Kokkos::is_initialized() )
1195  Kokkos::initialize();
1196 
1197  // 1-D Laplacian matrix
1198  GlobalOrdinal nrow = 50;
1199  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1200  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1201  RCP<const Tpetra_Map> map =
1202  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1203  nrow, comm);
1204  RCP<Tpetra_CrsGraph> graph =
1205  rcp(new Tpetra_CrsGraph(map, size_t(3)));
1206  Array<GlobalOrdinal> columnIndices(3);
1207  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1208  const size_t num_my_row = myGIDs.size();
1209  for (size_t i=0; i<num_my_row; ++i) {
1210  const GlobalOrdinal row = myGIDs[i];
1211  if (row == 0 || row == nrow-1) { // Boundary nodes
1212  columnIndices[0] = row;
1213  graph->insertGlobalIndices(row, columnIndices(0,1));
1214  }
1215  else { // Interior nodes
1216  columnIndices[0] = row-1;
1217  columnIndices[1] = row;
1218  columnIndices[2] = row+1;
1219  graph->insertGlobalIndices(row, columnIndices(0,3));
1220  }
1221  }
1222  graph->fillComplete();
1223  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1224 
1225  // Set values in matrix
1226  Array<Scalar> vals(3);
1227  Scalar a_val(VectorSize, BaseScalar(0.0));
1228  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1229  a_val.fastAccessCoeff(j) =
1230  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1231  }
1232  for (size_t i=0; i<num_my_row; ++i) {
1233  const GlobalOrdinal row = myGIDs[i];
1234  if (row == 0 || row == nrow-1) { // Boundary nodes
1235  columnIndices[0] = row;
1236  vals[0] = Scalar(1.0);
1237  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1238  }
1239  else {
1240  columnIndices[0] = row-1;
1241  columnIndices[1] = row;
1242  columnIndices[2] = row+1;
1243  vals[0] = Scalar(-1.0) * a_val;
1244  vals[1] = Scalar(2.0) * a_val;
1245  vals[2] = Scalar(-1.0) * a_val;
1246  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1247  }
1248  }
1249  matrix->fillComplete();
1250 
1251  // Fill RHS vector
1252  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1253  Scalar b_val;
1254  {
1255  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1256  b_val = Scalar(VectorSize, BaseScalar(0.0));
1257  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1258  b_val.fastAccessCoeff(j) =
1259  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1260  }
1261  for (size_t i=0; i<num_my_row; ++i) {
1262  const GlobalOrdinal row = myGIDs[i];
1263  if (row == 0 || row == nrow-1)
1264  b_view(i,0) = Scalar(0.0);
1265  else
1266  b_view(i,0) = -Scalar(b_val * h * h);
1267  }
1268  }
1269 
1270  // Create preconditioner
1271  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1272  RCP<OP> matrix_op = matrix;
1273  RCP<ParameterList> muelu_params =
1274  getParametersFromXmlFile("muelu_cheby.xml");
1275  RCP<OP> M =
1276  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1277 
1278  // Solve
1279  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1280 #if KOKKOS_VERSION >= 40799
1281  typedef KokkosKernels::ArithTraits<BaseScalar> BST;
1282 #else
1283  typedef Kokkos::ArithTraits<BaseScalar> BST;
1284 #endif
1285  typedef typename BST::mag_type base_mag_type;
1286  typedef typename Tpetra_Vector::mag_type mag_type;
1287  base_mag_type btol = 1e-9;
1288  mag_type tol = btol;
1289  int max_its = 1000;
1290  bool solved = Stokhos::PCG_Solve(*matrix, *x, *b, *M, tol, max_its,
1291  out.getOStream().get());
1292  TEST_EQUALITY_CONST( solved, true );
1293 
1294  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1295  // Teuchos::VERB_EXTREME);
1296 
1297  // Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1)
1298  btol = 1000*btol;
1299  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1300  Scalar val(VectorSize, BaseScalar(0.0));
1301  for (size_t i=0; i<num_my_row; ++i) {
1302  const GlobalOrdinal row = myGIDs[i];
1303  BaseScalar xx = row * h;
1304  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1305  val.fastAccessCoeff(j) =
1306  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0));
1307  }
1308  TEST_EQUALITY( x_view(i,0).size(), VectorSize );
1309 
1310  // Set small values to zero
1311  Scalar v = x_view(i,0);
1312  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1313  if (BST::magnitude(v.coeff(j)) < btol)
1314  v.fastAccessCoeff(j) = BaseScalar(0.0);
1315  if (BST::magnitude(val.coeff(j)) < btol)
1316  val.fastAccessCoeff(j) = BaseScalar(0.0);
1317  }
1318 
1319  for (LocalOrdinal j=0; j<VectorSize; ++j)
1320  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), btol);
1321  }
1322 
1323 }
1324 
1325 #else
1326 
1328  Tpetra_CrsMatrix_MP, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node ) {}
1329 
1330 #endif
1331 
1332 #if defined(HAVE_STOKHOS_BELOS)
1333 
1334 //
1335 // Test Belos GMRES solve for a simple banded upper-triangular matrix
1336 //
1338  Tpetra_CrsMatrix_MP, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1339 {
1340  using Teuchos::RCP;
1341  using Teuchos::rcp;
1342  using Teuchos::ArrayView;
1343  using Teuchos::Array;
1344  using Teuchos::ArrayRCP;
1345  using Teuchos::ParameterList;
1346 
1347  typedef typename Storage::value_type BaseScalar;
1349 
1350  typedef Teuchos::Comm<int> Tpetra_Comm;
1351  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1352  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
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::is_initialized() )
1358  Kokkos::initialize();
1359 
1360  // Build banded matrix
1361  GlobalOrdinal nrow = 10;
1362  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1363  RCP<const Tpetra_Map> map =
1364  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1365  nrow, comm);
1366  RCP<Tpetra_CrsGraph> graph =
1367  rcp(new Tpetra_CrsGraph(map, size_t(2)));
1368  Array<GlobalOrdinal> columnIndices(2);
1369  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1370  const size_t num_my_row = myGIDs.size();
1371  for (size_t i=0; i<num_my_row; ++i) {
1372  const GlobalOrdinal row = myGIDs[i];
1373  columnIndices[0] = row;
1374  size_t ncol = 1;
1375  if (row != nrow-1) {
1376  columnIndices[1] = row+1;
1377  ncol = 2;
1378  }
1379  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1380  }
1381  graph->fillComplete();
1382  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1383 
1384  // Set values in matrix
1385  Array<Scalar> vals(2);
1386  Scalar val(VectorSize, BaseScalar(0.0));
1387  for (size_t i=0; i<num_my_row; ++i) {
1388  const GlobalOrdinal row = myGIDs[i];
1389  columnIndices[0] = row;
1390  for (LocalOrdinal j=0; j<VectorSize; ++j)
1391  val.fastAccessCoeff(j) = j+1;
1392  vals[0] = val;
1393  size_t ncol = 1;
1394 
1395  if (row != nrow-1) {
1396  columnIndices[1] = row+1;
1397  for (LocalOrdinal j=0; j<VectorSize; ++j)
1398  val.fastAccessCoeff(j) = j+1;
1399  vals[1] = val;
1400  ncol = 2;
1401  }
1402  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1403  }
1404  matrix->fillComplete();
1405 
1406  // Fill RHS vector
1407  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1408  {
1409  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1410  for (size_t i=0; i<num_my_row; ++i) {
1411  b_view(i,0) = Scalar(1.0);
1412  }
1413  }
1414 
1415  // Solve
1417 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1418  typedef BaseScalar BelosScalar;
1419 #else
1420  typedef Scalar BelosScalar;
1421 #endif
1422  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1423  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1424  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1425  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1426  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1427  RCP<ParameterList> belosParams = rcp(new ParameterList);
1428  typename ST::magnitudeType tol = 1e-9;
1429  belosParams->set("Flexible Gmres", false);
1430  belosParams->set("Num Blocks", 100);
1431  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1432  belosParams->set("Maximum Iterations", 100);
1433  belosParams->set("Verbosity", 33);
1434  belosParams->set("Output Style", 1);
1435  belosParams->set("Output Frequency", 1);
1436  belosParams->set("Output Stream", out.getOStream());
1437  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1438  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1439  problem->setProblem();
1440  Belos::ReturnType ret = solver->solve();
1441  TEST_EQUALITY_CONST( ret, Belos::Converged );
1442 
1443  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1444  // Teuchos::VERB_EXTREME);
1445 
1446  // Check -- Correct answer is:
1447  // [ 0, 0, ..., 0 ]
1448  // [ 1, 1/2, ..., 1/VectorSize ]
1449  // [ 0, 0, ..., 0 ]
1450  // [ 1, 1/2, ..., 1/VectorSize ]
1451  // ....
1452  tol = 1000*tol;
1453  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1454  for (size_t i=0; i<num_my_row; ++i) {
1455  const GlobalOrdinal row = myGIDs[i];
1456  if (row % 2) {
1457  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1458  val.fastAccessCoeff(j) = BaseScalar(1.0) / BaseScalar(j+1);
1459  }
1460  }
1461  else
1462  val = Scalar(VectorSize, BaseScalar(0.0));
1463  TEST_EQUALITY( x_view(i,0).size(), VectorSize );
1464 
1465  // Set small values to zero
1466  Scalar v = x_view(i,0);
1467  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1468  if (ST::magnitude(v.coeff(j)) < tol)
1469  v.fastAccessCoeff(j) = BaseScalar(0.0);
1470  }
1471 
1472  for (LocalOrdinal j=0; j<VectorSize; ++j)
1473  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1474  }
1475 }
1476 
1477 //
1478 // Test Belos GMRES solve for a simple lower-triangular matrix with lucky breakdown with DGKS orthogonalization
1479 //
1481  Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1482 {
1483  using Teuchos::RCP;
1484  using Teuchos::rcp;
1485  using Teuchos::tuple;
1486  using Teuchos::ArrayView;
1487  using Teuchos::Array;
1488  using Teuchos::ArrayRCP;
1489  using Teuchos::ParameterList;
1490 
1491  typedef typename Storage::value_type BaseScalar;
1493 
1494  typedef Teuchos::Comm<int> Tpetra_Comm;
1495  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1496  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1497  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1498  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1499 
1500  // Ensure device is initialized
1501  if ( !Kokkos::is_initialized() )
1502  Kokkos::initialize();
1503 
1504  // Build diagonal matrix
1505  GlobalOrdinal nrow = 20;
1506  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1507  RCP<const Tpetra_Map> map =
1508  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1509  nrow, comm);
1510  RCP<Tpetra_CrsGraph> graph =
1511  rcp(new Tpetra_CrsGraph(map, size_t(1)));
1512  Array<GlobalOrdinal> columnIndices(1);
1513  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1514  const size_t num_my_row = myGIDs.size();
1515  for (size_t i=0; i<num_my_row; ++i) {
1516  const GlobalOrdinal row = myGIDs[i];
1517  columnIndices[0] = row;
1518  size_t ncol = 1;
1519  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1520  }
1521  graph->fillComplete();
1522  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1523 
1524  Array<Scalar> vals(1);
1525  for (size_t i=0; i<num_my_row; ++i) {
1526  const GlobalOrdinal row = myGIDs[i];
1527  columnIndices[0] = row;
1528  vals[0] = Scalar(row+1);
1529  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1530  }
1531  matrix->fillComplete();
1532 
1533  // Fill RHS vector:
1534  // [ 0, 0, ..., 0, 0, 1]
1535  // [ 0, 0, ..., 0, 2, 2]
1536  // [ 0, 0, ..., 3, 3, 3]
1537  // [ 0, 0, ..., 4, 4, 4]
1538  // ...
1539 
1540  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1541  {
1542  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1543  for (size_t i=0; i<num_my_row; ++i) {
1544  const GlobalOrdinal row = myGIDs[i];
1545  b_view(i,0) = Scalar(0.0);
1546  for (LocalOrdinal j=0; j<VectorSize; ++j)
1547  if (int(j+2+row-VectorSize) > 0)
1548  b_view(i,0).fastAccessCoeff(j) = BaseScalar(row+1);
1549  }
1550  }
1551 
1552  // Solve
1554 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1555  typedef BaseScalar BelosScalar;
1556 #else
1557  typedef Scalar BelosScalar;
1558 #endif
1559  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1560  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1561  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1562  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1563  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1564  RCP<ParameterList> belosParams = rcp(new ParameterList);
1565  typename ST::magnitudeType tol = 1e-9;
1566  belosParams->set("Flexible Gmres", false);
1567  belosParams->set("Num Blocks", 100);
1568  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1569  belosParams->set("Maximum Iterations", 100);
1570  belosParams->set("Verbosity", 33);
1571  belosParams->set("Output Style", 1);
1572  belosParams->set("Output Frequency", 1);
1573  belosParams->set("Output Stream", out.getOStream());
1574  belosParams->set("Orthogonalization","DGKS");
1575  RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1576  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1577  problem->setProblem();
1578  Belos::ReturnType ret = solver->solve();
1579  TEST_EQUALITY_CONST( ret, Belos::Converged );
1580 
1581 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1582  int numItersExpected = nrow;
1583  int numIters = solver->getNumIters();
1584  out << "numIters = " << numIters << std::endl;
1585  TEST_EQUALITY( numIters, numItersExpected);
1586 
1587  // Get and print number of ensemble iterations
1588  std::vector<int> ensemble_iterations =
1589  static_cast<const Belos::StatusTestImpResNorm<BelosScalar, MV, OP> *>(solver->getResidualStatusTest())->getEnsembleIterations();
1590  out << "Ensemble iterations = ";
1591  for (auto ensemble_iteration : ensemble_iterations)
1592  out << ensemble_iteration << " ";
1593  out << std::endl;
1594 
1595  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1596  if (int(j+1+nrow-VectorSize) > 0) {
1597  TEST_EQUALITY(int(j+1+nrow-VectorSize), ensemble_iterations[j]);
1598  }
1599  else {
1600  TEST_EQUALITY(int(0), ensemble_iterations[j]);
1601  }
1602  }
1603 #endif
1604 
1605  // Check -- Correct answer is:
1606  // [ 0, 0, ..., 0, 0, 1]
1607  // [ 0, 0, ..., 0, 1, 1]
1608  // [ 0, 0, ..., 1, 1, 1]
1609  // [ 0, 0, ..., 1, 1, 1]
1610  // ...
1611  tol = 1000*tol;
1612  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1613  for (size_t i=0; i<num_my_row; ++i) {
1614  const GlobalOrdinal row = myGIDs[i];
1615  Scalar v = x_view(i,0);
1616 
1617  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1618  if (ST::magnitude(v.coeff(j)) < tol)
1619  v.fastAccessCoeff(j) = BaseScalar(0.0);
1620  }
1621 
1622  Scalar val = Scalar(0.0);
1623 
1624  for (LocalOrdinal j=0; j<VectorSize; ++j)
1625  if (j+2+row-VectorSize > 0)
1626  val.fastAccessCoeff(j) = BaseScalar(1.0);
1627 
1628  for (LocalOrdinal j=0; j<VectorSize; ++j)
1629  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1630  }
1631 }
1632 
1633 //
1634 // Test Belos GMRES solve for a simple lower-triangular matrix with lucky breakdown with ICGS orthogonalization
1635 //
1637  Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1638 {
1639  using Teuchos::RCP;
1640  using Teuchos::rcp;
1641  using Teuchos::tuple;
1642  using Teuchos::ArrayView;
1643  using Teuchos::Array;
1644  using Teuchos::ArrayRCP;
1645  using Teuchos::ParameterList;
1646 
1647  typedef typename Storage::value_type BaseScalar;
1649 
1650  typedef Teuchos::Comm<int> Tpetra_Comm;
1651  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1652  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1653  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1654  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1655 
1656  // Ensure device is initialized
1657  if ( !Kokkos::is_initialized() )
1658  Kokkos::initialize();
1659 
1660  // Build diagonal matrix
1661  GlobalOrdinal nrow = 20;
1662  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1663  RCP<const Tpetra_Map> map =
1664  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1665  nrow, comm);
1666  RCP<Tpetra_CrsGraph> graph =
1667  rcp(new Tpetra_CrsGraph(map, size_t(1)));
1668  Array<GlobalOrdinal> columnIndices(1);
1669  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1670  const size_t num_my_row = myGIDs.size();
1671  for (size_t i=0; i<num_my_row; ++i) {
1672  const GlobalOrdinal row = myGIDs[i];
1673  columnIndices[0] = row;
1674  size_t ncol = 1;
1675  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1676  }
1677  graph->fillComplete();
1678  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1679 
1680  Array<Scalar> vals(1);
1681  for (size_t i=0; i<num_my_row; ++i) {
1682  const GlobalOrdinal row = myGIDs[i];
1683  columnIndices[0] = row;
1684  vals[0] = Scalar(row+1);
1685  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1686  }
1687  matrix->fillComplete();
1688 
1689  // Fill RHS vector:
1690  // [ 0, 0, ..., 0, 0, 1]
1691  // [ 0, 0, ..., 0, 2, 2]
1692  // [ 0, 0, ..., 3, 3, 3]
1693  // [ 0, 0, ..., 4, 4, 4]
1694  // ...
1695 
1696  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1697  {
1698  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1699  for (size_t i=0; i<num_my_row; ++i) {
1700  const GlobalOrdinal row = myGIDs[i];
1701  b_view(i,0) = Scalar(0.0);
1702  for (LocalOrdinal j=0; j<VectorSize; ++j)
1703  if (int(j+2+row-VectorSize) > 0)
1704  b_view(i,0).fastAccessCoeff(j) = BaseScalar(row+1);
1705  }
1706  }
1707 
1708  // Solve
1710 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1711  typedef BaseScalar BelosScalar;
1712 #else
1713  typedef Scalar BelosScalar;
1714 #endif
1715  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1716  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1717  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1718  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1719  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1720  RCP<ParameterList> belosParams = rcp(new ParameterList);
1721  typename ST::magnitudeType tol = 1e-9;
1722  belosParams->set("Flexible Gmres", false);
1723  belosParams->set("Num Blocks", 100);
1724  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1725  belosParams->set("Maximum Iterations", 100);
1726  belosParams->set("Verbosity", 33);
1727  belosParams->set("Output Style", 1);
1728  belosParams->set("Output Frequency", 1);
1729  belosParams->set("Output Stream", out.getOStream());
1730  belosParams->set("Orthogonalization","ICGS");
1731  RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1732  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1733  problem->setProblem();
1734  Belos::ReturnType ret = solver->solve();
1735  TEST_EQUALITY_CONST( ret, Belos::Converged );
1736 
1737 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1738  int numItersExpected = nrow;
1739  int numIters = solver->getNumIters();
1740  out << "numIters = " << numIters << std::endl;
1741  TEST_EQUALITY( numIters, numItersExpected);
1742 
1743  // Get and print number of ensemble iterations
1744  std::vector<int> ensemble_iterations =
1745  static_cast<const Belos::StatusTestImpResNorm<BelosScalar, MV, OP> *>(solver->getResidualStatusTest())->getEnsembleIterations();
1746  out << "Ensemble iterations = ";
1747  for (auto ensemble_iteration : ensemble_iterations)
1748  out << ensemble_iteration << " ";
1749  out << std::endl;
1750 
1751  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1752  if (int(j+1+nrow-VectorSize) > 0) {
1753  TEST_EQUALITY(int(j+1+nrow-VectorSize), ensemble_iterations[j]);
1754  }
1755  else {
1756  TEST_EQUALITY(int(0), ensemble_iterations[j]);
1757  }
1758  }
1759 #endif
1760 
1761  // Check -- Correct answer is:
1762  // [ 0, 0, ..., 0, 0, 1]
1763  // [ 0, 0, ..., 0, 1, 1]
1764  // [ 0, 0, ..., 1, 1, 1]
1765  // [ 0, 0, ..., 1, 1, 1]
1766  // ...
1767  tol = 1000*tol;
1768  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1769  for (size_t i=0; i<num_my_row; ++i) {
1770  const GlobalOrdinal row = myGIDs[i];
1771  Scalar v = x_view(i,0);
1772 
1773  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1774  if (ST::magnitude(v.coeff(j)) < tol)
1775  v.fastAccessCoeff(j) = BaseScalar(0.0);
1776  }
1777 
1778  Scalar val = Scalar(0.0);
1779 
1780  for (LocalOrdinal j=0; j<VectorSize; ++j)
1781  if (j+2+row-VectorSize > 0)
1782  val.fastAccessCoeff(j) = BaseScalar(1.0);
1783 
1784  for (LocalOrdinal j=0; j<VectorSize; ++j)
1785  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1786  }
1787 }
1788 
1789 //
1790 // Test Belos GMRES solve for a simple lower-triangular matrix with lucky breakdown with IMGS orthogonalization
1791 //
1793  Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1794 {
1795  using Teuchos::RCP;
1796  using Teuchos::rcp;
1797  using Teuchos::tuple;
1798  using Teuchos::ArrayView;
1799  using Teuchos::Array;
1800  using Teuchos::ArrayRCP;
1801  using Teuchos::ParameterList;
1802 
1803  typedef typename Storage::value_type BaseScalar;
1805 
1806  typedef Teuchos::Comm<int> Tpetra_Comm;
1807  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1808  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1809  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1810  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1811 
1812  // Ensure device is initialized
1813  if ( !Kokkos::is_initialized() )
1814  Kokkos::initialize();
1815 
1816  // Build diagonal matrix
1817  GlobalOrdinal nrow = 20;
1818  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1819  RCP<const Tpetra_Map> map =
1820  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1821  nrow, comm);
1822  RCP<Tpetra_CrsGraph> graph =
1823  rcp(new Tpetra_CrsGraph(map, size_t(1)));
1824  Array<GlobalOrdinal> columnIndices(1);
1825  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1826  const size_t num_my_row = myGIDs.size();
1827  for (size_t i=0; i<num_my_row; ++i) {
1828  const GlobalOrdinal row = myGIDs[i];
1829  columnIndices[0] = row;
1830  size_t ncol = 1;
1831  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1832  }
1833  graph->fillComplete();
1834  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1835 
1836  Array<Scalar> vals(1);
1837  for (size_t i=0; i<num_my_row; ++i) {
1838  const GlobalOrdinal row = myGIDs[i];
1839  columnIndices[0] = row;
1840  vals[0] = Scalar(row+1);
1841  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1842  }
1843  matrix->fillComplete();
1844 
1845  // Fill RHS vector:
1846  // [ 0, 0, ..., 0, 0, 1]
1847  // [ 0, 0, ..., 0, 2, 2]
1848  // [ 0, 0, ..., 3, 3, 3]
1849  // [ 0, 0, ..., 4, 4, 4]
1850  // ...
1851 
1852  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1853  {
1854  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1855  for (size_t i=0; i<num_my_row; ++i) {
1856  const GlobalOrdinal row = myGIDs[i];
1857  b_view(i,0) = Scalar(0.0);
1858  for (LocalOrdinal j=0; j<VectorSize; ++j)
1859  if (int(j+2+row-VectorSize) > 0)
1860  b_view(i,0).fastAccessCoeff(j) = BaseScalar(row+1);
1861  }
1862  }
1863 
1864  // Solve
1866 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1867  typedef BaseScalar BelosScalar;
1868 #else
1869  typedef Scalar BelosScalar;
1870 #endif
1871  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1872  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1873  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1874  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1875  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1876  RCP<ParameterList> belosParams = rcp(new ParameterList);
1877  typename ST::magnitudeType tol = 1e-9;
1878  belosParams->set("Flexible Gmres", false);
1879  belosParams->set("Num Blocks", 100);
1880  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1881  belosParams->set("Maximum Iterations", 100);
1882  belosParams->set("Verbosity", 33);
1883  belosParams->set("Output Style", 1);
1884  belosParams->set("Output Frequency", 1);
1885  belosParams->set("Output Stream", out.getOStream());
1886  belosParams->set("Orthogonalization","IMGS");
1887  RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1888  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1889  problem->setProblem();
1890  Belos::ReturnType ret = solver->solve();
1891  TEST_EQUALITY_CONST( ret, Belos::Converged );
1892 
1893 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1894  int numItersExpected = nrow;
1895  int numIters = solver->getNumIters();
1896  out << "numIters = " << numIters << std::endl;
1897  TEST_EQUALITY( numIters, numItersExpected);
1898 
1899  // Get and print number of ensemble iterations
1900  std::vector<int> ensemble_iterations =
1901  static_cast<const Belos::StatusTestImpResNorm<BelosScalar, MV, OP> *>(solver->getResidualStatusTest())->getEnsembleIterations();
1902  out << "Ensemble iterations = ";
1903  for (auto ensemble_iteration : ensemble_iterations)
1904  out << ensemble_iteration << " ";
1905  out << std::endl;
1906 
1907  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1908  if (int(j+1+nrow-VectorSize) > 0) {
1909  TEST_EQUALITY(int(j+1+nrow-VectorSize), ensemble_iterations[j]);
1910  }
1911  else {
1912  TEST_EQUALITY(int(0), ensemble_iterations[j]);
1913  }
1914  }
1915 #endif
1916 
1917  // Check -- Correct answer is:
1918  // [ 0, 0, ..., 0, 0, 1]
1919  // [ 0, 0, ..., 0, 1, 1]
1920  // [ 0, 0, ..., 1, 1, 1]
1921  // [ 0, 0, ..., 1, 1, 1]
1922  // ...
1923  tol = 1000*tol;
1924  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1925  for (size_t i=0; i<num_my_row; ++i) {
1926  const GlobalOrdinal row = myGIDs[i];
1927  Scalar v = x_view(i,0);
1928 
1929  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1930  if (ST::magnitude(v.coeff(j)) < tol)
1931  v.fastAccessCoeff(j) = BaseScalar(0.0);
1932  }
1933 
1934  Scalar val = Scalar(0.0);
1935 
1936  for (LocalOrdinal j=0; j<VectorSize; ++j)
1937  if (j+2+row-VectorSize > 0)
1938  val.fastAccessCoeff(j) = BaseScalar(1.0);
1939 
1940  for (LocalOrdinal j=0; j<VectorSize; ++j)
1941  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1942  }
1943 }
1944 
1945 #else
1946 
1948  Tpetra_CrsMatrix_MP, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1949 {}
1950 
1952  Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1953 {}
1954 
1956  Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1957 {}
1958 
1960  Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1961 {}
1962 
1963 #endif
1964 
1965 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
1966 
1967 //
1968 // Test Belos GMRES solve with Ifpack2 RILUK preconditioning for a
1969 // simple banded upper-triangular matrix
1970 //
1972  Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
1973 {
1974  using Teuchos::RCP;
1975  using Teuchos::rcp;
1976  using Teuchos::ArrayView;
1977  using Teuchos::Array;
1978  using Teuchos::ArrayRCP;
1979  using Teuchos::ParameterList;
1980 
1981  typedef typename Storage::value_type BaseScalar;
1983 
1984  typedef Teuchos::Comm<int> Tpetra_Comm;
1985  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1986  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1987  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1988  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1989 
1990  // Ensure device is initialized
1991  if ( !Kokkos::is_initialized() )
1992  Kokkos::initialize();
1993 
1994  // Build banded matrix
1995  GlobalOrdinal nrow = 10;
1996  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1997  RCP<const Tpetra_Map> map =
1998  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1999  nrow, comm);
2000  RCP<Tpetra_CrsGraph> graph =
2001  rcp(new Tpetra_CrsGraph(map, size_t(2)));
2002  Array<GlobalOrdinal> columnIndices(2);
2003  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2004  const size_t num_my_row = myGIDs.size();
2005  for (size_t i=0; i<num_my_row; ++i) {
2006  const GlobalOrdinal row = myGIDs[i];
2007  columnIndices[0] = row;
2008  size_t ncol = 1;
2009  if (row != nrow-1) {
2010  columnIndices[1] = row+1;
2011  ncol = 2;
2012  }
2013  graph->insertGlobalIndices(row, columnIndices(0,ncol));
2014  }
2015  graph->fillComplete();
2016  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
2017 
2018  // Set values in matrix
2019  Array<Scalar> vals(2);
2020  Scalar val(VectorSize, BaseScalar(0.0));
2021  for (size_t i=0; i<num_my_row; ++i) {
2022  const GlobalOrdinal row = myGIDs[i];
2023  columnIndices[0] = row;
2024  for (LocalOrdinal j=0; j<VectorSize; ++j)
2025  val.fastAccessCoeff(j) = j+1;
2026  vals[0] = val;
2027  size_t ncol = 1;
2028 
2029  if (row != nrow-1) {
2030  columnIndices[1] = row+1;
2031  for (LocalOrdinal j=0; j<VectorSize; ++j)
2032  val.fastAccessCoeff(j) = j+1;
2033  vals[1] = val;
2034  ncol = 2;
2035  }
2036  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
2037  }
2038  matrix->fillComplete();
2039 
2040  // Fill RHS vector
2041  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2042  {
2043  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
2044  for (size_t i=0; i<num_my_row; ++i) {
2045  b_view(i,0) = Scalar(1.0);
2046  }
2047  }
2048 
2049  // Create preconditioner
2050  typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
2051  Ifpack2::Factory factory;
2052  RCP<Prec> M = factory.create<Tpetra_CrsMatrix>("RILUK", matrix);
2053  M->initialize();
2054  M->compute();
2055 
2056  // Solve
2058 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2059  typedef BaseScalar BelosScalar;
2060 #else
2061  typedef Scalar BelosScalar;
2062 #endif
2063  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2064  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2065  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2066  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2067  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
2068  problem->setRightPrec(M);
2069  problem->setProblem();
2070  RCP<ParameterList> belosParams = rcp(new ParameterList);
2071  typename ST::magnitudeType tol = 1e-9;
2072  belosParams->set("Flexible Gmres", false);
2073  belosParams->set("Num Blocks", 100);
2074  belosParams->set("Convergence Tolerance", BelosScalar(tol));
2075  belosParams->set("Maximum Iterations", 100);
2076  belosParams->set("Verbosity", 33);
2077  belosParams->set("Output Style", 1);
2078  belosParams->set("Output Frequency", 1);
2079  belosParams->set("Output Stream", out.getOStream());
2080  //belosParams->set("Orthogonalization", "TSQR");
2081  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2082  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2083  Belos::ReturnType ret = solver->solve();
2084  TEST_EQUALITY_CONST( ret, Belos::Converged );
2085 
2086  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2087  // Teuchos::VERB_EXTREME);
2088 
2089  // Check -- Correct answer is:
2090  // [ 0, 0, ..., 0 ]
2091  // [ 1, 1/2, ..., 1/VectorSize ]
2092  // [ 0, 0, ..., 0 ]
2093  // [ 1, 1/2, ..., 1/VectorSize ]
2094  // ....
2095  tol = 1000*tol;
2096  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
2097  for (size_t i=0; i<num_my_row; ++i) {
2098  const GlobalOrdinal row = myGIDs[i];
2099  if (row % 2) {
2100  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2101  val.fastAccessCoeff(j) = BaseScalar(1.0) / BaseScalar(j+1);
2102  }
2103  }
2104  else
2105  val = Scalar(VectorSize, BaseScalar(0.0));
2106  TEST_EQUALITY( x_view(i,0).size(), VectorSize );
2107 
2108  // Set small values to zero
2109  Scalar v = x_view(i,0);
2110  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2111  if (ST::magnitude(v.coeff(j)) < tol)
2112  v.fastAccessCoeff(j) = BaseScalar(0.0);
2113  }
2114 
2115  for (LocalOrdinal j=0; j<VectorSize; ++j)
2116  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
2117  }
2118 }
2119 
2120 #else
2121 
2123  Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
2124 {}
2125 
2126 #endif
2127 
2128 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2)
2129 
2130 //
2131 // Test Belos CG solve with MueLu preconditioning for a 1-D Laplacian matrix
2132 //
2134  Tpetra_CrsMatrix_MP, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
2135 {
2136  using Teuchos::RCP;
2137  using Teuchos::rcp;
2138  using Teuchos::ArrayView;
2139  using Teuchos::Array;
2140  using Teuchos::ArrayRCP;
2141  using Teuchos::ParameterList;
2142  using Teuchos::getParametersFromXmlFile;
2143 
2144  typedef typename Storage::value_type BaseScalar;
2146 
2147  typedef Teuchos::Comm<int> Tpetra_Comm;
2148  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2149  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2150  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2151  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2152 
2153  // Ensure device is initialized
2154  if ( !Kokkos::is_initialized() )
2155  Kokkos::initialize();
2156 
2157  // 1-D Laplacian matrix
2158  GlobalOrdinal nrow = 50;
2159  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
2160  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2161  RCP<const Tpetra_Map> map =
2162  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2163  nrow, comm);
2164  RCP<Tpetra_CrsGraph> graph =
2165  rcp(new Tpetra_CrsGraph(map, size_t(3)));
2166  Array<GlobalOrdinal> columnIndices(3);
2167  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2168  const size_t num_my_row = myGIDs.size();
2169  for (size_t i=0; i<num_my_row; ++i) {
2170  const GlobalOrdinal row = myGIDs[i];
2171  if (row == 0 || row == nrow-1) { // Boundary nodes
2172  columnIndices[0] = row;
2173  graph->insertGlobalIndices(row, columnIndices(0,1));
2174  }
2175  else { // Interior nodes
2176  columnIndices[0] = row-1;
2177  columnIndices[1] = row;
2178  columnIndices[2] = row+1;
2179  graph->insertGlobalIndices(row, columnIndices(0,3));
2180  }
2181  }
2182  graph->fillComplete();
2183  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
2184 
2185  // Set values in matrix
2186  Array<Scalar> vals(3);
2187  Scalar a_val(VectorSize, BaseScalar(0.0));
2188  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2189  a_val.fastAccessCoeff(j) =
2190  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize);
2191  }
2192  for (size_t i=0; i<num_my_row; ++i) {
2193  const GlobalOrdinal row = myGIDs[i];
2194  if (row == 0 || row == nrow-1) { // Boundary nodes
2195  columnIndices[0] = row;
2196  vals[0] = Scalar(1.0);
2197  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2198  }
2199  else {
2200  columnIndices[0] = row-1;
2201  columnIndices[1] = row;
2202  columnIndices[2] = row+1;
2203  vals[0] = Scalar(-1.0) * a_val;
2204  vals[1] = Scalar(2.0) * a_val;
2205  vals[2] = Scalar(-1.0) * a_val;
2206  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2207  }
2208  }
2209  matrix->fillComplete();
2210 
2211  // Fill RHS vector
2212  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2213  Scalar b_val;
2214  {
2215  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
2216  b_val = Scalar(VectorSize, BaseScalar(0.0));
2217  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2218  b_val.fastAccessCoeff(j) =
2219  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize);
2220  }
2221  for (size_t i=0; i<num_my_row; ++i) {
2222  const GlobalOrdinal row = myGIDs[i];
2223  if (row == 0 || row == nrow-1)
2224  b_view(i,0) = Scalar(0.0);
2225  else
2226  b_view(i,0) = -Scalar(b_val * h * h);
2227  }
2228  }
2229 
2230  // Create preconditioner
2231  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2232  RCP<ParameterList> muelu_params =
2233  getParametersFromXmlFile("muelu_cheby.xml");
2234  RCP<OP> matrix_op = matrix;
2235  RCP<OP> M =
2236  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
2237 
2238  // Solve
2240 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2241  typedef BaseScalar BelosScalar;
2242 #else
2243  typedef Scalar BelosScalar;
2244 #endif
2245  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2246  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2247  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2248  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
2249  problem->setRightPrec(M);
2250  problem->setProblem();
2251  RCP<ParameterList> belosParams = rcp(new ParameterList);
2252  typename ST::magnitudeType tol = 1e-9;
2253  belosParams->set("Num Blocks", 100);
2254  belosParams->set("Convergence Tolerance", BelosScalar(tol));
2255  belosParams->set("Maximum Iterations", 100);
2256  belosParams->set("Verbosity", 33);
2257  belosParams->set("Output Style", 1);
2258  belosParams->set("Output Frequency", 1);
2259  belosParams->set("Output Stream", out.getOStream());
2260  // Turn off residual scaling so we can see some variation in the number
2261  // of iterations across the ensemble when not doing ensemble reductions
2262  belosParams->set("Implicit Residual Scaling", "None");
2263 
2264  RCP<Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true> > solver =
2265  rcp(new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true>(problem, belosParams));
2266  Belos::ReturnType ret = solver->solve();
2267  TEST_EQUALITY_CONST( ret, Belos::Converged );
2268 
2269 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
2270  // Get and print number of ensemble iterations
2271  std::vector<int> ensemble_iterations =
2272  solver->getResidualStatusTest()->getEnsembleIterations();
2273  out << "Ensemble iterations = ";
2274  for (int i=0; i<VectorSize; ++i)
2275  out << ensemble_iterations[i] << " ";
2276  out << std::endl;
2277 #endif
2278 
2279  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2280  // Teuchos::VERB_EXTREME);
2281 
2282  // Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1)
2283  tol = 1000*tol;
2284  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
2285  Scalar val(VectorSize, BaseScalar(0.0));
2286  for (size_t i=0; i<num_my_row; ++i) {
2287  const GlobalOrdinal row = myGIDs[i];
2288  BaseScalar xx = row * h;
2289  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2290  val.fastAccessCoeff(j) =
2291  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0));
2292  }
2293  TEST_EQUALITY( x_view(i,0).size(), VectorSize );
2294 
2295  // Set small values to zero
2296  Scalar v = x_view(i,0);
2297  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2298  if (ST::magnitude(v.coeff(j)) < tol)
2299  v.fastAccessCoeff(j) = BaseScalar(0.0);
2300  if (ST::magnitude(val.coeff(j)) < tol)
2301  val.fastAccessCoeff(j) = BaseScalar(0.0);
2302  }
2303 
2304  for (LocalOrdinal j=0; j<VectorSize; ++j)
2305  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
2306  }
2307 
2308 }
2309 
2310 #else
2311 
2313  Tpetra_CrsMatrix_MP, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
2314 {}
2315 
2316 #endif
2317 
2318 #if defined(HAVE_STOKHOS_AMESOS2)
2319 
2320 //
2321 // Test Amesos2 solve for a 1-D Laplacian matrix
2322 //
2323 #define TEST_AMESOS2_SOLVER(SolverName) \
2324  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_##SolverName, Storage, LocalOrdinal, GlobalOrdinal, Node) \
2325 { \
2326  using Teuchos::RCP; \
2327  using Teuchos::rcp; \
2328  using Teuchos::ArrayView; \
2329  using Teuchos::Array; \
2330  using Teuchos::ArrayRCP; \
2331  using Teuchos::ParameterList; \
2332  \
2333  typedef typename Storage::value_type BaseScalar; \
2334  typedef Sacado::MP::Vector<Storage> Scalar; \
2335  \
2336  typedef Teuchos::Comm<int> Tpetra_Comm; \
2337  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map; \
2338  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector; \
2339  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector; \
2340  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix; \
2341  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph; \
2342  \
2343  /* Ensure device is initialized */ \
2344  if ( !Kokkos::is_initialized() ) \
2345  Kokkos::initialize(); \
2346  \
2347  /* 1-D Laplacian matrix */ \
2348  GlobalOrdinal nrow = 50; \
2349  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1); \
2350  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm(); \
2351  RCP<const Tpetra_Map> map = \
2352  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>( \
2353  nrow, comm); \
2354  RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map, size_t(3)); \
2355  Array<GlobalOrdinal> columnIndices(3); \
2356  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList(); \
2357  const size_t num_my_row = myGIDs.size(); \
2358  for (size_t i=0; i<num_my_row; ++i) { \
2359  const GlobalOrdinal row = myGIDs[i]; \
2360  if (row == 0 || row == nrow-1) { /* Boundary nodes */ \
2361  columnIndices[0] = row; \
2362  graph->insertGlobalIndices(row, columnIndices(0,1)); \
2363  } \
2364  else { /* Interior nodes */ \
2365  columnIndices[0] = row-1; \
2366  columnIndices[1] = row; \
2367  columnIndices[2] = row+1; \
2368  graph->insertGlobalIndices(row, columnIndices(0,3)); \
2369  } \
2370  } \
2371  graph->fillComplete(); \
2372  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph)); \
2373  \
2374  /* Set values in matrix */ \
2375  Array<Scalar> vals(3); \
2376  Scalar a_val(VectorSize, BaseScalar(0.0)); \
2377  for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2378  a_val.fastAccessCoeff(j) = \
2379  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize); \
2380  } \
2381  for (size_t i=0; i<num_my_row; ++i) { \
2382  const GlobalOrdinal row = myGIDs[i]; \
2383  if (row == 0 || row == nrow-1) { /* Boundary nodes */ \
2384  columnIndices[0] = row; \
2385  vals[0] = Scalar(1.0); \
2386  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1)); \
2387  } \
2388  else { \
2389  columnIndices[0] = row-1; \
2390  columnIndices[1] = row; \
2391  columnIndices[2] = row+1; \
2392  vals[0] = Scalar(-1.0) * a_val; \
2393  vals[1] = Scalar(2.0) * a_val; \
2394  vals[2] = Scalar(-1.0) * a_val; \
2395  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3)); \
2396  } \
2397  } \
2398  matrix->fillComplete();\
2399  \
2400  /* Fill RHS vector */ \
2401  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map); \
2402  Scalar b_val; \
2403  { \
2404  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll); \
2405  b_val = Scalar(VectorSize, BaseScalar(0.0)); \
2406  for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2407  b_val.fastAccessCoeff(j) = \
2408  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize); \
2409  } \
2410  for (size_t i=0; i<num_my_row; ++i) { \
2411  const GlobalOrdinal row = myGIDs[i]; \
2412  if (row == 0 || row == nrow-1) \
2413  b_view(i,0) = Scalar(0.0); \
2414  else \
2415  b_view(i,0) = -Scalar(b_val * h * h); \
2416  } \
2417  } \
2418  \
2419  /* Solve */ \
2420  typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver; \
2421  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map); \
2422  std::string solver_name(#SolverName); \
2423  out << "Solving linear system with " << solver_name << std::endl; \
2424  RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>( \
2425  solver_name, matrix, x, b); \
2426  solver->solve(); \
2427  \
2428  /* Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1) */ \
2429  solver = Teuchos::null; /* Delete solver to eliminate live device views of x */ \
2430  typedef Teuchos::ScalarTraits<BaseScalar> ST; \
2431  typename ST::magnitudeType tol = 1e-9; \
2432  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly); \
2433  Scalar val(VectorSize, BaseScalar(0.0)); \
2434  for (size_t i=0; i<num_my_row; ++i) { \
2435  const GlobalOrdinal row = myGIDs[i]; \
2436  BaseScalar xx = row * h; \
2437  for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2438  val.fastAccessCoeff(j) = \
2439  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0)); \
2440  } \
2441  TEST_EQUALITY( x_view(i,0).size(), VectorSize ); \
2442  \
2443  /* Set small values to zero */ \
2444  Scalar v = x_view(i,0); \
2445  for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2446  if (ST::magnitude(v.coeff(j)) < tol) \
2447  v.fastAccessCoeff(j) = BaseScalar(0.0); \
2448  if (ST::magnitude(val.coeff(j)) < tol) \
2449  val.fastAccessCoeff(j) = BaseScalar(0.0); \
2450  } \
2451  \
2452  for (LocalOrdinal j=0; j<VectorSize; ++j) \
2453  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol); \
2454  } \
2455 }
2456 
2457 #endif // defined(HAVE_STOKHOS_AMESOS2)
2458 
2459 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_BASKER)
2460  TEST_AMESOS2_SOLVER(basker)
2461 #else
2462  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_basker, Storage, LocalOrdinal, GlobalOrdinal, Node) {}
2463 #endif
2464 
2465 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_KLU2)
2466  TEST_AMESOS2_SOLVER(klu2)
2467 #else
2469 #endif
2470 
2471 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_SUPERLUDIST)
2472  TEST_AMESOS2_SOLVER(superlu_dist)
2473 #else
2474  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_superlu_dist, Storage, LocalOrdinal, GlobalOrdinal, Node) {}
2475 #endif
2476 
2477 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_SUPERLUMT)
2478  TEST_AMESOS2_SOLVER(superlu_mt)
2479 #else
2480  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_superlu_mt, Storage, LocalOrdinal, GlobalOrdinal, Node) {}
2481 #endif
2482 
2483 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_SUPERLU)
2484  TEST_AMESOS2_SOLVER(superlu)
2485 #else
2486  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_superlu, Storage, LocalOrdinal, GlobalOrdinal, Node) {}
2487 #endif
2488 
2489 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_PARDISO_MKL)
2490  TEST_AMESOS2_SOLVER(pardisomkl)
2491 #else
2492  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_pardisomkl, Storage, LocalOrdinal, GlobalOrdinal, Node) {}
2493 #endif
2494 
2495 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_LAPACK)
2496  TEST_AMESOS2_SOLVER(lapack)
2497 #else
2498  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_lapack, Storage, LocalOrdinal, GlobalOrdinal, Node) {}
2499 #endif
2500 
2501 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2502  TEST_AMESOS2_SOLVER(cholmod)
2503 #else
2504  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_cholmod, Storage, LocalOrdinal, GlobalOrdinal, Node) {}
2505 #endif
2506 
2507 #define CRSMATRIX_MP_VECTOR_TESTS_SLGN(S, LO, GO, N) \
2508  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorAdd, S, LO, GO, N ) \
2509  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorDot, S, LO, GO, N ) \
2510  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorAdd, S, LO, GO, N ) \
2511  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDot, S, LO, GO, N ) \
2512  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDotSub, S, LO, GO, N ) \
2513  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixVectorMultiply, S, LO, GO, N ) \
2514  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2515  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, WrappedDualView, S, LO, GO, N ) \
2516  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Flatten, S, LO, GO, N ) \
2517  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimpleCG, S, LO, GO, N ) \
2518  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimplePCG_Muelu, S, LO, GO, N ) \
2519  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES, S, LO, GO, N ) \
2520  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, S, LO, GO, N ) \
2521  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, S, LO, GO, N ) \
2522  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, S, LO, GO, N ) \
2523  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, S, LO, GO, N ) \
2524  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosCG_Muelu, S, LO, GO, N ) \
2525  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_basker, S, LO, GO, N ) \
2526  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_klu2, S, LO, GO, N ) \
2527  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_superlu_dist, S, LO, GO, N ) \
2528  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_superlu_mt, S, LO, GO, N ) \
2529  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_superlu, S, LO, GO, N ) \
2530  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_pardisomkl, S, LO, GO, N ) \
2531  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_lapack, S, LO, GO, N ) \
2532  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_cholmod, S, LO, GO, N )
2533 
2534 #define CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) \
2535  typedef Stokhos::DeviceForNode<N>::type Device; \
2536  typedef Stokhos::StaticFixedStorage<int,double,VectorSize,Device::execution_space> SFS; \
2537  using default_global_ordinal_type = ::Tpetra::Map<>::global_ordinal_type; \
2538  using default_local_ordinal_type = ::Tpetra::Map<>::local_ordinal_type; \
2539  CRSMATRIX_MP_VECTOR_TESTS_SLGN(SFS, default_local_ordinal_type, default_global_ordinal_type, N)
2540 
2541 #define CRSMATRIX_MP_VECTOR_TESTS_N(N) \
2542  CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N)
2543 
2544 // Disabling testing of dynamic storage -- we don't really need it
2545  // typedef Stokhos::DynamicStorage<int,double,Device> DS;
2546  // using default_global_ordinal_type = ::Tpetra::Map<>::global_ordinal_type;
2547  // using default_local_ordinal_type = ::Tpetra::Map<>::local_ordinal_type;
2548  // CRSMATRIX_MP_VECTOR_TESTS_SLGN(DS, default_global_ordinal_type, default_local_ordinal_type, N)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &cijk_graph, const CijkType &cijk_dev)
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
scalar generate_multi_vector_coefficient(const ordinal nFEM, const ordinal nVec, const ordinal nStoch, const ordinal iColFEM, const ordinal iVec, const ordinal iStoch)
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
#define TEST_EQUALITY_CONST(v1, v2)
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_flat_mp_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_range_map, const LocalOrdinal block_size)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j)-expr2.val(j)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
expr val()
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType Node
BelosGMRES
#define TEST_EQUALITY(v1, v2)