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