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