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