Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_TpetraCrsMatrixMPVectorUnitTest.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stokhos Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include "Teuchos_UnitTestHelpers.hpp"
44 #include "Stokhos_Ensemble_Sizes.hpp"
45 
46 // Teuchos
47 #include "Teuchos_XMLParameterListCoreHelpers.hpp"
48 
49 // Tpetra
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
64 #include "BelosLinearProblem.hpp"
65 #include "BelosPseudoBlockGmresSolMgr.hpp"
66 #include "BelosPseudoBlockCGSolMgr.hpp"
67 #endif
68 
69 // Ifpack2 preconditioner
70 #ifdef HAVE_STOKHOS_IFPACK2
72 #include "Ifpack2_Factory.hpp"
73 #endif
74 
75 // MueLu preconditioner
76 #ifdef HAVE_STOKHOS_MUELU
78 #include "MueLu_CreateTpetraPreconditioner.hpp"
79 #endif
80 
81 // Amesos2 solver
82 #ifdef HAVE_STOKHOS_AMESOS2
84 #include "Amesos2_Factory.hpp"
85 #endif
86 
87 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 //
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;
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->getNodeElementList();
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  ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
160  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
161  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
162  for (size_t i=0; i<num_my_row; ++i) {
163  const GlobalOrdinal row = myGIDs[i];
164  for (LocalOrdinal j=0; j<VectorSize; ++j) {
165  val1.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
166  val2.fastAccessCoeff(j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
167  }
168  x1_view[i] = val1;
169  x2_view[i] = val2;
170  }
171  x1_view = Teuchos::null;
172  x2_view = Teuchos::null;
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  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
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].size(), VectorSize );
195  for (LocalOrdinal j=0; j<VectorSize; ++j)
196  TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(j), val.fastAccessCoeff(j), tol );
197  }
198 }
199 
200 //
201 // Test vector dot product
202 //
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;
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->getNodeElementList();
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  ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
239  ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
240  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
241  for (size_t i=0; i<num_my_row; ++i) {
242  const GlobalOrdinal row = myGIDs[i];
243  for (LocalOrdinal j=0; j<VectorSize; ++j) {
244  val1.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
245  val2.fastAccessCoeff(j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, VectorSize, row, j);
246  }
247  x1_view[i] = val1;
248  x2_view[i] = val2;
249  }
250  x1_view = Teuchos::null;
251  x2_view = Teuchos::null;
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 //
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;
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->getNodeElementList();
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  ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
341  ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
342  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
343  for (size_t i=0; i<num_my_row; ++i) {
344  const GlobalOrdinal row = myGIDs[i];
345  for (size_t j=0; j<ncol; ++j) {
346  for (LocalOrdinal k=0; k<VectorSize; ++k) {
347  BaseScalar v =
348  generate_multi_vector_coefficient<BaseScalar,size_t>(
349  nrow, ncol, VectorSize, row, j, k);
350  val1.fastAccessCoeff(k) = v;
351  val2.fastAccessCoeff(k) = 0.12345 * v;
352  }
353  x1_view[j][i] = val1;
354  x2_view[j][i] = val2;
355  }
356  }
357  x1_view = Teuchos::null;
358  x2_view = Teuchos::null;
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  ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
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[j][i].size(), VectorSize );
382  for (LocalOrdinal k=0; k<VectorSize; ++k)
383  TEST_FLOATING_EQUALITY( y_view[j][i].fastAccessCoeff(k),
384  val.fastAccessCoeff(k), tol );
385  }
386  }
387 }
388 
389 //
390 // Test multi-vector dot product
391 //
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;
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->getNodeElementList();
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  ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
429  ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
430  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
431  for (size_t i=0; i<num_my_row; ++i) {
432  const GlobalOrdinal row = myGIDs[i];
433  for (size_t j=0; j<ncol; ++j) {
434  for (LocalOrdinal k=0; k<VectorSize; ++k) {
435  BaseScalar v =
436  generate_multi_vector_coefficient<BaseScalar,size_t>(
437  nrow, ncol, VectorSize, row, j, k);
438  val1.fastAccessCoeff(k) = v;
439  val2.fastAccessCoeff(k) = 0.12345 * v;
440  }
441  x1_view[j][i] = val1;
442  x2_view[j][i] = val2;
443  }
444  }
445  x1_view = Teuchos::null;
446  x2_view = Teuchos::null;
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 //
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;
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->getNodeElementList();
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  ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
544  ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
545  Scalar val1(VectorSize, BaseScalar(0.0)), val2(VectorSize, BaseScalar(0.0));
546  for (size_t i=0; i<num_my_row; ++i) {
547  const GlobalOrdinal row = myGIDs[i];
548  for (size_t j=0; j<ncol; ++j) {
549  for (LocalOrdinal k=0; k<VectorSize; ++k) {
550  BaseScalar v =
551  generate_multi_vector_coefficient<BaseScalar,size_t>(
552  nrow, ncol, VectorSize, row, j, k);
553  val1.fastAccessCoeff(k) = v;
554  val2.fastAccessCoeff(k) = 0.12345 * v;
555  }
556  x1_view[j][i] = val1;
557  x2_view[j][i] = val2;
558  }
559  }
560  x1_view = Teuchos::null;
561  x2_view = Teuchos::null;
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 //
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;
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), Tpetra::StaticProfile));
661  Array<GlobalOrdinal> columnIndices(2);
662  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
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  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
704  for (size_t i=0; i<num_my_row; ++i) {
705  const GlobalOrdinal row = myGIDs[i];
706  for (LocalOrdinal j=0; j<VectorSize; ++j)
707  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
708  nrow, VectorSize, row, j);
709  x_view[i] = val;
710  }
711  x_view = Teuchos::null;
712 
713  // matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
714  // Teuchos::VERB_EXTREME);
715 
716  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
717  // Teuchos::VERB_EXTREME);
718 
719  // Multiply
720  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
721  matrix->apply(*x, *y);
722 
723  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
724  // Teuchos::VERB_EXTREME);
725 
726  // Check
727  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
728  BaseScalar tol = 1.0e-14;
729  for (size_t i=0; i<num_my_row; ++i) {
730  const GlobalOrdinal row = myGIDs[i];
731  for (LocalOrdinal j=0; j<VectorSize; ++j) {
732  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
733  nrow, VectorSize, row, j);
734  val.fastAccessCoeff(j) = v*v;
735  }
736  if (row != nrow-1) {
737  for (LocalOrdinal j=0; j<VectorSize; ++j) {
738  BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
739  nrow, VectorSize, row+1, j);
740  val.fastAccessCoeff(j) += v*v;
741  }
742  }
743  TEST_EQUALITY( y_view[i].size(), VectorSize );
744  for (LocalOrdinal j=0; j<VectorSize; ++j)
745  TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(j), val.fastAccessCoeff(j), tol );
746  }
747 }
748 
749 //
750 // Test matrix-multi-vector multiplication for a simple banded upper-triangular matrix
751 //
753  Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, Storage, LocalOrdinal, GlobalOrdinal, Node )
754 {
755  using Teuchos::RCP;
756  using Teuchos::rcp;
757  using Teuchos::ArrayView;
758  using Teuchos::Array;
759  using Teuchos::ArrayRCP;
760 
761  typedef typename Storage::value_type BaseScalar;
763 
764  typedef Teuchos::Comm<int> Tpetra_Comm;
765  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
766  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
767  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
768  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
769 
770  // Ensure device is initialized
771  if ( !Kokkos::is_initialized() )
772  Kokkos::initialize();
773 
774  // Build banded matrix
775  GlobalOrdinal nrow = 10;
776  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
777  RCP<const Tpetra_Map> map =
778  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
779  nrow, comm);
780  RCP<Tpetra_CrsGraph> graph =
781  rcp(new Tpetra_CrsGraph(map, size_t(2), Tpetra::StaticProfile));
782  Array<GlobalOrdinal> columnIndices(2);
783  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
784  const size_t num_my_row = myGIDs.size();
785  for (size_t i=0; i<num_my_row; ++i) {
786  const GlobalOrdinal row = myGIDs[i];
787  columnIndices[0] = row;
788  size_t ncol = 1;
789  if (row != nrow-1) {
790  columnIndices[1] = row+1;
791  ncol = 2;
792  }
793  graph->insertGlobalIndices(row, columnIndices(0,ncol));
794  }
795  graph->fillComplete();
796  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
797 
798  // Set values in matrix
799  Array<Scalar> vals(2);
800  Scalar val(VectorSize, BaseScalar(0.0));
801  for (size_t i=0; i<num_my_row; ++i) {
802  const GlobalOrdinal row = myGIDs[i];
803  columnIndices[0] = row;
804  for (LocalOrdinal j=0; j<VectorSize; ++j)
805  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
806  nrow, VectorSize, row, j);
807  vals[0] = val;
808  size_t ncol = 1;
809 
810  if (row != nrow-1) {
811  columnIndices[1] = row+1;
812  for (LocalOrdinal j=0; j<VectorSize; ++j)
813  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
814  nrow, VectorSize, row+1, j);
815  vals[1] = val;
816  ncol = 2;
817  }
818  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
819  }
820  matrix->fillComplete();
821 
822  // Fill multi-vector
823  size_t ncol = 5;
824  RCP<Tpetra_MultiVector> x = Tpetra::createMultiVector<Scalar>(map, ncol);
825  ArrayRCP< ArrayRCP<Scalar> > x_view = x->get2dViewNonConst();
826  for (size_t i=0; i<num_my_row; ++i) {
827  const GlobalOrdinal row = myGIDs[i];
828  for (size_t j=0; j<ncol; ++j) {
829  for (LocalOrdinal k=0; k<VectorSize; ++k) {
830  BaseScalar v =
831  generate_multi_vector_coefficient<BaseScalar,size_t>(
832  nrow, ncol, VectorSize, row, j, k);
833  val.fastAccessCoeff(k) = v;
834  }
835  x_view[j][i] = val;
836  }
837  }
838  x_view = Teuchos::null;
839 
840  // matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
841  // Teuchos::VERB_EXTREME);
842 
843  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
844  // Teuchos::VERB_EXTREME);
845 
846  // Multiply
847  RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
848  matrix->apply(*x, *y);
849 
850  // y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
851  // Teuchos::VERB_EXTREME);
852 
853  // Check
854  ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
855  BaseScalar tol = 1.0e-14;
856  for (size_t i=0; i<num_my_row; ++i) {
857  const GlobalOrdinal row = myGIDs[i];
858  for (size_t j=0; j<ncol; ++j) {
859  for (LocalOrdinal k=0; k<VectorSize; ++k) {
860  BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
861  nrow, VectorSize, row, k);
862  BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
863  nrow, ncol, VectorSize, row, j, k);
864  val.fastAccessCoeff(k) = v1*v2;
865  }
866  if (row != nrow-1) {
867  for (LocalOrdinal k=0; k<VectorSize; ++k) {
868  BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
869  nrow, VectorSize, row+1, k);
870  BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
871  nrow, ncol, VectorSize, row+1, j, k);
872  val.fastAccessCoeff(k) += v1*v2;
873  }
874  }
875  TEST_EQUALITY( y_view[j][i].size(), VectorSize );
876  for (LocalOrdinal k=0; k<VectorSize; ++k)
877  TEST_FLOATING_EQUALITY( y_view[j][i].fastAccessCoeff(k),
878  val.fastAccessCoeff(k), tol );
879  }
880  }
881 }
882 
883 //
884 // Test flattening MP::Vector matrix
885 //
887  Tpetra_CrsMatrix_MP, Flatten, Storage, LocalOrdinal, GlobalOrdinal, Node )
888 {
889  using Teuchos::RCP;
890  using Teuchos::rcp;
891  using Teuchos::ArrayView;
892  using Teuchos::Array;
893  using Teuchos::ArrayRCP;
894 
895  typedef typename Storage::value_type BaseScalar;
897 
898  typedef Teuchos::Comm<int> Tpetra_Comm;
899  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
900  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
901  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
902  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
903 
904  typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
905  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
906 
907  // Ensure device is initialized
908  if ( !Kokkos::is_initialized() )
909  Kokkos::initialize();
910 
911  // Build banded matrix
912  GlobalOrdinal nrow = 10;
913  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
914  RCP<const Tpetra_Map> map =
915  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
916  nrow, comm);
917  RCP<Tpetra_CrsGraph> graph =
918  rcp(new Tpetra_CrsGraph(map, size_t(2), Tpetra::StaticProfile));
919  Array<GlobalOrdinal> columnIndices(2);
920  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
921  const size_t num_my_row = myGIDs.size();
922  for (size_t i=0; i<num_my_row; ++i) {
923  const GlobalOrdinal row = myGIDs[i];
924  columnIndices[0] = row;
925  size_t ncol = 1;
926  if (row != nrow-1) {
927  columnIndices[1] = row+1;
928  ncol = 2;
929  }
930  graph->insertGlobalIndices(row, columnIndices(0,ncol));
931  }
932  graph->fillComplete();
933  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
934 
935  // Set values in matrix
936  Array<Scalar> vals(2);
937  Scalar val(VectorSize, BaseScalar(0.0));
938  for (size_t i=0; i<num_my_row; ++i) {
939  const GlobalOrdinal row = myGIDs[i];
940  columnIndices[0] = row;
941  for (LocalOrdinal j=0; j<VectorSize; ++j)
942  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
943  nrow, VectorSize, row, j);
944  vals[0] = val;
945  size_t ncol = 1;
946 
947  if (row != nrow-1) {
948  columnIndices[1] = row+1;
949  for (LocalOrdinal j=0; j<VectorSize; ++j)
950  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
951  nrow, VectorSize, row+1, j);
952  vals[1] = val;
953  ncol = 2;
954  }
955  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
956  }
957  matrix->fillComplete();
958 
959  // Fill vector
960  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
961  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
962  for (size_t i=0; i<num_my_row; ++i) {
963  const GlobalOrdinal row = myGIDs[i];
964  for (LocalOrdinal j=0; j<VectorSize; ++j)
965  val.fastAccessCoeff(j) = generate_vector_coefficient<BaseScalar,size_t>(
966  nrow, VectorSize, row, j);
967  x_view[i] = val;
968  }
969 
970  // Multiply
971  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
972  matrix->apply(*x, *y);
973 
974  // Flatten matrix
975  RCP<const Tpetra_Map> flat_x_map, flat_y_map;
976  RCP<const Tpetra_CrsGraph> flat_graph =
977  Stokhos::create_flat_mp_graph(*graph, flat_x_map, flat_y_map, VectorSize);
978  RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
979  Stokhos::create_flat_matrix(*matrix, flat_graph, VectorSize);
980 
981  // Multiply with flattened matix
982  RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
983  RCP<Flat_Tpetra_Vector> flat_x =
984  Stokhos::create_flat_vector_view(*x, flat_x_map);
985  RCP<Flat_Tpetra_Vector> flat_y =
986  Stokhos::create_flat_vector_view(*y2, flat_y_map);
987  flat_matrix->apply(*flat_x, *flat_y);
988 
989  // flat_y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
990  // Teuchos::VERB_EXTREME);
991 
992  // Check
993  BaseScalar tol = 1.0e-14;
994  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
995  ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
996  for (size_t i=0; i<num_my_row; ++i) {
997  TEST_EQUALITY( y_view[i].size(), VectorSize );
998  TEST_EQUALITY( y2_view[i].size(), VectorSize );
999  for (LocalOrdinal j=0; j<VectorSize; ++j)
1000  TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(j),
1001  y2_view[i].fastAccessCoeff(j), tol );
1002  }
1003 }
1004 
1005 //
1006 // Test interaction between Tpetra WrappedDualView and MP::Vector
1007 //
1009  Tpetra_CrsMatrix_MP, WrappedDualView, Storage, LocalOrdinal, GlobalOrdinal, Node )
1010 {
1011  //BMK 6-2021: This test is required because a View of MP::Vector has slightly different behavior than a typical Kokkos::View.
1012  //If you construct a Kokkos::View with a label and 0 extent, it gets a non-null allocation.
1013  //But for View<MP::Vector>, the same constructor produces a null data pointer but
1014  //an active reference counting node (use_count() > 0).
1015  //This test makes sure that Tpetra WrappedDualView works correctly with a View where data() == nullptr but use_count() > 0.
1016  using Teuchos::RCP;
1017  using Teuchos::rcp;
1018  using Teuchos::ArrayView;
1019  using Teuchos::Array;
1020  using Teuchos::ArrayRCP;
1021 
1022  typedef typename Storage::value_type BaseScalar;
1024 
1025  using DualViewType = Kokkos::DualView<Scalar*, typename Node::device_type>;
1026  using WDV = Tpetra::Details::WrappedDualView<DualViewType>;
1027  using values_view = typename DualViewType::t_dev;
1028 
1029  // Ensure device is initialized
1030  if ( !Kokkos::is_initialized() )
1031  Kokkos::initialize();
1032 
1033  WDV wdv;
1034  {
1035  values_view myView("emptyTestView", 0);
1036  wdv = WDV(myView);
1037  }
1038  size_t use_h = wdv.getHostView(Tpetra::Access::ReadOnly).use_count();
1039  size_t use_d = wdv.getDeviceView(Tpetra::Access::ReadOnly).use_count();
1040  //The WrappedDualView is now the only object holding references to the host and device views,
1041  //so they should have identical use counts.
1042  TEST_EQUALITY(use_h, use_d);
1043 }
1044 
1045 //
1046 // Test simple CG solve without preconditioning for a 1-D Laplacian matrix
1047 //
1049  Tpetra_CrsMatrix_MP, SimpleCG, Storage, LocalOrdinal, GlobalOrdinal, Node )
1050 {
1051  using Teuchos::RCP;
1052  using Teuchos::rcp;
1053  using Teuchos::ArrayView;
1054  using Teuchos::Array;
1055  using Teuchos::ArrayRCP;
1056  using Teuchos::ParameterList;
1057 
1058  typedef typename Storage::value_type BaseScalar;
1060 
1061  typedef Teuchos::Comm<int> Tpetra_Comm;
1062  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1063  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1064  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1065  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1066 
1067  // Ensure device is initialized
1068  if ( !Kokkos::is_initialized() )
1069  Kokkos::initialize();
1070 
1071  // 1-D Laplacian matrix
1072  GlobalOrdinal nrow = 50;
1073  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1074  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1075  RCP<const Tpetra_Map> map =
1076  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1077  nrow, comm);
1078  RCP<Tpetra_CrsGraph> graph =
1079  rcp(new Tpetra_CrsGraph(map, size_t(3), Tpetra::StaticProfile));
1080  Array<GlobalOrdinal> columnIndices(3);
1081  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1082  const size_t num_my_row = myGIDs.size();
1083  for (size_t i=0; i<num_my_row; ++i) {
1084  const GlobalOrdinal row = myGIDs[i];
1085  if (row == 0 || row == nrow-1) { // Boundary nodes
1086  columnIndices[0] = row;
1087  graph->insertGlobalIndices(row, columnIndices(0,1));
1088  }
1089  else { // Interior nodes
1090  columnIndices[0] = row-1;
1091  columnIndices[1] = row;
1092  columnIndices[2] = row+1;
1093  graph->insertGlobalIndices(row, columnIndices(0,3));
1094  }
1095  }
1096  graph->fillComplete();
1097  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1098 
1099  // Set values in matrix
1100  Array<Scalar> vals(3);
1101  Scalar a_val(VectorSize, BaseScalar(0.0));
1102  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1103  a_val.fastAccessCoeff(j) =
1104  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1105  }
1106  for (size_t i=0; i<num_my_row; ++i) {
1107  const GlobalOrdinal row = myGIDs[i];
1108  if (row == 0 || row == nrow-1) { // Boundary nodes
1109  columnIndices[0] = row;
1110  vals[0] = Scalar(1.0);
1111  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1112  }
1113  else {
1114  columnIndices[0] = row-1;
1115  columnIndices[1] = row;
1116  columnIndices[2] = row+1;
1117  vals[0] = Scalar(-1.0) * a_val;
1118  vals[1] = Scalar(2.0) * a_val;
1119  vals[2] = Scalar(-1.0) * a_val;
1120  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1121  }
1122  }
1123  matrix->fillComplete();
1124 
1125  matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1126  Teuchos::VERB_EXTREME);
1127 
1128  // Fill RHS vector
1129  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1130  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1131  Scalar b_val(VectorSize, BaseScalar(0.0));
1132  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1133  b_val.fastAccessCoeff(j) =
1134  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1135  }
1136  for (size_t i=0; i<num_my_row; ++i) {
1137  const GlobalOrdinal row = myGIDs[i];
1138  if (row == 0 || row == nrow-1)
1139  b_view[i] = Scalar(0.0);
1140  else
1141  b_view[i] = -Scalar(b_val * h * h);
1142  }
1143 
1144  // Solve
1145  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1146  typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1147  typedef typename BST::mag_type base_mag_type;
1148  typedef typename Tpetra_Vector::mag_type mag_type;
1149  base_mag_type btol = 1e-9;
1150  mag_type tol = btol;
1151  int max_its = 1000;
1152  bool solved = Stokhos::CG_Solve(*matrix, *x, *b, tol, max_its,
1153  out.getOStream().get());
1154  TEST_EQUALITY_CONST( solved, true );
1155 
1156  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1157  // Teuchos::VERB_EXTREME);
1158 
1159  // Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1)
1160  btol = 1000*btol;
1161  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1162  Scalar val(VectorSize, BaseScalar(0.0));
1163  for (size_t i=0; i<num_my_row; ++i) {
1164  const GlobalOrdinal row = myGIDs[i];
1165  BaseScalar xx = row * h;
1166  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1167  val.fastAccessCoeff(j) =
1168  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0));
1169  }
1170  TEST_EQUALITY( x_view[i].size(), VectorSize );
1171 
1172  // Set small values to zero
1173  Scalar v = x_view[i];
1174  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1175  if (BST::abs(v.coeff(j)) < btol)
1176  v.fastAccessCoeff(j) = BaseScalar(0.0);
1177  if (BST::abs(val.coeff(j)) < btol)
1178  val.fastAccessCoeff(j) = BaseScalar(0.0);
1179  }
1180 
1181  for (LocalOrdinal j=0; j<VectorSize; ++j)
1182  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), btol);
1183  }
1184 
1185 }
1186 
1187 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2)
1188 
1189 //
1190 // Test simple CG solve with MueLu preconditioning for a 1-D Laplacian matrix
1191 //
1193  Tpetra_CrsMatrix_MP, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
1194 {
1195  using Teuchos::RCP;
1196  using Teuchos::rcp;
1197  using Teuchos::ArrayView;
1198  using Teuchos::Array;
1199  using Teuchos::ArrayRCP;
1200  using Teuchos::ParameterList;
1201  using Teuchos::getParametersFromXmlFile;
1202 
1203  typedef typename Storage::value_type BaseScalar;
1205 
1206  typedef Teuchos::Comm<int> Tpetra_Comm;
1207  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1208  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1209  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1210  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1211 
1212  // Ensure device is initialized
1213  if ( !Kokkos::is_initialized() )
1214  Kokkos::initialize();
1215 
1216  // 1-D Laplacian matrix
1217  GlobalOrdinal nrow = 50;
1218  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
1219  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1220  RCP<const Tpetra_Map> map =
1221  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1222  nrow, comm);
1223  RCP<Tpetra_CrsGraph> graph =
1224  rcp(new Tpetra_CrsGraph(map, size_t(3), Tpetra::StaticProfile));
1225  Array<GlobalOrdinal> columnIndices(3);
1226  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1227  const size_t num_my_row = myGIDs.size();
1228  for (size_t i=0; i<num_my_row; ++i) {
1229  const GlobalOrdinal row = myGIDs[i];
1230  if (row == 0 || row == nrow-1) { // Boundary nodes
1231  columnIndices[0] = row;
1232  graph->insertGlobalIndices(row, columnIndices(0,1));
1233  }
1234  else { // Interior nodes
1235  columnIndices[0] = row-1;
1236  columnIndices[1] = row;
1237  columnIndices[2] = row+1;
1238  graph->insertGlobalIndices(row, columnIndices(0,3));
1239  }
1240  }
1241  graph->fillComplete();
1242  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1243 
1244  // Set values in matrix
1245  Array<Scalar> vals(3);
1246  Scalar a_val(VectorSize, BaseScalar(0.0));
1247  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1248  a_val.fastAccessCoeff(j) =
1249  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1250  }
1251  for (size_t i=0; i<num_my_row; ++i) {
1252  const GlobalOrdinal row = myGIDs[i];
1253  if (row == 0 || row == nrow-1) { // Boundary nodes
1254  columnIndices[0] = row;
1255  vals[0] = Scalar(1.0);
1256  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1257  }
1258  else {
1259  columnIndices[0] = row-1;
1260  columnIndices[1] = row;
1261  columnIndices[2] = row+1;
1262  vals[0] = Scalar(-1.0) * a_val;
1263  vals[1] = Scalar(2.0) * a_val;
1264  vals[2] = Scalar(-1.0) * a_val;
1265  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1266  }
1267  }
1268  matrix->fillComplete();
1269 
1270  // Fill RHS vector
1271  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1272  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1273  Scalar b_val(VectorSize, BaseScalar(0.0));
1274  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1275  b_val.fastAccessCoeff(j) =
1276  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize);
1277  }
1278  for (size_t i=0; i<num_my_row; ++i) {
1279  const GlobalOrdinal row = myGIDs[i];
1280  if (row == 0 || row == nrow-1)
1281  b_view[i] = Scalar(0.0);
1282  else
1283  b_view[i] = -Scalar(b_val * h * h);
1284  }
1285 
1286  // Create preconditioner
1287  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1288  RCP<OP> matrix_op = matrix;
1289  RCP<ParameterList> muelu_params =
1290  getParametersFromXmlFile("muelu_cheby.xml");
1291  RCP<OP> M =
1292  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1293 
1294  // Solve
1295  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1296  typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1297  typedef typename BST::mag_type base_mag_type;
1298  typedef typename Tpetra_Vector::mag_type mag_type;
1299  base_mag_type btol = 1e-9;
1300  mag_type tol = btol;
1301  int max_its = 1000;
1302  bool solved = Stokhos::PCG_Solve(*matrix, *x, *b, *M, tol, max_its,
1303  out.getOStream().get());
1304  TEST_EQUALITY_CONST( solved, true );
1305 
1306  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1307  // Teuchos::VERB_EXTREME);
1308 
1309  // Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1)
1310  btol = 1000*btol;
1311  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1312  Scalar val(VectorSize, BaseScalar(0.0));
1313  for (size_t i=0; i<num_my_row; ++i) {
1314  const GlobalOrdinal row = myGIDs[i];
1315  BaseScalar xx = row * h;
1316  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1317  val.fastAccessCoeff(j) =
1318  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0));
1319  }
1320  TEST_EQUALITY( x_view[i].size(), VectorSize );
1321 
1322  // Set small values to zero
1323  Scalar v = x_view[i];
1324  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1325  if (BST::magnitude(v.coeff(j)) < btol)
1326  v.fastAccessCoeff(j) = BaseScalar(0.0);
1327  if (BST::magnitude(val.coeff(j)) < btol)
1328  val.fastAccessCoeff(j) = BaseScalar(0.0);
1329  }
1330 
1331  for (LocalOrdinal j=0; j<VectorSize; ++j)
1332  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), btol);
1333  }
1334 
1335 }
1336 
1337 #else
1338 
1340  Tpetra_CrsMatrix_MP, SimplePCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node ) {}
1341 
1342 #endif
1343 
1344 #if defined(HAVE_STOKHOS_BELOS)
1345 
1346 //
1347 // Test Belos GMRES solve for a simple banded upper-triangular matrix
1348 //
1350  Tpetra_CrsMatrix_MP, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1351 {
1352  using Teuchos::RCP;
1353  using Teuchos::rcp;
1354  using Teuchos::ArrayView;
1355  using Teuchos::Array;
1356  using Teuchos::ArrayRCP;
1357  using Teuchos::ParameterList;
1358 
1359  typedef typename Storage::value_type BaseScalar;
1361 
1362  typedef Teuchos::Comm<int> Tpetra_Comm;
1363  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1364  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1365  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1366  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1367 
1368  // Ensure device is initialized
1369  if ( !Kokkos::is_initialized() )
1370  Kokkos::initialize();
1371 
1372  // Build banded matrix
1373  GlobalOrdinal nrow = 10;
1374  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1375  RCP<const Tpetra_Map> map =
1376  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1377  nrow, comm);
1378  RCP<Tpetra_CrsGraph> graph =
1379  rcp(new Tpetra_CrsGraph(map, size_t(2), Tpetra::StaticProfile));
1380  Array<GlobalOrdinal> columnIndices(2);
1381  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1382  const size_t num_my_row = myGIDs.size();
1383  for (size_t i=0; i<num_my_row; ++i) {
1384  const GlobalOrdinal row = myGIDs[i];
1385  columnIndices[0] = row;
1386  size_t ncol = 1;
1387  if (row != nrow-1) {
1388  columnIndices[1] = row+1;
1389  ncol = 2;
1390  }
1391  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1392  }
1393  graph->fillComplete();
1394  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1395 
1396  // Set values in matrix
1397  Array<Scalar> vals(2);
1398  Scalar val(VectorSize, BaseScalar(0.0));
1399  for (size_t i=0; i<num_my_row; ++i) {
1400  const GlobalOrdinal row = myGIDs[i];
1401  columnIndices[0] = row;
1402  for (LocalOrdinal j=0; j<VectorSize; ++j)
1403  val.fastAccessCoeff(j) = j+1;
1404  vals[0] = val;
1405  size_t ncol = 1;
1406 
1407  if (row != nrow-1) {
1408  columnIndices[1] = row+1;
1409  for (LocalOrdinal j=0; j<VectorSize; ++j)
1410  val.fastAccessCoeff(j) = j+1;
1411  vals[1] = val;
1412  ncol = 2;
1413  }
1414  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1415  }
1416  matrix->fillComplete();
1417 
1418  // Fill RHS vector
1419  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1420  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1421  for (size_t i=0; i<num_my_row; ++i) {
1422  b_view[i] = Scalar(1.0);
1423  }
1424 
1425  // Solve
1426  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1427 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1428  typedef BaseScalar BelosScalar;
1429 #else
1430  typedef Scalar BelosScalar;
1431 #endif
1432  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1433  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1434  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1435  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1436  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1437  RCP<ParameterList> belosParams = rcp(new ParameterList);
1438  typename ST::magnitudeType tol = 1e-9;
1439  belosParams->set("Flexible Gmres", false);
1440  belosParams->set("Num Blocks", 100);
1441  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1442  belosParams->set("Maximum Iterations", 100);
1443  belosParams->set("Verbosity", 33);
1444  belosParams->set("Output Style", 1);
1445  belosParams->set("Output Frequency", 1);
1446  belosParams->set("Output Stream", out.getOStream());
1447  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1448  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1449  problem->setProblem();
1450  Belos::ReturnType ret = solver->solve();
1451  TEST_EQUALITY_CONST( ret, Belos::Converged );
1452 
1453  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
1454  // Teuchos::VERB_EXTREME);
1455 
1456  // Check -- Correct answer is:
1457  // [ 0, 0, ..., 0 ]
1458  // [ 1, 1/2, ..., 1/VectorSize ]
1459  // [ 0, 0, ..., 0 ]
1460  // [ 1, 1/2, ..., 1/VectorSize ]
1461  // ....
1462  tol = 1000*tol;
1463  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1464  for (size_t i=0; i<num_my_row; ++i) {
1465  const GlobalOrdinal row = myGIDs[i];
1466  if (row % 2) {
1467  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1468  val.fastAccessCoeff(j) = BaseScalar(1.0) / BaseScalar(j+1);
1469  }
1470  }
1471  else
1472  val = Scalar(VectorSize, BaseScalar(0.0));
1473  TEST_EQUALITY( x_view[i].size(), VectorSize );
1474 
1475  // Set small values to zero
1476  Scalar v = x_view[i];
1477  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1478  if (ST::magnitude(v.coeff(j)) < tol)
1479  v.fastAccessCoeff(j) = BaseScalar(0.0);
1480  }
1481 
1482  for (LocalOrdinal j=0; j<VectorSize; ++j)
1483  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1484  }
1485 }
1486 
1487 //
1488 // Test Belos GMRES solve for a simple lower-triangular matrix with lucky breakdown with DGKS orthogonalization
1489 //
1491  Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1492 {
1493  using Teuchos::RCP;
1494  using Teuchos::rcp;
1495  using Teuchos::tuple;
1496  using Teuchos::ArrayView;
1497  using Teuchos::Array;
1498  using Teuchos::ArrayRCP;
1499  using Teuchos::ParameterList;
1500 
1501  typedef typename Storage::value_type BaseScalar;
1503 
1504  typedef Teuchos::Comm<int> Tpetra_Comm;
1505  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1506  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1507  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1508  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1509 
1510  // Ensure device is initialized
1511  if ( !Kokkos::is_initialized() )
1512  Kokkos::initialize();
1513 
1514  // Build diagonal matrix
1515  GlobalOrdinal nrow = 20;
1516  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1517  RCP<const Tpetra_Map> map =
1518  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1519  nrow, comm);
1520  RCP<Tpetra_CrsGraph> graph =
1521  rcp(new Tpetra_CrsGraph(map, size_t(1), Tpetra::StaticProfile));
1522  Array<GlobalOrdinal> columnIndices(1);
1523  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1524  const size_t num_my_row = myGIDs.size();
1525  for (size_t i=0; i<num_my_row; ++i) {
1526  const GlobalOrdinal row = myGIDs[i];
1527  columnIndices[0] = row;
1528  size_t ncol = 1;
1529  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1530  }
1531  graph->fillComplete();
1532  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1533 
1534  Array<Scalar> vals(1);
1535  for (size_t i=0; i<num_my_row; ++i) {
1536  const GlobalOrdinal row = myGIDs[i];
1537  columnIndices[0] = row;
1538  vals[0] = Scalar(row+1);
1539  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1540  }
1541  matrix->fillComplete();
1542 
1543  // Fill RHS vector:
1544  // [ 0, 0, ..., 0, 0, 1]
1545  // [ 0, 0, ..., 0, 2, 2]
1546  // [ 0, 0, ..., 3, 3, 3]
1547  // [ 0, 0, ..., 4, 4, 4]
1548  // ...
1549 
1550  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1551  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1552  for (size_t i=0; i<num_my_row; ++i) {
1553  const GlobalOrdinal row = myGIDs[i];
1554  b_view[i] = Scalar(0.0);
1555  for (LocalOrdinal j=0; j<VectorSize; ++j)
1556  if (int(j+2+row-VectorSize) > 0)
1557  b_view[i].fastAccessCoeff(j) = BaseScalar(row+1);
1558  }
1559 
1560  // Solve
1561  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1562 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1563  typedef BaseScalar BelosScalar;
1564 #else
1565  typedef Scalar BelosScalar;
1566 #endif
1567  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1568  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1569  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1570  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1571  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1572  RCP<ParameterList> belosParams = rcp(new ParameterList);
1573  typename ST::magnitudeType tol = 1e-9;
1574  belosParams->set("Flexible Gmres", false);
1575  belosParams->set("Num Blocks", 100);
1576  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1577  belosParams->set("Maximum Iterations", 100);
1578  belosParams->set("Verbosity", 33);
1579  belosParams->set("Output Style", 1);
1580  belosParams->set("Output Frequency", 1);
1581  belosParams->set("Output Stream", out.getOStream());
1582  belosParams->set("Orthogonalization","DGKS");
1583  RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1584  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1585  problem->setProblem();
1586  Belos::ReturnType ret = solver->solve();
1587  TEST_EQUALITY_CONST( ret, Belos::Converged );
1588 
1589 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1590  int numItersExpected = nrow;
1591  int numIters = solver->getNumIters();
1592  out << "numIters = " << numIters << std::endl;
1593  TEST_EQUALITY( numIters, numItersExpected);
1594 
1595  // Get and print number of ensemble iterations
1596  std::vector<int> ensemble_iterations =
1597  static_cast<const Belos::StatusTestImpResNorm<BelosScalar, MV, OP> *>(solver->getResidualStatusTest())->getEnsembleIterations();
1598  out << "Ensemble iterations = ";
1599  for (auto ensemble_iteration : ensemble_iterations)
1600  out << ensemble_iteration << " ";
1601  out << std::endl;
1602 
1603  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1604  if (int(j+1+nrow-VectorSize) > 0) {
1605  TEST_EQUALITY(int(j+1+nrow-VectorSize), ensemble_iterations[j]);
1606  }
1607  else {
1608  TEST_EQUALITY(int(0), ensemble_iterations[j]);
1609  }
1610  }
1611 #endif
1612 
1613  // Check -- Correct answer is:
1614  // [ 0, 0, ..., 0, 0, 1]
1615  // [ 0, 0, ..., 0, 1, 1]
1616  // [ 0, 0, ..., 1, 1, 1]
1617  // [ 0, 0, ..., 1, 1, 1]
1618  // ...
1619  tol = 1000*tol;
1620  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1621  for (size_t i=0; i<num_my_row; ++i) {
1622  const GlobalOrdinal row = myGIDs[i];
1623  Scalar v = x_view[i];
1624 
1625  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1626  if (ST::magnitude(v.coeff(j)) < tol)
1627  v.fastAccessCoeff(j) = BaseScalar(0.0);
1628  }
1629 
1630  Scalar val = Scalar(0.0);
1631 
1632  for (LocalOrdinal j=0; j<VectorSize; ++j)
1633  if (j+2+row-VectorSize > 0)
1634  val.fastAccessCoeff(j) = BaseScalar(1.0);
1635 
1636  for (LocalOrdinal j=0; j<VectorSize; ++j)
1637  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1638  }
1639 }
1640 
1641 //
1642 // Test Belos GMRES solve for a simple lower-triangular matrix with lucky breakdown with ICGS orthogonalization
1643 //
1645  Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1646 {
1647  using Teuchos::RCP;
1648  using Teuchos::rcp;
1649  using Teuchos::tuple;
1650  using Teuchos::ArrayView;
1651  using Teuchos::Array;
1652  using Teuchos::ArrayRCP;
1653  using Teuchos::ParameterList;
1654 
1655  typedef typename Storage::value_type BaseScalar;
1657 
1658  typedef Teuchos::Comm<int> Tpetra_Comm;
1659  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1660  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1661  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1662  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1663 
1664  // Ensure device is initialized
1665  if ( !Kokkos::is_initialized() )
1666  Kokkos::initialize();
1667 
1668  // Build diagonal matrix
1669  GlobalOrdinal nrow = 20;
1670  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1671  RCP<const Tpetra_Map> map =
1672  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1673  nrow, comm);
1674  RCP<Tpetra_CrsGraph> graph =
1675  rcp(new Tpetra_CrsGraph(map, size_t(1), Tpetra::StaticProfile));
1676  Array<GlobalOrdinal> columnIndices(1);
1677  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1678  const size_t num_my_row = myGIDs.size();
1679  for (size_t i=0; i<num_my_row; ++i) {
1680  const GlobalOrdinal row = myGIDs[i];
1681  columnIndices[0] = row;
1682  size_t ncol = 1;
1683  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1684  }
1685  graph->fillComplete();
1686  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1687 
1688  Array<Scalar> vals(1);
1689  for (size_t i=0; i<num_my_row; ++i) {
1690  const GlobalOrdinal row = myGIDs[i];
1691  columnIndices[0] = row;
1692  vals[0] = Scalar(row+1);
1693  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1694  }
1695  matrix->fillComplete();
1696 
1697  // Fill RHS vector:
1698  // [ 0, 0, ..., 0, 0, 1]
1699  // [ 0, 0, ..., 0, 2, 2]
1700  // [ 0, 0, ..., 3, 3, 3]
1701  // [ 0, 0, ..., 4, 4, 4]
1702  // ...
1703 
1704  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1705  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1706  for (size_t i=0; i<num_my_row; ++i) {
1707  const GlobalOrdinal row = myGIDs[i];
1708  b_view[i] = Scalar(0.0);
1709  for (LocalOrdinal j=0; j<VectorSize; ++j)
1710  if (int(j+2+row-VectorSize) > 0)
1711  b_view[i].fastAccessCoeff(j) = BaseScalar(row+1);
1712  }
1713 
1714  // Solve
1715  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1716 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1717  typedef BaseScalar BelosScalar;
1718 #else
1719  typedef Scalar BelosScalar;
1720 #endif
1721  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1722  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1723  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1724  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1725  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1726  RCP<ParameterList> belosParams = rcp(new ParameterList);
1727  typename ST::magnitudeType tol = 1e-9;
1728  belosParams->set("Flexible Gmres", false);
1729  belosParams->set("Num Blocks", 100);
1730  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1731  belosParams->set("Maximum Iterations", 100);
1732  belosParams->set("Verbosity", 33);
1733  belosParams->set("Output Style", 1);
1734  belosParams->set("Output Frequency", 1);
1735  belosParams->set("Output Stream", out.getOStream());
1736  belosParams->set("Orthogonalization","ICGS");
1737  RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1738  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1739  problem->setProblem();
1740  Belos::ReturnType ret = solver->solve();
1741  TEST_EQUALITY_CONST( ret, Belos::Converged );
1742 
1743 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1744  int numItersExpected = nrow;
1745  int numIters = solver->getNumIters();
1746  out << "numIters = " << numIters << std::endl;
1747  TEST_EQUALITY( numIters, numItersExpected);
1748 
1749  // Get and print number of ensemble iterations
1750  std::vector<int> ensemble_iterations =
1751  static_cast<const Belos::StatusTestImpResNorm<BelosScalar, MV, OP> *>(solver->getResidualStatusTest())->getEnsembleIterations();
1752  out << "Ensemble iterations = ";
1753  for (auto ensemble_iteration : ensemble_iterations)
1754  out << ensemble_iteration << " ";
1755  out << std::endl;
1756 
1757  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1758  if (int(j+1+nrow-VectorSize) > 0) {
1759  TEST_EQUALITY(int(j+1+nrow-VectorSize), ensemble_iterations[j]);
1760  }
1761  else {
1762  TEST_EQUALITY(int(0), ensemble_iterations[j]);
1763  }
1764  }
1765 #endif
1766 
1767  // Check -- Correct answer is:
1768  // [ 0, 0, ..., 0, 0, 1]
1769  // [ 0, 0, ..., 0, 1, 1]
1770  // [ 0, 0, ..., 1, 1, 1]
1771  // [ 0, 0, ..., 1, 1, 1]
1772  // ...
1773  tol = 1000*tol;
1774  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1775  for (size_t i=0; i<num_my_row; ++i) {
1776  const GlobalOrdinal row = myGIDs[i];
1777  Scalar v = x_view[i];
1778 
1779  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1780  if (ST::magnitude(v.coeff(j)) < tol)
1781  v.fastAccessCoeff(j) = BaseScalar(0.0);
1782  }
1783 
1784  Scalar val = Scalar(0.0);
1785 
1786  for (LocalOrdinal j=0; j<VectorSize; ++j)
1787  if (j+2+row-VectorSize > 0)
1788  val.fastAccessCoeff(j) = BaseScalar(1.0);
1789 
1790  for (LocalOrdinal j=0; j<VectorSize; ++j)
1791  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1792  }
1793 }
1794 
1795 //
1796 // Test Belos GMRES solve for a simple lower-triangular matrix with lucky breakdown with IMGS orthogonalization
1797 //
1799  Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1800 {
1801  using Teuchos::RCP;
1802  using Teuchos::rcp;
1803  using Teuchos::tuple;
1804  using Teuchos::ArrayView;
1805  using Teuchos::Array;
1806  using Teuchos::ArrayRCP;
1807  using Teuchos::ParameterList;
1808 
1809  typedef typename Storage::value_type BaseScalar;
1811 
1812  typedef Teuchos::Comm<int> Tpetra_Comm;
1813  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1814  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1815  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1816  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1817 
1818  // Ensure device is initialized
1819  if ( !Kokkos::is_initialized() )
1820  Kokkos::initialize();
1821 
1822  // Build diagonal matrix
1823  GlobalOrdinal nrow = 20;
1824  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1825  RCP<const Tpetra_Map> map =
1826  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1827  nrow, comm);
1828  RCP<Tpetra_CrsGraph> graph =
1829  rcp(new Tpetra_CrsGraph(map, size_t(1), Tpetra::StaticProfile));
1830  Array<GlobalOrdinal> columnIndices(1);
1831  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1832  const size_t num_my_row = myGIDs.size();
1833  for (size_t i=0; i<num_my_row; ++i) {
1834  const GlobalOrdinal row = myGIDs[i];
1835  columnIndices[0] = row;
1836  size_t ncol = 1;
1837  graph->insertGlobalIndices(row, columnIndices(0,ncol));
1838  }
1839  graph->fillComplete();
1840  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
1841 
1842  Array<Scalar> vals(1);
1843  for (size_t i=0; i<num_my_row; ++i) {
1844  const GlobalOrdinal row = myGIDs[i];
1845  columnIndices[0] = row;
1846  vals[0] = Scalar(row+1);
1847  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1848  }
1849  matrix->fillComplete();
1850 
1851  // Fill RHS vector:
1852  // [ 0, 0, ..., 0, 0, 1]
1853  // [ 0, 0, ..., 0, 2, 2]
1854  // [ 0, 0, ..., 3, 3, 3]
1855  // [ 0, 0, ..., 4, 4, 4]
1856  // ...
1857 
1858  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1859  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1860  for (size_t i=0; i<num_my_row; ++i) {
1861  const GlobalOrdinal row = myGIDs[i];
1862  b_view[i] = Scalar(0.0);
1863  for (LocalOrdinal j=0; j<VectorSize; ++j)
1864  if (int(j+2+row-VectorSize) > 0)
1865  b_view[i].fastAccessCoeff(j) = BaseScalar(row+1);
1866  }
1867 
1868  // Solve
1869  typedef Teuchos::ScalarTraits<BaseScalar> ST;
1870 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1871  typedef BaseScalar BelosScalar;
1872 #else
1873  typedef Scalar BelosScalar;
1874 #endif
1875  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1876  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1877  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1878  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1879  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
1880  RCP<ParameterList> belosParams = rcp(new ParameterList);
1881  typename ST::magnitudeType tol = 1e-9;
1882  belosParams->set("Flexible Gmres", false);
1883  belosParams->set("Num Blocks", 100);
1884  belosParams->set("Convergence Tolerance", BelosScalar(tol));
1885  belosParams->set("Maximum Iterations", 100);
1886  belosParams->set("Verbosity", 33);
1887  belosParams->set("Output Style", 1);
1888  belosParams->set("Output Frequency", 1);
1889  belosParams->set("Output Stream", out.getOStream());
1890  belosParams->set("Orthogonalization","IMGS");
1891  RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1892  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1893  problem->setProblem();
1894  Belos::ReturnType ret = solver->solve();
1895  TEST_EQUALITY_CONST( ret, Belos::Converged );
1896 
1897 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1898  int numItersExpected = nrow;
1899  int numIters = solver->getNumIters();
1900  out << "numIters = " << numIters << std::endl;
1901  TEST_EQUALITY( numIters, numItersExpected);
1902 
1903  // Get and print number of ensemble iterations
1904  std::vector<int> ensemble_iterations =
1905  static_cast<const Belos::StatusTestImpResNorm<BelosScalar, MV, OP> *>(solver->getResidualStatusTest())->getEnsembleIterations();
1906  out << "Ensemble iterations = ";
1907  for (auto ensemble_iteration : ensemble_iterations)
1908  out << ensemble_iteration << " ";
1909  out << std::endl;
1910 
1911  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1912  if (int(j+1+nrow-VectorSize) > 0) {
1913  TEST_EQUALITY(int(j+1+nrow-VectorSize), ensemble_iterations[j]);
1914  }
1915  else {
1916  TEST_EQUALITY(int(0), ensemble_iterations[j]);
1917  }
1918  }
1919 #endif
1920 
1921  // Check -- Correct answer is:
1922  // [ 0, 0, ..., 0, 0, 1]
1923  // [ 0, 0, ..., 0, 1, 1]
1924  // [ 0, 0, ..., 1, 1, 1]
1925  // [ 0, 0, ..., 1, 1, 1]
1926  // ...
1927  tol = 1000*tol;
1928  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1929  for (size_t i=0; i<num_my_row; ++i) {
1930  const GlobalOrdinal row = myGIDs[i];
1931  Scalar v = x_view[i];
1932 
1933  for (LocalOrdinal j=0; j<VectorSize; ++j) {
1934  if (ST::magnitude(v.coeff(j)) < tol)
1935  v.fastAccessCoeff(j) = BaseScalar(0.0);
1936  }
1937 
1938  Scalar val = Scalar(0.0);
1939 
1940  for (LocalOrdinal j=0; j<VectorSize; ++j)
1941  if (j+2+row-VectorSize > 0)
1942  val.fastAccessCoeff(j) = BaseScalar(1.0);
1943 
1944  for (LocalOrdinal j=0; j<VectorSize; ++j)
1945  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
1946  }
1947 }
1948 
1949 #else
1950 
1952  Tpetra_CrsMatrix_MP, BelosGMRES, Storage, LocalOrdinal, GlobalOrdinal, Node )
1953 {}
1954 
1956  Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1957 {}
1958 
1960  Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1961 {}
1962 
1964  Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, Storage, LocalOrdinal, GlobalOrdinal, Node )
1965 {}
1966 
1967 #endif
1968 
1969 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
1970 
1971 //
1972 // Test Belos GMRES solve with Ifpack2 RILUK preconditioning for a
1973 // simple banded upper-triangular matrix
1974 //
1976  Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
1977 {
1978  using Teuchos::RCP;
1979  using Teuchos::rcp;
1980  using Teuchos::ArrayView;
1981  using Teuchos::Array;
1982  using Teuchos::ArrayRCP;
1983  using Teuchos::ParameterList;
1984 
1985  typedef typename Storage::value_type BaseScalar;
1987 
1988  typedef Teuchos::Comm<int> Tpetra_Comm;
1989  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1990  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1991  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1992  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1993 
1994  // Ensure device is initialized
1995  if ( !Kokkos::is_initialized() )
1996  Kokkos::initialize();
1997 
1998  // Build banded matrix
1999  GlobalOrdinal nrow = 10;
2000  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2001  RCP<const Tpetra_Map> map =
2002  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2003  nrow, comm);
2004  RCP<Tpetra_CrsGraph> graph =
2005  rcp(new Tpetra_CrsGraph(map, size_t(2), Tpetra::StaticProfile));
2006  Array<GlobalOrdinal> columnIndices(2);
2007  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2008  const size_t num_my_row = myGIDs.size();
2009  for (size_t i=0; i<num_my_row; ++i) {
2010  const GlobalOrdinal row = myGIDs[i];
2011  columnIndices[0] = row;
2012  size_t ncol = 1;
2013  if (row != nrow-1) {
2014  columnIndices[1] = row+1;
2015  ncol = 2;
2016  }
2017  graph->insertGlobalIndices(row, columnIndices(0,ncol));
2018  }
2019  graph->fillComplete();
2020  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
2021 
2022  // Set values in matrix
2023  Array<Scalar> vals(2);
2024  Scalar val(VectorSize, BaseScalar(0.0));
2025  for (size_t i=0; i<num_my_row; ++i) {
2026  const GlobalOrdinal row = myGIDs[i];
2027  columnIndices[0] = row;
2028  for (LocalOrdinal j=0; j<VectorSize; ++j)
2029  val.fastAccessCoeff(j) = j+1;
2030  vals[0] = val;
2031  size_t ncol = 1;
2032 
2033  if (row != nrow-1) {
2034  columnIndices[1] = row+1;
2035  for (LocalOrdinal j=0; j<VectorSize; ++j)
2036  val.fastAccessCoeff(j) = j+1;
2037  vals[1] = val;
2038  ncol = 2;
2039  }
2040  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
2041  }
2042  matrix->fillComplete();
2043 
2044  // Fill RHS vector
2045  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2046  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2047  for (size_t i=0; i<num_my_row; ++i) {
2048  b_view[i] = Scalar(1.0);
2049  }
2050 
2051  // Create preconditioner
2052  typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
2053  Ifpack2::Factory factory;
2054  RCP<Prec> M = factory.create<Tpetra_CrsMatrix>("RILUK", matrix);
2055  M->initialize();
2056  M->compute();
2057 
2058  // Solve
2059  typedef Teuchos::ScalarTraits<BaseScalar> ST;
2060 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2061  typedef BaseScalar BelosScalar;
2062 #else
2063  typedef Scalar BelosScalar;
2064 #endif
2065  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2066  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2067  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2068  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2069  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
2070  problem->setRightPrec(M);
2071  problem->setProblem();
2072  RCP<ParameterList> belosParams = rcp(new ParameterList);
2073  typename ST::magnitudeType tol = 1e-9;
2074  belosParams->set("Flexible Gmres", false);
2075  belosParams->set("Num Blocks", 100);
2076  belosParams->set("Convergence Tolerance", BelosScalar(tol));
2077  belosParams->set("Maximum Iterations", 100);
2078  belosParams->set("Verbosity", 33);
2079  belosParams->set("Output Style", 1);
2080  belosParams->set("Output Frequency", 1);
2081  belosParams->set("Output Stream", out.getOStream());
2082  //belosParams->set("Orthogonalization", "TSQR");
2083  RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2084  rcp(new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2085  Belos::ReturnType ret = solver->solve();
2086  TEST_EQUALITY_CONST( ret, Belos::Converged );
2087 
2088  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2089  // Teuchos::VERB_EXTREME);
2090 
2091  // Check -- Correct answer is:
2092  // [ 0, 0, ..., 0 ]
2093  // [ 1, 1/2, ..., 1/VectorSize ]
2094  // [ 0, 0, ..., 0 ]
2095  // [ 1, 1/2, ..., 1/VectorSize ]
2096  // ....
2097  tol = 1000*tol;
2098  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2099  for (size_t i=0; i<num_my_row; ++i) {
2100  const GlobalOrdinal row = myGIDs[i];
2101  if (row % 2) {
2102  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2103  val.fastAccessCoeff(j) = BaseScalar(1.0) / BaseScalar(j+1);
2104  }
2105  }
2106  else
2107  val = Scalar(VectorSize, BaseScalar(0.0));
2108  TEST_EQUALITY( x_view[i].size(), VectorSize );
2109 
2110  // Set small values to zero
2111  Scalar v = x_view[i];
2112  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2113  if (ST::magnitude(v.coeff(j)) < tol)
2114  v.fastAccessCoeff(j) = BaseScalar(0.0);
2115  }
2116 
2117  for (LocalOrdinal j=0; j<VectorSize; ++j)
2118  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
2119  }
2120 }
2121 
2122 #else
2123 
2125  Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, Storage, LocalOrdinal, GlobalOrdinal, Node )
2126 {}
2127 
2128 #endif
2129 
2130 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2)
2131 
2132 //
2133 // Test Belos CG solve with MueLu preconditioning for a 1-D Laplacian matrix
2134 //
2136  Tpetra_CrsMatrix_MP, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
2137 {
2138  using Teuchos::RCP;
2139  using Teuchos::rcp;
2140  using Teuchos::ArrayView;
2141  using Teuchos::Array;
2142  using Teuchos::ArrayRCP;
2143  using Teuchos::ParameterList;
2144  using Teuchos::getParametersFromXmlFile;
2145 
2146  typedef typename Storage::value_type BaseScalar;
2148 
2149  typedef Teuchos::Comm<int> Tpetra_Comm;
2150  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2151  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2152  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2153  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2154 
2155  // Ensure device is initialized
2156  if ( !Kokkos::is_initialized() )
2157  Kokkos::initialize();
2158 
2159  // 1-D Laplacian matrix
2160  GlobalOrdinal nrow = 50;
2161  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
2162  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2163  RCP<const Tpetra_Map> map =
2164  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2165  nrow, comm);
2166  RCP<Tpetra_CrsGraph> graph =
2167  rcp(new Tpetra_CrsGraph(map, size_t(3), Tpetra::StaticProfile));
2168  Array<GlobalOrdinal> columnIndices(3);
2169  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2170  const size_t num_my_row = myGIDs.size();
2171  for (size_t i=0; i<num_my_row; ++i) {
2172  const GlobalOrdinal row = myGIDs[i];
2173  if (row == 0 || row == nrow-1) { // Boundary nodes
2174  columnIndices[0] = row;
2175  graph->insertGlobalIndices(row, columnIndices(0,1));
2176  }
2177  else { // Interior nodes
2178  columnIndices[0] = row-1;
2179  columnIndices[1] = row;
2180  columnIndices[2] = row+1;
2181  graph->insertGlobalIndices(row, columnIndices(0,3));
2182  }
2183  }
2184  graph->fillComplete();
2185  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
2186 
2187  // Set values in matrix
2188  Array<Scalar> vals(3);
2189  Scalar a_val(VectorSize, BaseScalar(0.0));
2190  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2191  a_val.fastAccessCoeff(j) =
2192  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize);
2193  }
2194  for (size_t i=0; i<num_my_row; ++i) {
2195  const GlobalOrdinal row = myGIDs[i];
2196  if (row == 0 || row == nrow-1) { // Boundary nodes
2197  columnIndices[0] = row;
2198  vals[0] = Scalar(1.0);
2199  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2200  }
2201  else {
2202  columnIndices[0] = row-1;
2203  columnIndices[1] = row;
2204  columnIndices[2] = row+1;
2205  vals[0] = Scalar(-1.0) * a_val;
2206  vals[1] = Scalar(2.0) * a_val;
2207  vals[2] = Scalar(-1.0) * a_val;
2208  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2209  }
2210  }
2211  matrix->fillComplete();
2212 
2213  // Fill RHS vector
2214  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2215  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2216  Scalar b_val(VectorSize, BaseScalar(0.0));
2217  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2218  b_val.fastAccessCoeff(j) =
2219  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize);
2220  }
2221  for (size_t i=0; i<num_my_row; ++i) {
2222  const GlobalOrdinal row = myGIDs[i];
2223  if (row == 0 || row == nrow-1)
2224  b_view[i] = Scalar(0.0);
2225  else
2226  b_view[i] = -Scalar(b_val * h * h);
2227  }
2228 
2229  // Create preconditioner
2230  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2231  RCP<ParameterList> muelu_params =
2232  getParametersFromXmlFile("muelu_cheby.xml");
2233  RCP<OP> matrix_op = matrix;
2234  RCP<OP> M =
2235  MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
2236 
2237  // Solve
2238  typedef Teuchos::ScalarTraits<BaseScalar> ST;
2239 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2240  typedef BaseScalar BelosScalar;
2241 #else
2242  typedef Scalar BelosScalar;
2243 #endif
2244  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2245  typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2246  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2247  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
2248  problem->setRightPrec(M);
2249  problem->setProblem();
2250  RCP<ParameterList> belosParams = rcp(new ParameterList);
2251  typename ST::magnitudeType tol = 1e-9;
2252  belosParams->set("Num Blocks", 100);
2253  belosParams->set("Convergence Tolerance", BelosScalar(tol));
2254  belosParams->set("Maximum Iterations", 100);
2255  belosParams->set("Verbosity", 33);
2256  belosParams->set("Output Style", 1);
2257  belosParams->set("Output Frequency", 1);
2258  belosParams->set("Output Stream", out.getOStream());
2259  // Turn off residual scaling so we can see some variation in the number
2260  // of iterations across the ensemble when not doing ensemble reductions
2261  belosParams->set("Implicit Residual Scaling", "None");
2262 
2263  RCP<Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true> > solver =
2264  rcp(new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true>(problem, belosParams));
2265  Belos::ReturnType ret = solver->solve();
2266  TEST_EQUALITY_CONST( ret, Belos::Converged );
2267 
2268 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
2269  // Get and print number of ensemble iterations
2270  std::vector<int> ensemble_iterations =
2271  solver->getResidualStatusTest()->getEnsembleIterations();
2272  out << "Ensemble iterations = ";
2273  for (int i=0; i<VectorSize; ++i)
2274  out << ensemble_iterations[i] << " ";
2275  out << std::endl;
2276 #endif
2277 
2278  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2279  // Teuchos::VERB_EXTREME);
2280 
2281  // Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1)
2282  tol = 1000*tol;
2283  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2284  Scalar val(VectorSize, BaseScalar(0.0));
2285  for (size_t i=0; i<num_my_row; ++i) {
2286  const GlobalOrdinal row = myGIDs[i];
2287  BaseScalar xx = row * h;
2288  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2289  val.fastAccessCoeff(j) =
2290  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0));
2291  }
2292  TEST_EQUALITY( x_view[i].size(), VectorSize );
2293 
2294  // Set small values to zero
2295  Scalar v = x_view[i];
2296  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2297  if (ST::magnitude(v.coeff(j)) < tol)
2298  v.fastAccessCoeff(j) = BaseScalar(0.0);
2299  if (ST::magnitude(val.coeff(j)) < tol)
2300  val.fastAccessCoeff(j) = BaseScalar(0.0);
2301  }
2302 
2303  for (LocalOrdinal j=0; j<VectorSize; ++j)
2304  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
2305  }
2306 
2307 }
2308 
2309 #else
2310 
2312  Tpetra_CrsMatrix_MP, BelosCG_Muelu, Storage, LocalOrdinal, GlobalOrdinal, Node )
2313 {}
2314 
2315 #endif
2316 
2317 #if defined(HAVE_STOKHOS_AMESOS2)
2318 
2319 //
2320 // Test Amesos2 solve for a 1-D Laplacian matrix
2321 //
2323  Tpetra_CrsMatrix_MP, Amesos2, Storage, LocalOrdinal, GlobalOrdinal, Node )
2324 {
2325  using Teuchos::RCP;
2326  using Teuchos::rcp;
2327  using Teuchos::ArrayView;
2328  using Teuchos::Array;
2329  using Teuchos::ArrayRCP;
2330  using Teuchos::ParameterList;
2331 
2332  typedef typename Storage::value_type BaseScalar;
2334 
2335  typedef Teuchos::Comm<int> Tpetra_Comm;
2336  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2337  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2338  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
2339  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2340  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2341 
2342  // Ensure device is initialized
2343  if ( !Kokkos::is_initialized() )
2344  Kokkos::initialize();
2345 
2346  // 1-D Laplacian matrix
2347  GlobalOrdinal nrow = 50;
2348  BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1);
2349  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2350  RCP<const Tpetra_Map> map =
2351  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2352  nrow, comm);
2353  RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map, size_t(3));
2354  Array<GlobalOrdinal> columnIndices(3);
2355  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2356  const size_t num_my_row = myGIDs.size();
2357  for (size_t i=0; i<num_my_row; ++i) {
2358  const GlobalOrdinal row = myGIDs[i];
2359  if (row == 0 || row == nrow-1) { // Boundary nodes
2360  columnIndices[0] = row;
2361  graph->insertGlobalIndices(row, columnIndices(0,1));
2362  }
2363  else { // Interior nodes
2364  columnIndices[0] = row-1;
2365  columnIndices[1] = row;
2366  columnIndices[2] = row+1;
2367  graph->insertGlobalIndices(row, columnIndices(0,3));
2368  }
2369  }
2370  graph->fillComplete();
2371  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
2372 
2373  // Set values in matrix
2374  Array<Scalar> vals(3);
2375  Scalar a_val(VectorSize, BaseScalar(0.0));
2376  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2377  a_val.fastAccessCoeff(j) =
2378  BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize);
2379  }
2380  for (size_t i=0; i<num_my_row; ++i) {
2381  const GlobalOrdinal row = myGIDs[i];
2382  if (row == 0 || row == nrow-1) { // Boundary nodes
2383  columnIndices[0] = row;
2384  vals[0] = Scalar(1.0);
2385  matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2386  }
2387  else {
2388  columnIndices[0] = row-1;
2389  columnIndices[1] = row;
2390  columnIndices[2] = row+1;
2391  vals[0] = Scalar(-1.0) * a_val;
2392  vals[1] = Scalar(2.0) * a_val;
2393  vals[2] = Scalar(-1.0) * a_val;
2394  matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2395  }
2396  }
2397  matrix->fillComplete();
2398 
2399  // Fill RHS vector
2400  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2401  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2402  Scalar b_val(VectorSize, BaseScalar(0.0));
2403  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2404  b_val.fastAccessCoeff(j) =
2405  BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize);
2406  }
2407  for (size_t i=0; i<num_my_row; ++i) {
2408  const GlobalOrdinal row = myGIDs[i];
2409  if (row == 0 || row == nrow-1)
2410  b_view[i] = Scalar(0.0);
2411  else
2412  b_view[i] = -Scalar(b_val * h * h);
2413  }
2414 
2415  // Solve
2416  typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2417  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2418  std::string solver_name;
2419 #if defined(HAVE_AMESOS2_BASKER)
2420  solver_name = "basker";
2421 #elif defined(HAVE_AMESOS2_KLU2)
2422  solver_name = "klu2";
2423 #elif defined(HAVE_AMESOS2_SUPERLUDIST)
2424  solver_name = "superlu_dist";
2425 #elif defined(HAVE_AMESOS2_SUPERLUMT)
2426  solver_name = "superlu_mt";
2427 #elif defined(HAVE_AMESOS2_SUPERLU)
2428  solver_name = "superlu";
2429 #elif defined(HAVE_AMESOS2_PARDISO_MKL)
2430  solver_name = "pardisomkl";
2431 #elif defined(HAVE_AMESOS2_LAPACK)
2432  solver_name = "lapack";
2433 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2434  solver_name = "lapack";
2435 #else
2436  // if there are no solvers, we just return as a successful test
2437  success = true;
2438  return;
2439 #endif
2440  out << "Solving linear system with " << solver_name << std::endl;
2441  RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2442  solver_name, matrix, x, b);
2443  solver->solve();
2444 
2445  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
2446  // Teuchos::VERB_EXTREME);
2447 
2448  // Check -- For a*y'' = b, correct answer is y = 0.5 *(b/a) * x * (x-1)
2449  typedef Teuchos::ScalarTraits<BaseScalar> ST;
2450  typename ST::magnitudeType tol = 1e-9;
2451  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2452  Scalar val(VectorSize, BaseScalar(0.0));
2453  for (size_t i=0; i<num_my_row; ++i) {
2454  const GlobalOrdinal row = myGIDs[i];
2455  BaseScalar xx = row * h;
2456  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2457  val.fastAccessCoeff(j) =
2458  BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0));
2459  }
2460  TEST_EQUALITY( x_view[i].size(), VectorSize );
2461 
2462  // Set small values to zero
2463  Scalar v = x_view[i];
2464  for (LocalOrdinal j=0; j<VectorSize; ++j) {
2465  if (ST::magnitude(v.coeff(j)) < tol)
2466  v.fastAccessCoeff(j) = BaseScalar(0.0);
2467  if (ST::magnitude(val.coeff(j)) < tol)
2468  val.fastAccessCoeff(j) = BaseScalar(0.0);
2469  }
2470 
2471  for (LocalOrdinal j=0; j<VectorSize; ++j)
2472  TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol);
2473  }
2474 }
2475 
2476 #else
2477 
2479  Tpetra_CrsMatrix_MP, Amesos2, Storage, LocalOrdinal, GlobalOrdinal, Node )
2480 {}
2481 
2482 #endif
2483 
2484 #define CRSMATRIX_MP_VECTOR_TESTS_SLGN(S, LO, GO, N) \
2485  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorAdd, S, LO, GO, N ) \
2486  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorDot, S, LO, GO, N ) \
2487  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorAdd, S, LO, GO, N ) \
2488  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDot, S, LO, GO, N ) \
2489  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDotSub, S, LO, GO, N ) \
2490  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixVectorMultiply, S, LO, GO, N ) \
2491  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2492  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, WrappedDualView, S, LO, GO, N ) \
2493  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Flatten, S, LO, GO, N ) \
2494  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimpleCG, S, LO, GO, N ) \
2495  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimplePCG_Muelu, S, LO, GO, N ) \
2496  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES, S, LO, GO, N ) \
2497  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, S, LO, GO, N ) \
2498  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, S, LO, GO, N ) \
2499  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, S, LO, GO, N ) \
2500  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, S, LO, GO, N ) \
2501  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosCG_Muelu, S, LO, GO, N ) \
2502  TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2, S, LO, GO, N )
2503 
2504 #define CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) \
2505  typedef Stokhos::DeviceForNode<N>::type Device; \
2506  typedef Stokhos::StaticFixedStorage<int,double,VectorSize,Device::execution_space> SFS; \
2507  using default_global_ordinal_type = ::Tpetra::Map<>::global_ordinal_type; \
2508  using default_local_ordinal_type = ::Tpetra::Map<>::local_ordinal_type; \
2509  CRSMATRIX_MP_VECTOR_TESTS_SLGN(SFS, default_local_ordinal_type, default_global_ordinal_type, N)
2510 
2511 #define CRSMATRIX_MP_VECTOR_TESTS_N(N) \
2512  CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N)
2513 
2514 // Disabling testing of dynamic storage -- we don't really need it
2515  // typedef Stokhos::DynamicStorage<int,double,Device> DS;
2516  // using default_global_ordinal_type = ::Tpetra::Map<>::global_ordinal_type;
2517  // using default_local_ordinal_type = ::Tpetra::Map<>::local_ordinal_type;
2518  // CRSMATRIX_MP_VECTOR_TESTS_SLGN(DS, default_global_ordinal_type, default_local_ordinal_type, N)
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
scalar generate_multi_vector_coefficient(const ordinal nFEM, const ordinal nVec, const ordinal nStoch, const ordinal iColFEM, const ordinal iVec, const ordinal iStoch)
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j) - expr2.val(j)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec_const, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
bool PCG_Solve(const Matrix &A, Vector &x, const Vector &b, const Prec &M, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
KokkosClassic::DefaultNode::DefaultNodeType Node
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_flat_mp_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_range_map, const LocalOrdinal block_size)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const LocalOrdinal block_size)
expr val()
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)
BelosGMRES
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)