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