Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_KokkosArrayKernelsUnitTest.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#ifndef STOKHOS_KOKKOS_CORE_KERNELS_UNIT_TEST_HPP
43#define STOKHOS_KOKKOS_CORE_KERNELS_UNIT_TEST_HPP
44
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestHelpers.hpp"
47#include "Teuchos_TestForException.hpp"
48
49#include "Stokhos_Epetra.hpp"
51#include "EpetraExt_BlockUtility.h"
53
54#ifdef HAVE_MPI
55#include "Epetra_MpiComm.h"
56#else
57#include "Epetra_SerialComm.h"
58#endif
59
60#include "Kokkos_Macros.hpp"
61
62#include "Stokhos_Update.hpp"
63#include "Stokhos_CrsMatrix.hpp"
75
76#ifdef HAVE_STOKHOS_KOKKOSLINALG
77#include "KokkosSparse_CrsMatrix.hpp"
78#include "KokkosSparse_spmv.hpp"
79#include "KokkosBlas1_update.hpp"
80#endif
81
83
84using Teuchos::rcp;
85using Teuchos::RCP;
86using Teuchos::Array;
87using Teuchos::ParameterList;
88
89template< typename IntType >
90inline
91IntType map_fem_graph_coord( const IntType & N ,
92 const IntType & i ,
93 const IntType & j ,
94 const IntType & k )
95{
96 return k + N * ( j + N * i );
97}
98
99template < typename ordinal >
100inline
101ordinal generate_fem_graph( ordinal N ,
102 std::vector< std::vector<ordinal> > & graph )
103{
104 graph.resize( N * N * N , std::vector<ordinal>() );
105
106 ordinal total = 0 ;
107
108 for ( int i = 0 ; i < (int) N ; ++i ) {
109 for ( int j = 0 ; j < (int) N ; ++j ) {
110 for ( int k = 0 ; k < (int) N ; ++k ) {
111
112 const ordinal row = map_fem_graph_coord((int)N,i,j,k);
113
114 graph[row].reserve(27);
115
116 for ( int ii = -1 ; ii < 2 ; ++ii ) {
117 for ( int jj = -1 ; jj < 2 ; ++jj ) {
118 for ( int kk = -1 ; kk < 2 ; ++kk ) {
119 if ( 0 <= i + ii && i + ii < (int) N &&
120 0 <= j + jj && j + jj < (int) N &&
121 0 <= k + kk && k + kk < (int) N ) {
122 ordinal col = map_fem_graph_coord((int)N,i+ii,j+jj,k+kk);
123
124 graph[row].push_back(col);
125 }
126 }}}
127 total += graph[row].size();
128 }}}
129
130 return total ;
131}
132
133template <typename scalar, typename ordinal>
134inline
135scalar generate_matrix_coefficient( const ordinal nFEM ,
136 const ordinal nStoch ,
137 const ordinal iRowFEM ,
138 const ordinal iColFEM ,
139 const ordinal iStoch )
140{
141 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
142 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
143
144 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
145
146 return A_fem + A_stoch ;
147 //return 1.0;
148}
149
150template <typename scalar, typename ordinal>
151inline
152scalar generate_vector_coefficient( const ordinal nFEM ,
153 const ordinal nStoch ,
154 const ordinal iColFEM ,
155 const ordinal iStoch )
156{
157 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
158 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
159 return X_fem + X_stoch ;
160 //return 1.0;
161}
162
163template <typename Device>
165 typedef double value_type ;
168 //typedef Stokhos::CompletePolynomialBasis<int,value_type> product_basis_type;
172
175 std::vector< std::vector<int> > fem_graph ;
176 RCP< product_basis_type> basis;
177 RCP<Cijk_type> Cijk;
178 RCP<const Epetra_Comm> globalComm;
179 RCP<Stokhos::EpetraVectorOrthogPoly> sg_x, sg_y;
180 RCP<Stokhos::ProductEpetraVector> sg_y_commuted;
181 Teuchos::Array<int> perm, inv_perm;
182
183 // Can't be a constructor because MPI will not be initialized
184 void setup(int p_ = 5, int d_ = 2) {
185
186 p = p_;
187 d = d_;
188 nGrid = 3;
189 rel_tol = 1e-12;
190 abs_tol = 1e-12;
191
192 // Create a communicator for Epetra objects
193#ifdef HAVE_MPI
194 globalComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
195#else
196 globalComm = rcp(new Epetra_SerialComm);
197#endif
198 //int MyPID = globalComm->MyPID();
199
200 //------------------------------
201 // Generate FEM graph:
202
205
206 // Create Stochastic Galerkin basis and expansion
207 Array< RCP<const abstract_basis_type> > bases(d);
208 for (int i=0; i<d; i++)
209 bases[i] = rcp(new basis_type(p,true));
210 basis = rcp(new product_basis_type(bases, 1e-12));
211 stoch_length = basis->size();
212 Cijk = basis->computeTripleProductTensor();
213
214 // Create stochastic parallel distribution
215 ParameterList parallelParams;
216 RCP<Stokhos::ParallelData> sg_parallel_data =
217 rcp(new Stokhos::ParallelData(basis, Cijk, globalComm, parallelParams));
218 RCP<const EpetraExt::MultiComm> sg_comm =
219 sg_parallel_data->getMultiComm();
220 RCP<const Epetra_Comm> app_comm =
221 sg_parallel_data->getSpatialComm();
222 RCP<const Stokhos::EpetraSparse3Tensor> epetraCijk =
223 sg_parallel_data->getEpetraCijk();
224 RCP<const Epetra_BlockMap> stoch_row_map =
225 epetraCijk->getStochasticRowMap();
226
227 //------------------------------
228
229 // Generate Epetra objects
230 RCP<const Epetra_Map> x_map =
231 rcp(new Epetra_Map(fem_length, 0, *app_comm));
232 RCP<const Epetra_Map> sg_x_map =
233 rcp(EpetraExt::BlockUtility::GenerateBlockMap(
234 *x_map, *stoch_row_map, *sg_comm));
235
236 RCP<Epetra_CrsGraph> graph = rcp(new Epetra_CrsGraph(Copy, *x_map, 27));
237 int *my_GIDs = x_map->MyGlobalElements();
238 int num_my_GIDs = x_map->NumMyElements();
239 for (int i=0; i<num_my_GIDs; ++i) {
240 int row = my_GIDs[i];
241 int num_indices = fem_graph[row].size();
242 int *indices = &(fem_graph[row][0]);
243 graph->InsertGlobalIndices(row, num_indices, indices);
244 }
245 graph->FillComplete();
246
247 RCP<ParameterList> sg_op_params = rcp(new ParameterList);
248 RCP<Stokhos::MatrixFreeOperator> sg_A =
249 rcp(new Stokhos::MatrixFreeOperator(sg_comm, basis, epetraCijk,
250 x_map, x_map, sg_x_map, sg_x_map,
251 sg_op_params));
252 RCP<Epetra_BlockMap> sg_A_overlap_map =
253 rcp(new Epetra_LocalMap(
254 stoch_length, 0, *(sg_parallel_data->getStochasticComm())));
255 RCP< Stokhos::EpetraOperatorOrthogPoly > A_sg_blocks =
257 basis, sg_A_overlap_map, x_map, x_map, sg_x_map, sg_comm));
258 for (int i=0; i<stoch_length; i++) {
259 RCP<Epetra_CrsMatrix> A = rcp(new Epetra_CrsMatrix(Copy, *graph));
260
261 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < fem_length ; ++iRowFEM ) {
262 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
263 const int iColFEM = fem_graph[iRowFEM][iRowEntryFEM] ;
264
265 double v = generate_matrix_coefficient<double>(
266 fem_length , stoch_length , iRowFEM , iColFEM , i );
267 A->ReplaceGlobalValues(iRowFEM, 1, &v, &iColFEM);
268 }
269 }
270 A->FillComplete();
271 A_sg_blocks->setCoeffPtr(i, A);
272 }
273 sg_A->setupOperator(A_sg_blocks);
274
276 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
278 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
279
280 // Initialize vectors
281 for (int iColFEM=0; iColFEM < fem_length; ++iColFEM ) {
282 for (int iColStoch=0 ; iColStoch < stoch_length; ++iColStoch ) {
283 (*sg_x)[iColStoch][iColFEM] = generate_vector_coefficient<double>(
284 fem_length , stoch_length , iColFEM , iColStoch );
285 }
286 }
287 sg_y->init(0.0);
288
289 // Apply operator
290 sg_A->Apply( *(sg_x->getBlockVector()), *(sg_y->getBlockVector()) );
291
292 // Transpose y to commuted layout
294 x_map, stoch_row_map, sg_comm));
295 for (int block=0; block<stoch_length; ++block) {
296 for (int i=0; i<fem_length; ++i)
297 (*sg_y_commuted)[i][block] = (*sg_y)[block][i];
298 }
299
300 typedef typename Kokkos::ViewTraits<value_type,Device,void,void>::host_mirror_space host_device;
302 typedef Stokhos::StochasticProductTensor< value_type , tensor_type ,
303 host_device > stoch_tensor_type ;
304
305 stoch_tensor_type tensor =
306 Stokhos::create_stochastic_product_tensor< tensor_type >( *basis, *Cijk );
307 stoch_length_aligned = tensor.aligned_dimension();
308
309 perm.resize(stoch_length);
310 inv_perm.resize(stoch_length);
311 for (int i=0; i<stoch_length; ++i) {
313 for (int j=0; j<d; ++j)
314 term[j] = tensor.bases_degree(i,j);
315 int idx = basis->index(term);
316 perm[idx] = i;
317 inv_perm[i] = idx;
318 }
319 }
320
321 template <typename vec_type>
322 bool test_original(const std::vector<vec_type>& y,
323 Teuchos::FancyOStream& out) const {
324 bool success = true;
325 for (int block=0; block<stoch_length; ++block) {
326 for (int i=0; i<fem_length; ++i) {
327 double diff = std::abs( (*sg_y)[block][i] - y[block][i] );
328 double tol = rel_tol*std::abs((*sg_y)[block][i]) + abs_tol;
329 bool s = diff < tol;
330 if (!s) {
331 out << "y_expected[" << block << "][" << i << "] - "
332 << "y[" << block << "][" << i << "] = " << (*sg_y)[block][i]
333 << " - " << y[block][i] << " == "
334 << diff << " < " << tol << " : failed"
335 << std::endl;
336 }
337 success = success && s;
338 }
339 }
340
341 return success;
342 }
343
344 template <typename multi_vec_type>
345 bool test_original(const multi_vec_type& y,
346 Teuchos::FancyOStream& out) const {
347 bool success = true;
348 for (int block=0; block<stoch_length; ++block) {
349 for (int i=0; i<fem_length; ++i) {
350 double diff = std::abs( (*sg_y)[block][i] - y(i,block) );
351 double tol = rel_tol*std::abs((*sg_y)[block][i]) + abs_tol;
352 bool s = diff < tol;
353 if (!s) {
354 out << "y_expected[" << block << "][" << i << "] - "
355 << "y(" << i << "," << block << ") = " << (*sg_y)[block][i]
356 << " - " << y(i,block) << " == "
357 << diff << " < " << tol << " : failed"
358 << std::endl;
359 }
360 success = success && s;
361 }
362 }
363
364 return success;
365 }
366
367 template <typename vec_type>
368 bool test_commuted_perm(const vec_type& y,
369 Teuchos::FancyOStream& out) const {
370 bool success = true;
371 for (int block=0; block<stoch_length; ++block) {
372 int b = perm[block];
373 for (int i=0; i<fem_length; ++i) {
374 double diff = std::abs( (*sg_y)[block][i] - y(b,i) );
375 double tol = rel_tol*std::abs((*sg_y)[block][i]) + abs_tol;
376 bool s = diff < tol;
377 if (!s) {
378 out << "y_expected[" << block << "][" << i << "] - "
379 << "y(" << b << "," << i << ") = " << (*sg_y)[block][i]
380 << " - " << y(b,i) << " == "
381 << diff << " < " << tol << " : failed"
382 << std::endl;
383 }
384 success = success && s;
385 }
386 }
387
388 return success;
389 }
390
391 template <typename vec_type>
392 bool test_commuted(const vec_type& y,
393 Teuchos::FancyOStream& out) const {
394 bool success = true;
395 for (int block=0; block<stoch_length; ++block) {
396 int b = block;
397 for (int i=0; i<fem_length; ++i) {
398 double diff = std::abs( (*sg_y)[block][i] - y(b,i) );
399 double tol = rel_tol*std::abs((*sg_y)[block][i]) + abs_tol;
400 bool s = diff < tol;
401 if (!s) {
402 out << "y_expected[" << block << "][" << i << "] - "
403 << "y(" << b << "," << i << ") = " << (*sg_y)[block][i] << " - "
404 << y(b,i) << " == "
405 << diff << " < " << tol << " : failed"
406 << std::endl;
407 }
408 success = success && s;
409 }
410 }
411
412 return success;
413 }
414
415 template <typename vec_type>
416 bool test_commuted_flat(const vec_type& y,
417 Teuchos::FancyOStream& out) const {
418 bool success = true;
419 for (int block=0; block<stoch_length; ++block) {
420 int b = block;
421 for (int i=0; i<fem_length; ++i) {
422 double diff = std::abs( (*sg_y)[block][i] - y(i*stoch_length+b) );
423 double tol = rel_tol*std::abs((*sg_y)[block][i]) + abs_tol;
424 bool s = diff < tol;
425 if (!s) {
426 out << "y_expected[" << block << "][" << i << "] - "
427 << "y(" << b << "," << i << ") == "
428 << diff << " < " << tol << " : failed"
429 << std::endl;
430 }
431 success = success && s;
432 }
433 }
434
435 return success;
436 }
437
438 template <typename vec_type>
439 bool test_original_flat(const vec_type& y,
440 Teuchos::FancyOStream& out) const {
441 bool success = true;
442 for (int block=0; block<stoch_length; ++block) {
443 int b = block;
444 for (int i=0; i<fem_length; ++i) {
445 double diff = std::abs( (*sg_y)[block][i] - y(b*fem_length+i) );
446 double tol = rel_tol*std::abs((*sg_y)[block][i]) + abs_tol;
447 bool s = diff < tol;
448 if (!s) {
449 out << "y_expected[" << block << "][" << i << "] - "
450 << "y(" << b << "," << i << ") == "
451 << diff << " < " << tol << " : failed"
452 << std::endl;
453 }
454 success = success && s;
455 }
456 }
457
458 return success;
459 }
460
461};
462
463template <typename value_type, typename Device, typename SparseMatOps>
465 Teuchos::FancyOStream& out) {
466 typedef Stokhos::CrsMatrix<value_type,Device> matrix_type;
467 typedef typename matrix_type::values_type matrix_values_type;
468 typedef typename matrix_type::graph_type matrix_graph_type;
469 typedef Kokkos::View<value_type*,Device> vec_type;
470
471 //------------------------------
472
473 std::vector<matrix_type> matrix( setup.stoch_length ) ;
474 std::vector<vec_type> x( setup.stoch_length ) ;
475 std::vector<vec_type> y( setup.stoch_length ) ;
476 std::vector<vec_type> tmp( setup.stoch_length ) ;
477
478 for (int block=0; block<setup.stoch_length; ++block) {
479 matrix[block].graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
480 std::string("testing") , setup.fem_graph );
481
482 matrix[block].values =
483 matrix_values_type( "matrix" , setup.fem_graph_length );
484
485 x[block] = vec_type( "x" , setup.fem_length );
486 y[block] = vec_type( "y" , setup.fem_length );
487 tmp[block] = vec_type( "tmp" , setup.fem_length );
488
489 typename matrix_values_type::HostMirror hM =
490 Kokkos::create_mirror( matrix[block].values );
491
492 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < setup.fem_length ; ++iRowFEM ) {
493 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
494 const int iColFEM = setup.fem_graph[iRowFEM][iRowEntryFEM] ;
495
496 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
497 setup.fem_length , setup.stoch_length , iRowFEM , iColFEM , block );
498 }
499 }
500
501 Kokkos::deep_copy( matrix[block].values , hM );
502
503 typename vec_type::HostMirror hx =
504 Kokkos::create_mirror( x[block] );
505 typename vec_type::HostMirror hy =
506 Kokkos::create_mirror( y[block] );
507
508 for ( int i = 0 ; i < setup.fem_length ; ++i ) {
509 hx(i) = generate_vector_coefficient<value_type>(
510 setup.fem_length , setup.stoch_length , i , block );
511 hy(i) = 0.0;
512 }
513
514 Kokkos::deep_copy( x[block] , hx );
515 Kokkos::deep_copy( y[block] , hy );
516 }
517
518 // Original matrix-free multiply algorithm using a block apply
519 SparseMatOps smo;
521 setup.Cijk->k_begin();
523 setup.Cijk->k_end();
524 for (typename UnitTestSetup<Device>::Cijk_type::k_iterator k_it=k_begin;
525 k_it!=k_end; ++k_it) {
526 int nj = setup.Cijk->num_j(k_it);
527 if (nj > 0) {
528 int k = index(k_it);
530 setup.Cijk->j_begin(k_it);
532 setup.Cijk->j_end(k_it);
533 std::vector<vec_type> xx(nj), yy(nj);
534 int jdx = 0;
535 for (typename UnitTestSetup<Device>::Cijk_type::kj_iterator j_it = j_begin;
536 j_it != j_end;
537 ++j_it) {
538 int j = index(j_it);
539 xx[jdx] = x[j];
540 yy[jdx] = tmp[j];
541 jdx++;
542 }
543 Stokhos::multiply( matrix[k] , xx , yy, smo );
544 jdx = 0;
545 for (typename UnitTestSetup<Device>::Cijk_type::kj_iterator j_it = j_begin;
546 j_it != j_end; ++j_it) {
548 setup.Cijk->i_begin(j_it);
550 setup.Cijk->i_end(j_it);
551 for (typename UnitTestSetup<Device>::Cijk_type::kji_iterator i_it = i_begin;
552 i_it != i_end;
553 ++i_it) {
554 int i = index(i_it);
555 value_type c = value(i_it);
556 Stokhos::update( value_type(1.0) , y[i] , c , yy[jdx] );
557 }
558 jdx++;
559 }
560 }
561 }
562
563 std::vector<typename vec_type::HostMirror> hy(setup.stoch_length);
564 for (int i=0; i<setup.stoch_length; ++i) {
565 hy[i] = Kokkos::create_mirror( y[i] );
566 Kokkos::deep_copy( hy[i] , y[i] );
567 }
568 bool success = setup.test_original(hy, out);
569
570 return success;
571}
572
573template <typename value_type, typename Device, typename SparseMatOps>
575 Teuchos::FancyOStream& out) {
576 typedef Stokhos::CrsMatrix<value_type,Device> matrix_type;
577 typedef typename matrix_type::values_type matrix_values_type;
578 typedef typename matrix_type::graph_type matrix_graph_type;
579 typedef Kokkos::View<value_type*, Kokkos::LayoutLeft, Device, Kokkos::MemoryUnmanaged> vec_type;
580 typedef Kokkos::View<value_type**, Kokkos::LayoutLeft, Device> multi_vec_type;
581
582 //------------------------------
583
584 std::vector<matrix_type> matrix( setup.stoch_length ) ;
585 multi_vec_type x( "x", setup.fem_length, setup.stoch_length ) ;
586 multi_vec_type y( "y", setup.fem_length, setup.stoch_length ) ;
587 multi_vec_type tmp_x( "tmp_x", setup.fem_length, setup.stoch_length ) ;
588 multi_vec_type tmp_y( "tmp_y", setup.fem_length, setup.stoch_length ) ;
589
590 typename multi_vec_type::HostMirror hx = Kokkos::create_mirror( x );
591 typename multi_vec_type::HostMirror hy = Kokkos::create_mirror( y );
592
593 for (int block=0; block<setup.stoch_length; ++block) {
594 matrix[block].graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
595 std::string("testing") , setup.fem_graph );
596
597 matrix[block].values =
598 matrix_values_type( "matrix" , setup.fem_graph_length );
599
600 typename matrix_values_type::HostMirror hM =
601 Kokkos::create_mirror( matrix[block].values );
602
603 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < setup.fem_length ; ++iRowFEM ) {
604 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
605 const int iColFEM = setup.fem_graph[iRowFEM][iRowEntryFEM] ;
606
607 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
608 setup.fem_length , setup.stoch_length , iRowFEM , iColFEM , block );
609 }
610 }
611
612 Kokkos::deep_copy( matrix[block].values , hM );
613
614 for ( int i = 0 ; i < setup.fem_length ; ++i ) {
615 hx(i, block) = generate_vector_coefficient<value_type>(
616 setup.fem_length , setup.stoch_length , i , block );
617 hy(i, block) = 0.0;
618 }
619
620 }
621
622 Kokkos::deep_copy( x , hx );
623 Kokkos::deep_copy( y , hy );
624
625 // Original matrix-free multiply algorithm using a block apply
626 typedef typename UnitTestSetup<Device>::Cijk_type::k_iterator k_iterator;
627 typedef typename UnitTestSetup<Device>::Cijk_type::kj_iterator kj_iterator;
628 typedef typename UnitTestSetup<Device>::Cijk_type::kji_iterator kji_iterator;
629 SparseMatOps smo;
630 k_iterator k_begin = setup.Cijk->k_begin();
631 k_iterator k_end = setup.Cijk->k_end();
632 for (k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
633 unsigned nj = setup.Cijk->num_j(k_it);
634 if (nj > 0) {
635 int k = index(k_it);
636 kj_iterator j_begin = setup.Cijk->j_begin(k_it);
637 kj_iterator j_end = setup.Cijk->j_end(k_it);
638 std::vector<int> j_indices(nj);
639 unsigned jdx = 0;
640 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
641 int j = index(j_it);
642 vec_type xx = Kokkos::subview( x, Kokkos::ALL(), j );
643 vec_type tt = Kokkos::subview( tmp_x, Kokkos::ALL(), jdx++ );
644 Kokkos::deep_copy(tt, xx);
645 }
646 multi_vec_type tmp_x_view =
647 Kokkos::subview( tmp_x, Kokkos::ALL(),
648 std::make_pair(0u,nj));
649 multi_vec_type tmp_y_view =
650 Kokkos::subview( tmp_y, Kokkos::ALL(),
651 std::make_pair(0u,nj));
652 Stokhos::multiply( matrix[k] , tmp_x_view , tmp_y_view, smo );
653 jdx = 0;
654 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
655 vec_type tmp_y_view =
656 Kokkos::subview( tmp_y, Kokkos::ALL(), jdx++ );
657 kji_iterator i_begin = setup.Cijk->i_begin(j_it);
658 kji_iterator i_end = setup.Cijk->i_end(j_it);
659 for (kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
660 int i = index(i_it);
661 value_type c = value(i_it);
662 vec_type y_view = Kokkos::subview( y, Kokkos::ALL(), i );
663 Stokhos::update( value_type(1.0) , y_view , c , tmp_y_view );
664 }
665 }
666 }
667 }
668
669 Kokkos::deep_copy( hy , y );
670 bool success = setup.test_original(hy, out);
671
672 return success;
673}
674
675#ifdef HAVE_STOKHOS_KOKKOSLINALG
676
677template <typename value_type, typename Device>
678bool test_crs_matrix_free_kokkos(const UnitTestSetup<Device>& setup,
679 Teuchos::FancyOStream& out) {
680 typedef int ordinal_type;
681 typedef KokkosSparse::CrsMatrix<value_type,ordinal_type,Device> matrix_type;
682 typedef typename matrix_type::values_type matrix_values_type;
683 typedef typename matrix_type::StaticCrsGraphType matrix_graph_type;
684 typedef Kokkos::View<value_type*, Kokkos::LayoutLeft, Device, Kokkos::MemoryUnmanaged> vec_type;
685 typedef Kokkos::View<value_type**, Kokkos::LayoutLeft, Device> multi_vec_type;
686
687 //------------------------------
688
689 std::vector<matrix_type> matrix( setup.stoch_length ) ;
690 multi_vec_type x( "x", setup.fem_length, setup.stoch_length ) ;
691 multi_vec_type y( "y", setup.fem_length, setup.stoch_length ) ;
692 multi_vec_type tmp_x( "tmp_x", setup.fem_length, setup.stoch_length ) ;
693 multi_vec_type tmp_y( "tmp_y", setup.fem_length, setup.stoch_length ) ;
694
695 typename multi_vec_type::HostMirror hx = Kokkos::create_mirror( x );
696 typename multi_vec_type::HostMirror hy = Kokkos::create_mirror( y );
697
698 for (int block=0; block<setup.stoch_length; ++block) {
699 matrix_graph_type matrix_graph =
700 Kokkos::create_staticcrsgraph<matrix_graph_type>(
701 std::string("test crs graph"), setup.fem_graph);
702 matrix_values_type matrix_values =
703 matrix_values_type( "matrix" , setup.fem_graph_length );
704 typename matrix_values_type::HostMirror hM =
705 Kokkos::create_mirror( matrix_values );
706
707 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < setup.fem_length ; ++iRowFEM ) {
708 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
709 const int iColFEM = setup.fem_graph[iRowFEM][iRowEntryFEM] ;
710
711 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
712 setup.fem_length , setup.stoch_length , iRowFEM , iColFEM , block );
713 }
714 }
715
716 Kokkos::deep_copy( matrix_values , hM );
717 matrix[block] = matrix_type("matrix", setup.fem_length, matrix_values,
718 matrix_graph);
719
720 for ( int i = 0 ; i < setup.fem_length ; ++i ) {
721 hx(i, block) = generate_vector_coefficient<value_type>(
722 setup.fem_length , setup.stoch_length , i , block );
723 hy(i, block) = 0.0;
724 }
725
726 }
727
728 Kokkos::deep_copy( x , hx );
729 Kokkos::deep_copy( y , hy );
730
731 // Original matrix-free multiply algorithm using a block apply
732 typedef typename UnitTestSetup<Device>::Cijk_type::k_iterator k_iterator;
733 typedef typename UnitTestSetup<Device>::Cijk_type::kj_iterator kj_iterator;
734 typedef typename UnitTestSetup<Device>::Cijk_type::kji_iterator kji_iterator;
735 k_iterator k_begin = setup.Cijk->k_begin();
736 k_iterator k_end = setup.Cijk->k_end();
737 for (k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
738 int nj = setup.Cijk->num_j(k_it);
739 if (nj > 0) {
740 int k = index(k_it);
741 kj_iterator j_begin = setup.Cijk->j_begin(k_it);
742 kj_iterator j_end = setup.Cijk->j_end(k_it);
743 unsigned jdx = 0;
744 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
745 int j = index(j_it);
746 vec_type xx = Kokkos::subview( x, Kokkos::ALL(), j );
747 vec_type tt = Kokkos::subview( tmp_x, Kokkos::ALL(), jdx++ );
748 Kokkos::deep_copy(tt, xx);
749 }
750 multi_vec_type tmp_x_view =
751 Kokkos::subview( tmp_x, Kokkos::ALL(),
752 std::make_pair(0u,jdx));
753 multi_vec_type tmp_y_view =
754 Kokkos::subview( tmp_y, Kokkos::ALL(),
755 std::make_pair(0u,jdx));
756 KokkosSparse::spmv( "N", value_type(1.0), matrix[k] , tmp_x_view , value_type(0.0) , tmp_y_view );
757 jdx = 0;
758 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
759 vec_type tmp_y_view =
760 Kokkos::subview( tmp_y, Kokkos::ALL(), jdx++ );
761 kji_iterator i_begin = setup.Cijk->i_begin(j_it);
762 kji_iterator i_end = setup.Cijk->i_end(j_it);
763 for (kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
764 int i = index(i_it);
765 value_type c = value(i_it);
766 vec_type y_view = Kokkos::subview( y, Kokkos::ALL(), i );
767 //Stokhos::update( value_type(1.0) , y_view , c , tmp_y_view );
768 KokkosBlas::update(c, tmp_y_view, value_type(1.0), y_view, value_type(0.0), y_view);
769 }
770 }
771 }
772 }
773
774 Kokkos::deep_copy( hy , y );
775 bool success = setup.test_original(hy, out);
776
777 return success;
778}
779
780#endif
781
782template< typename ScalarType , class Device >
783bool
785 Teuchos::FancyOStream& out)
786{
787 typedef ScalarType value_type ;
788 typedef Kokkos::View< value_type**, Kokkos::LayoutLeft , Device > block_vector_type ;
789
790 //------------------------------
791
792 typedef Stokhos::BlockCrsMatrix< Stokhos::SymmetricDiagonalSpec< Device > , value_type , Device > matrix_type ;
793
794 typedef typename matrix_type::graph_type graph_type ;
795
796 // Convert sparse Cijk to dense for faster assembly
797 typedef typename UnitTestSetup<Device>::Cijk_type::k_iterator k_iterator;
798 typedef typename UnitTestSetup<Device>::Cijk_type::kj_iterator kj_iterator;
799 typedef typename UnitTestSetup<Device>::Cijk_type::kji_iterator kji_iterator;
800 Stokhos::Dense3Tensor<int,double> dense_cijk(setup.stoch_length);
801 for (k_iterator k_it=setup.Cijk->k_begin();
802 k_it!=setup.Cijk->k_end(); ++k_it) {
803 int k = index(k_it);
804 for (kj_iterator j_it = setup.Cijk->j_begin(k_it);
805 j_it != setup.Cijk->j_end(k_it); ++j_it) {
806 int j = index(j_it);
807 for (kji_iterator i_it = setup.Cijk->i_begin(j_it);
808 i_it != setup.Cijk->i_end(j_it); ++i_it) {
809 int i = index(i_it);
810 double c = value(i_it);
811 dense_cijk(i,j,k) = c;
812 }
813 }
814 }
815
816 //------------------------------
817 // Generate input multivector:
818
819 block_vector_type x =
820 block_vector_type( "x" , setup.stoch_length , setup.fem_length );
821 block_vector_type y =
822 block_vector_type( "y" , setup.stoch_length , setup.fem_length );
823
824 typename block_vector_type::HostMirror hx = Kokkos::create_mirror( x );
825
826 for ( int iColFEM = 0 ; iColFEM < setup.fem_length ; ++iColFEM ) {
827 for ( int iColStoch = 0 ; iColStoch < setup.stoch_length ; ++iColStoch ) {
828 hx(iColStoch,iColFEM) =
829 generate_vector_coefficient<ScalarType>(
830 setup.fem_length , setup.stoch_length , iColFEM , iColStoch );
831 }
832 }
833
834 Kokkos::deep_copy( x , hx );
835
836 //------------------------------
837 // Generate CRS matrix of blocks with symmetric diagonal storage
838
839 matrix_type matrix ;
840
841 matrix.block =
843 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
844 std::string("test crs graph") , setup.fem_graph );
845 matrix.values = block_vector_type(
846 "matrix" , matrix.block.matrix_size() , setup.fem_graph_length );
847
848 {
849 typename block_vector_type::HostMirror hM =
850 Kokkos::create_mirror( matrix.values );
851
852 for ( int iRowStoch = 0 ; iRowStoch < setup.stoch_length ; ++iRowStoch ) {
853 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < setup.fem_length ; ++iRowFEM ) {
854
855 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
856 const int iColFEM = setup.fem_graph[iRowFEM][iRowEntryFEM];
857
858 for ( int iColStoch = 0 ; iColStoch < setup.stoch_length ; ++iColStoch ) {
859
860 const size_t offset =
861 matrix.block.matrix_offset( iRowStoch , iColStoch );
862
863 ScalarType value = 0 ;
864
865 for ( int k = 0 ; k < setup.stoch_length ; ++k ) {
866 value += dense_cijk( iRowStoch , iColStoch , k ) *
867 generate_matrix_coefficient<ScalarType>(
868 setup.fem_length , setup.stoch_length , iRowFEM , iColFEM , k );
869 }
870
871 hM( offset , iEntryFEM ) = value ;
872 }
873 }
874
875 }
876 }
877
878 Kokkos::deep_copy( matrix.values , hM );
879 }
880
881 //------------------------------
882
883 Stokhos::multiply( matrix , x , y );
884
885 typename block_vector_type::HostMirror hy = Kokkos::create_mirror( y );
886 Kokkos::deep_copy( hy , y );
887
888 bool success = setup.test_commuted(hy, out);
889
890 return success;
891}
892
893template< typename ScalarType , class Device >
894bool
896 Teuchos::FancyOStream& out)
897{
898 typedef ScalarType value_type ;
899 typedef Kokkos::View< value_type* , Device > vector_type ;
900
901 //------------------------------
902
903 typedef Stokhos::CrsMatrix<value_type,Device> matrix_type;
904 typedef typename matrix_type::values_type matrix_values_type;
905 typedef typename matrix_type::graph_type matrix_graph_type;
906
907 // Build stochastic graph
908 std::vector< std::vector< int > > stoch_graph( setup.stoch_length );
909 Teuchos::RCP<Epetra_CrsGraph> cijk_graph = Stokhos::sparse3Tensor2CrsGraph(
910 *setup.basis, *setup.Cijk, *setup.globalComm);
911 for ( int i = 0 ; i < setup.stoch_length ; ++i ) {
912 int len = cijk_graph->NumGlobalIndices(i);
913 stoch_graph[i].resize(len);
914 int len2;
915 cijk_graph->ExtractGlobalRowCopy(i, len, len2, &stoch_graph[i][0]);
916 }
917
918 // Convert sparse Cijk to dense for faster assembly in debug builds
919 typedef typename UnitTestSetup<Device>::Cijk_type::k_iterator k_iterator;
920 typedef typename UnitTestSetup<Device>::Cijk_type::kj_iterator kj_iterator;
921 typedef typename UnitTestSetup<Device>::Cijk_type::kji_iterator kji_iterator;
922 Stokhos::Dense3Tensor<int,double> dense_cijk(setup.stoch_length);
923 for (k_iterator k_it=setup.Cijk->k_begin();
924 k_it!=setup.Cijk->k_end(); ++k_it) {
925 int k = index(k_it);
926 for (kj_iterator j_it = setup.Cijk->j_begin(k_it);
927 j_it != setup.Cijk->j_end(k_it); ++j_it) {
928 int j = index(j_it);
929 for (kji_iterator i_it = setup.Cijk->i_begin(j_it);
930 i_it != setup.Cijk->i_end(j_it); ++i_it) {
931 int i = index(i_it);
932 double c = value(i_it);
933 dense_cijk(i,j,k) = c;
934 }
935 }
936 }
937
938 //------------------------------
939 // Generate flattened graph with FEM outer and stochastic inner
940
941 const int flat_length = setup.fem_length * setup.stoch_length ;
942
943 std::vector< std::vector<int> > flat_graph( flat_length );
944
945 for ( int iOuterRow = 0 ; iOuterRow < setup.fem_length ; ++iOuterRow ) {
946
947 const size_t iOuterRowNZ = setup.fem_graph[iOuterRow].size();
948
949 for ( int iInnerRow = 0 ; iInnerRow < setup.stoch_length ; ++iInnerRow ) {
950
951 const size_t iInnerRowNZ = stoch_graph[ iInnerRow ].size(); ;
952 const int iFlatRowNZ = iOuterRowNZ * iInnerRowNZ ;
953 const int iFlatRow = iInnerRow + iOuterRow * setup.stoch_length ;
954
955 flat_graph[iFlatRow].resize( iFlatRowNZ );
956
957 int iFlatEntry = 0 ;
958
959 for ( size_t iOuterEntry = 0 ; iOuterEntry < iOuterRowNZ ; ++iOuterEntry ) {
960
961 const int iOuterCol = setup.fem_graph[iOuterRow][iOuterEntry];
962
963 for ( size_t iInnerEntry = 0 ; iInnerEntry < iInnerRowNZ ; ++iInnerEntry ) {
964
965 const int iInnerCol = stoch_graph[iInnerRow][iInnerEntry] ;
966 const int iFlatColumn = iInnerCol + iOuterCol * setup.stoch_length ;
967
968 flat_graph[iFlatRow][iFlatEntry] = iFlatColumn ;
969
970 ++iFlatEntry ;
971 }
972 }
973 }
974 }
975
976 //------------------------------
977
978 vector_type x = vector_type( "x" , flat_length );
979 vector_type y = vector_type( "y" , flat_length );
980
981 typename vector_type::HostMirror hx = Kokkos::create_mirror( x );
982
983 for ( int iColFEM = 0 ; iColFEM < setup.fem_length ; ++iColFEM ) {
984 for ( int iColStoch = 0 ; iColStoch < setup.stoch_length ; ++iColStoch ) {
985 hx(iColStoch + iColFEM*setup.stoch_length) =
986 generate_vector_coefficient<ScalarType>(
987 setup.fem_length , setup.stoch_length , iColFEM , iColStoch );
988 }
989 }
990
991 Kokkos::deep_copy( x , hx );
992
993 //------------------------------
994
995 matrix_type matrix ;
996
997 matrix.graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
998 std::string("testing") , flat_graph );
999
1000 const size_t flat_graph_length = matrix.graph.entries.extent(0);
1001
1002 matrix.values = matrix_values_type( "matrix" , flat_graph_length );
1003 {
1004 typename matrix_values_type::HostMirror hM =
1005 Kokkos::create_mirror( matrix.values );
1006
1007 for ( int iRow = 0 , iEntry = 0 ; iRow < flat_length ; ++iRow ) {
1008 const int iRowFEM = iRow / setup.stoch_length ;
1009 const int iRowStoch = iRow % setup.stoch_length ;
1010
1011 for ( size_t iRowEntry = 0 ; iRowEntry < flat_graph[ iRow ].size() ; ++iRowEntry , ++iEntry ) {
1012 const int iCol = flat_graph[ iRow ][ iRowEntry ];
1013 const int iColFEM = iCol / setup.stoch_length ;
1014 const int iColStoch = iCol % setup.stoch_length ;
1015
1016 double value = 0 ;
1017 for ( int k = 0 ; k < setup.stoch_length ; ++k ) {
1018 const double A_fem_k =
1019 generate_matrix_coefficient<ScalarType>(
1020 setup.fem_length , setup.stoch_length , iRowFEM, iColFEM, k );
1021 value += dense_cijk(iRowStoch,iColStoch,k) * A_fem_k ;
1022 }
1023 hM( iEntry ) = value ;
1024 }
1025 }
1026
1027 Kokkos::deep_copy( matrix.values , hM );
1028 }
1029
1030 //------------------------------
1031
1032
1033 Stokhos::multiply( matrix , x , y );
1034
1035 typename vector_type::HostMirror hy = Kokkos::create_mirror( y );
1036 Kokkos::deep_copy( hy , y );
1037
1038 bool success = setup.test_commuted_flat(hy, out);
1039 return success;
1040}
1041
1042template< typename ScalarType , class Device >
1043bool
1045 Teuchos::FancyOStream& out)
1046{
1047 typedef ScalarType value_type ;
1048 typedef Kokkos::View< value_type* , Device > vector_type ;
1049
1050 //------------------------------
1051
1052 typedef Stokhos::CrsMatrix<value_type,Device> matrix_type;
1053 typedef typename matrix_type::values_type matrix_values_type;
1054 typedef typename matrix_type::graph_type matrix_graph_type;
1055
1056 // Build stochastic graph
1057 std::vector< std::vector< int > > stoch_graph( setup.stoch_length );
1058 Teuchos::RCP<Epetra_CrsGraph> cijk_graph = Stokhos::sparse3Tensor2CrsGraph(
1059 *setup.basis, *setup.Cijk, *setup.globalComm);
1060 for ( int i = 0 ; i < setup.stoch_length ; ++i ) {
1061 int len = cijk_graph->NumGlobalIndices(i);
1062 stoch_graph[i].resize(len);
1063 int len2;
1064 cijk_graph->ExtractGlobalRowCopy(i, len, len2, &stoch_graph[i][0]);
1065 }
1066
1067 // Convert sparse Cijk to dense for faster assembly in debug builds
1068 typedef typename UnitTestSetup<Device>::Cijk_type::k_iterator k_iterator;
1069 typedef typename UnitTestSetup<Device>::Cijk_type::kj_iterator kj_iterator;
1070 typedef typename UnitTestSetup<Device>::Cijk_type::kji_iterator kji_iterator;
1071 Stokhos::Dense3Tensor<int,double> dense_cijk(setup.stoch_length);
1072 for (k_iterator k_it=setup.Cijk->k_begin();
1073 k_it!=setup.Cijk->k_end(); ++k_it) {
1074 int k = index(k_it);
1075 for (kj_iterator j_it = setup.Cijk->j_begin(k_it);
1076 j_it != setup.Cijk->j_end(k_it); ++j_it) {
1077 int j = index(j_it);
1078 for (kji_iterator i_it = setup.Cijk->i_begin(j_it);
1079 i_it != setup.Cijk->i_end(j_it); ++i_it) {
1080 int i = index(i_it);
1081 double c = value(i_it);
1082 dense_cijk(i,j,k) = c;
1083 }
1084 }
1085 }
1086
1087 //------------------------------
1088 // Generate flattened graph with stochastic outer and FEM inner
1089
1090 const size_t flat_length = setup.fem_length * setup.stoch_length ;
1091
1092 std::vector< std::vector<int> > flat_graph( flat_length );
1093
1094 for ( int iOuterRow = 0 ; iOuterRow < setup.stoch_length ; ++iOuterRow ) {
1095
1096 const size_t iOuterRowNZ = stoch_graph[iOuterRow].size();
1097
1098 for ( int iInnerRow = 0 ; iInnerRow < setup.fem_length ; ++iInnerRow ) {
1099
1100 const size_t iInnerRowNZ = setup.fem_graph[iInnerRow].size();
1101 const int iFlatRowNZ = iOuterRowNZ * iInnerRowNZ ;
1102 const int iFlatRow = iInnerRow + iOuterRow * setup.fem_length ;
1103
1104 flat_graph[iFlatRow].resize( iFlatRowNZ );
1105
1106 int iFlatEntry = 0 ;
1107
1108 for ( size_t iOuterEntry = 0 ; iOuterEntry < iOuterRowNZ ; ++iOuterEntry ) {
1109
1110 const int iOuterCol = stoch_graph[ iOuterRow ][ iOuterEntry ];
1111
1112 for ( size_t iInnerEntry = 0 ; iInnerEntry < iInnerRowNZ ; ++iInnerEntry ) {
1113
1114 const int iInnerCol = setup.fem_graph[ iInnerRow][iInnerEntry];
1115 const int iFlatColumn = iInnerCol + iOuterCol * setup.fem_length ;
1116
1117 flat_graph[iFlatRow][iFlatEntry] = iFlatColumn ;
1118 ++iFlatEntry ;
1119 }
1120 }
1121 }
1122 }
1123
1124 //------------------------------
1125
1126 vector_type x = vector_type( "x" , flat_length );
1127 vector_type y = vector_type( "y" , flat_length );
1128
1129 typename vector_type::HostMirror hx = Kokkos::create_mirror( x );
1130
1131 for ( size_t iCol = 0 ; iCol < flat_length ; ++iCol ) {
1132 const int iColStoch = iCol / setup.fem_length ;
1133 const int iColFEM = iCol % setup.fem_length ;
1134
1135 hx(iCol) = generate_vector_coefficient<ScalarType>(
1136 setup.fem_length , setup.stoch_length , iColFEM , iColStoch );
1137 }
1138
1139 Kokkos::deep_copy( x , hx );
1140
1141 //------------------------------
1142
1143 matrix_type matrix ;
1144
1145 matrix.graph = Kokkos::create_staticcrsgraph<matrix_graph_type>( std::string("testing") , flat_graph );
1146
1147 const size_t flat_graph_length = matrix.graph.entries.extent(0);
1148
1149 matrix.values = matrix_values_type( "matrix" , flat_graph_length );
1150 {
1151 typename matrix_values_type::HostMirror hM =
1152 Kokkos::create_mirror( matrix.values );
1153
1154 for ( size_t iRow = 0 , iEntry = 0 ; iRow < flat_length ; ++iRow ) {
1155 const int iRowStoch = iRow / setup.fem_length ;
1156 const int iRowFEM = iRow % setup.fem_length ;
1157
1158 for ( size_t iRowEntry = 0 ; iRowEntry < flat_graph[ iRow ].size() ; ++iRowEntry , ++iEntry ) {
1159 const int iCol = flat_graph[ iRow ][ iRowEntry ];
1160 const int iColStoch = iCol / setup.fem_length ;
1161 const int iColFEM = iCol % setup.fem_length ;
1162
1163 double value = 0 ;
1164 for ( int k = 0 ; k < setup.stoch_length ; ++k ) {
1165 const double A_fem_k =
1166 generate_matrix_coefficient<ScalarType>(
1167 setup.fem_length , setup.stoch_length ,
1168 iRowFEM , iColFEM , k );
1169 value += dense_cijk(iRowStoch,iColStoch,k) * A_fem_k ;
1170 }
1171 hM( iEntry ) = value ;
1172
1173 }
1174
1175 }
1176
1177 Kokkos::deep_copy( matrix.values , hM );
1178 }
1179
1180 Stokhos::multiply( matrix , x , y );
1181
1182 typename vector_type::HostMirror hy = Kokkos::create_mirror( y );
1183 Kokkos::deep_copy( hy , y );
1184
1185 bool success = setup.test_original_flat(hy, out);
1186 return success;
1187}
1188
1189template< typename ScalarType , typename TensorType, class Device >
1192 Teuchos::FancyOStream& out,
1193 const Teuchos::ParameterList& params = Teuchos::ParameterList()) {
1194 typedef ScalarType value_type ;
1195 typedef Kokkos::View< value_type** , Kokkos::LayoutLeft ,
1196 Device > block_vector_type ;
1197
1199
1201 typedef typename matrix_type::graph_type graph_type ;
1202
1203 //------------------------------
1204 // Generate input multivector:
1205
1206 block_vector_type x =
1207 block_vector_type( "x" , setup.stoch_length_aligned , setup.fem_length );
1208 block_vector_type y =
1209 block_vector_type( "y" , setup.stoch_length_aligned , setup.fem_length );
1210
1211 typename block_vector_type::HostMirror hx =
1213
1214 for ( int iColFEM = 0 ; iColFEM < setup.fem_length ; ++iColFEM ) {
1215 for ( int iColStoch = 0 ; iColStoch < setup.stoch_length ; ++iColStoch ) {
1216 hx(setup.perm[iColStoch],iColFEM) =
1217 generate_vector_coefficient<ScalarType>(
1218 setup.fem_length , setup.stoch_length , iColFEM , iColStoch );
1219 }
1220 }
1221
1222 Kokkos::deep_copy( x , hx );
1223
1224 //------------------------------
1225
1226 matrix_type matrix ;
1227
1228 matrix.block =
1229 Stokhos::create_stochastic_product_tensor< TensorType >( *setup.basis,
1230 *setup.Cijk,
1231 params);
1232
1233 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1234 std::string("test crs graph") , setup.fem_graph );
1235
1236 matrix.values = block_vector_type(
1237 "matrix" , setup.stoch_length_aligned , setup.fem_graph_length );
1238
1239 typename block_vector_type::HostMirror hM =
1240 Kokkos::create_mirror( matrix.values );
1241
1242 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < setup.fem_length ; ++iRowFEM ) {
1243 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1244 const int iColFEM = setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1245
1246 for ( int k = 0 ; k < setup.stoch_length ; ++k ) {
1247 hM(setup.perm[k],iEntryFEM) =
1248 generate_matrix_coefficient<ScalarType>(
1249 setup.fem_length , setup.stoch_length , iRowFEM , iColFEM , k );
1250 //hM(k,iEntryFEM) = 1.0;
1251 }
1252 }
1253 }
1254
1255 Kokkos::deep_copy( matrix.values , hM );
1256
1257 //------------------------------
1258
1259 Stokhos::multiply( matrix , x , y );
1260
1261 typename block_vector_type::HostMirror hy = Kokkos::create_mirror( y );
1262 Kokkos::deep_copy( hy , y );
1263
1264 bool success = setup.test_commuted_perm(hy, out);
1265 //bool success = true;
1266 return success;
1267}
1268
1269template< typename ScalarType , class Device , int BlockSize >
1271 Teuchos::FancyOStream& out,
1272 const bool symmetric) {
1273 typedef ScalarType value_type ;
1274 typedef Kokkos::View< value_type** , Kokkos::LayoutLeft ,
1275 Device > block_vector_type ;
1278
1280 typedef typename matrix_type::graph_type graph_type ;
1281
1282 // Build tensor
1283 matrix_type matrix ;
1284
1285 Teuchos::ParameterList params;
1286 params.set("Symmetric",symmetric);
1287 matrix.block =
1288 Stokhos::create_stochastic_product_tensor< TensorType >( *setup.basis,
1289 *setup.Cijk,
1290 params );
1291 int aligned_stoch_length = matrix.block.tensor().aligned_dimension();
1292
1293 //------------------------------
1294 // Generate input multivector:
1295
1296 block_vector_type x =
1297 block_vector_type( "x" , aligned_stoch_length , setup.fem_length );
1298 block_vector_type y =
1299 block_vector_type( "y" , aligned_stoch_length , setup.fem_length );
1300
1301 typename block_vector_type::HostMirror hx =
1303
1304 for ( int iColFEM = 0 ; iColFEM < setup.fem_length ; ++iColFEM ) {
1305 for ( int iColStoch = 0 ; iColStoch < setup.stoch_length ; ++iColStoch ) {
1306 hx(iColStoch,iColFEM) =
1307 generate_vector_coefficient<ScalarType>(
1308 setup.fem_length , setup.stoch_length , iColFEM , iColStoch );
1309 //hx(iColStoch,iColFEM) = 1.0;
1310 }
1311 }
1312
1313 Kokkos::deep_copy( x , hx );
1314
1315 //------------------------------
1316
1317 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1318 std::string("test crs graph") , setup.fem_graph );
1319
1320 matrix.values = block_vector_type(
1321 "matrix" , aligned_stoch_length , setup.fem_graph_length );
1322
1323 typename block_vector_type::HostMirror hM =
1324 Kokkos::create_mirror( matrix.values );
1325
1326 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < setup.fem_length ; ++iRowFEM ) {
1327 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1328 const int iColFEM = setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1329
1330 for ( int k = 0 ; k < setup.stoch_length ; ++k ) {
1331 hM(k,iEntryFEM) =
1332 generate_matrix_coefficient<ScalarType>(
1333 setup.fem_length , setup.stoch_length , iRowFEM , iColFEM , k );
1334 //hM(k,iEntryFEM) = 1.0;
1335 }
1336 }
1337 }
1338
1339 Kokkos::deep_copy( matrix.values , hM );
1340
1341 //------------------------------
1342
1343 Stokhos::multiply( matrix , x , y );
1344
1345 typename block_vector_type::HostMirror hy = Kokkos::create_mirror( y );
1346 Kokkos::deep_copy( hy , y );
1347
1348 bool success = setup.test_commuted(hy, out);
1349 //bool success = true;
1350 return success;
1351}
1352
1353template< typename ScalarType , class Device >
1355 Teuchos::FancyOStream& out) {
1356 typedef ScalarType value_type ;
1357 typedef int ordinal_type;
1358 typedef Kokkos::View< value_type** , Kokkos::LayoutLeft ,
1359 Device > block_vector_type ;
1362
1364 typedef typename matrix_type::graph_type graph_type ;
1365
1366 //------------------------------
1367 // Generate input multivector:
1368
1369 block_vector_type x =
1370 block_vector_type( "x" , setup.stoch_length , setup.fem_length );
1371 block_vector_type y =
1372 block_vector_type( "y" , setup.stoch_length , setup.fem_length );
1373
1374 typename block_vector_type::HostMirror hx =
1376
1377 for ( int iColFEM = 0 ; iColFEM < setup.fem_length ; ++iColFEM ) {
1378 for ( int iColStoch = 0 ; iColStoch < setup.stoch_length ; ++iColStoch ) {
1379 hx(iColStoch,iColFEM) =
1380 generate_vector_coefficient<ScalarType>(
1381 setup.fem_length , setup.stoch_length , iColFEM , iColStoch );
1382 }
1383 }
1384
1385 Kokkos::deep_copy( x , hx );
1386
1387 //------------------------------
1388
1389 matrix_type matrix ;
1390
1391 /*
1392 typedef UnitTestSetup<Device>::abstract_basis_type abstract_basis_type;
1393 typedef UnitTestSetup<Device>::basis_type basis_type;
1394 typedef Stokhos::LexographicLess< Stokhos::MultiIndex<int> > less_type;
1395 typedef Stokhos::TotalOrderBasis<ordinal_type,value_type,less_type> product_basis_type;
1396 Teuchos::Array< Teuchos::RCP<const abstract_basis_type> > bases(setup.d);
1397 for (int i=0; i<setup.d; i++)
1398 bases[i] = rcp(new basis_type(setup.p,true));
1399 product_basis_type basis(bases, 1e-12);
1400 */
1401 const bool symmetric = true;
1402 Teuchos::RCP< Stokhos::LTBSparse3Tensor<ordinal_type, value_type> > Cijk =
1404
1405 matrix.block =
1406 Stokhos::create_stochastic_product_tensor< TensorType >( *setup.basis,
1407 *Cijk );
1408
1409 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1410 std::string("test crs graph") , setup.fem_graph );
1411
1412 matrix.values = block_vector_type(
1413 "matrix" , setup.stoch_length , setup.fem_graph_length );
1414
1415 typename block_vector_type::HostMirror hM =
1416 Kokkos::create_mirror( matrix.values );
1417
1418 for ( int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM < setup.fem_length ; ++iRowFEM ) {
1419 for ( size_t iRowEntryFEM = 0 ; iRowEntryFEM < setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1420 const int iColFEM = setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1421
1422 for ( int k = 0 ; k < setup.stoch_length ; ++k ) {
1423 hM(k,iEntryFEM) =
1424 generate_matrix_coefficient<ScalarType>(
1425 setup.fem_length , setup.stoch_length , iRowFEM , iColFEM , k );
1426 }
1427 }
1428 }
1429
1430 Kokkos::deep_copy( matrix.values , hM );
1431
1432 //------------------------------
1433
1434 Stokhos::multiply( matrix , x , y );
1435
1436 typename block_vector_type::HostMirror hy = Kokkos::create_mirror( y );
1437 Kokkos::deep_copy( hy , y );
1438
1439 bool success = setup.test_commuted(hy, out);
1440 return success;
1441}
1442
1443}
1444
1445#endif
Copy
UnitTestSetup< Kokkos::Cuda > setup
CRS matrix of dense blocks.
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
Data structure storing a dense 3-tensor C(i,j,k).
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
Legendre polynomial basis.
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
A comparison functor implementing a strict weak ordering based lexographic ordering.
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
An Epetra operator representing the block stochastic Galerkin operator.
A multidimensional index.
Abstract base class for 1-D orthogonal polynomials.
A container class for products of Epetra_Vector's.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
Bases defined by combinatorial product of polynomial bases.
Symmetric diagonal storage for a dense matrix.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Sacado::MP::Vector< storage_type > vec_type
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< XD, XP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< YD, YP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< ZD, ZP... > >::value >::type update(const typename Kokkos::View< XD, XP... >::array_type::non_const_value_type &alpha, const Kokkos::View< XD, XP... > &x, const typename Kokkos::View< YD, YP... >::array_type::non_const_value_type &beta, const Kokkos::View< YD, YP... > &y, const typename Kokkos::View< ZD, ZP... >::array_type::non_const_value_type &gamma, const Kokkos::View< ZD, ZP... > &z)
bool test_crs_flat_original(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
bool test_crs_product_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out, const Teuchos::ParameterList &params=Teuchos::ParameterList())
bool test_crs_matrix_free_view(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
bool test_crs_dense_block(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
bool test_crs_flat_commuted(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
bool test_linear_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out, const bool symmetric)
bool test_crs_matrix_free(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
bool test_lexo_block_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< InputType, InputP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< OutputType, OutputP... > >::value >::type spmv(const char mode[], const AlphaType &a, const MatrixType &A, const Kokkos::View< InputType, InputP... > &x, const BetaType &b, const Kokkos::View< OutputType, OutputP... > &y, const RANK_ONE)
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTBBlockLeaf(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Teuchos::RCP< Epetra_CrsGraph > sparse3Tensor2CrsGraph(const Stokhos::OrthogPolyBasis< ordinal_type, value_type > &basis, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, const Epetra_Comm &comm)
Build an Epetra_CrsGraph from a sparse 3 tensor.
void update(const ValueType &alpha, VectorType &x, const ValueType &beta, const VectorType &y)
void multiply(const CrsMatrix< MatrixValue, Device, Layout > &A, const InputMultiVectorType &x, OutputMultiVectorType &y, const std::vector< OrdinalType > &col_indices, SingleColumnMultivectorMultiply)
Stokhos::TotalOrderBasis< int, value_type, less_type > product_basis_type
bool test_original(const multi_vec_type &y, Teuchos::FancyOStream &out) const
bool test_commuted_perm(const vec_type &y, Teuchos::FancyOStream &out) const
bool test_commuted_flat(const vec_type &y, Teuchos::FancyOStream &out) const
Stokhos::Sparse3Tensor< int, value_type > Cijk_type
bool test_original(const std::vector< vec_type > &y, Teuchos::FancyOStream &out) const
bool test_original_flat(const vec_type &y, Teuchos::FancyOStream &out) const
bool test_commuted(const vec_type &y, Teuchos::FancyOStream &out) const
Stokhos::LexographicLess< Stokhos::MultiIndex< int > > less_type
Stokhos::LegendreBasis< int, value_type > basis_type
Stokhos::OneDOrthogPolyBasis< int, value_type > abstract_basis_type
Bi-directional iterator for traversing a sparse array.