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