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