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