Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_Tpetra_Utilities_UQ_PCE.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_TPETRA_UTILITIES_UQ_PCE_HPP
43#define STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP
44
47#include "Tpetra_Map.hpp"
48#include "Tpetra_MultiVector.hpp"
49#include "Tpetra_CrsGraph.hpp"
50#include "Tpetra_CrsMatrix.hpp"
51
52namespace Stokhos {
53
54
55 // Build a CRS graph from a sparse Cijk tensor
56 template <typename LocalOrdinal, typename GlobalOrdinal, typename Device,
57 typename CijkType>
58 Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,
59 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
60 create_cijk_crs_graph(const CijkType& cijk_dev,
61 const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
62 const size_t matrix_pce_size) {
63 using Teuchos::RCP;
64 using Teuchos::arrayView;
65
66 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
67 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
68 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
69
70 // Code below accesses cijk entries on the host, so make sure it is
71 // accessible there
72 auto cijk = create_mirror_view(cijk_dev);
73 deep_copy(cijk, cijk_dev);
74
75 const size_t pce_sz = cijk.dimension();
76 RCP<const Map> map =
77 Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(pce_sz, comm);
78 RCP<Graph> graph;
79 if (matrix_pce_size == 1) {
80 graph = Tpetra::createCrsGraph(map, 1);
81 // Mean-based case -- graph is diagonal
82 for (size_t i=0; i<pce_sz; ++i) {
83 const GlobalOrdinal row = i;
84 graph->insertGlobalIndices(row, arrayView(&row, 1));
85 }
86 }
87 else {
88 // General case
89
90 // Get max num entries
91 size_t max_num_entry = 0;
92 for (size_t i=0; i<pce_sz; ++i) {
93 const size_t num_entry = cijk.num_entry(i);
94 max_num_entry = (num_entry > max_num_entry) ? num_entry : max_num_entry;
95 }
96 max_num_entry *= 2; // 1 entry each for j, k coord
97 graph = Tpetra::createCrsGraph(map, max_num_entry);
98
99 for (size_t i=0; i<pce_sz; ++i) {
100 const GlobalOrdinal row = i;
101 const size_t num_entry = cijk.num_entry(i);
102 const size_t entry_beg = cijk.entry_begin(i);
103 const size_t entry_end = entry_beg + num_entry;
104 for (size_t entry = entry_beg; entry < entry_end; ++entry) {
105 const GlobalOrdinal j = cijk.coord(entry,0);
106 const GlobalOrdinal k = cijk.coord(entry,1);
107 graph->insertGlobalIndices(row, arrayView(&j, 1));
108 graph->insertGlobalIndices(row, arrayView(&k, 1));
109 }
110 }
111 }
112 graph->fillComplete();
113 return graph;
114 }
115
116 // Create a flattened graph for a graph from a matrix with the
117 // UQ::PCE scalar type
118 // If flat_domain_map and/or flat_range_map are null, they will be computed,
119 // otherwise they will be used as-is.
120 template <typename LocalOrdinal, typename GlobalOrdinal, typename Device,
121 typename CijkType>
122 Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,
123 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
125 const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& graph,
126 const CijkType& cijk,
127 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_domain_map,
128 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_range_map,
129 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& cijk_graph,
130 const size_t matrix_pce_size) {
131 using Teuchos::ArrayView;
132 using Teuchos::ArrayRCP;
133 using Teuchos::Array;
134 using Teuchos::RCP;
135 using Teuchos::rcp;
136
137 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
138 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
139 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
140
141 const LocalOrdinal block_size = cijk.dimension();
142
143 // Build domain map if necessary
144 if (flat_domain_map == Teuchos::null)
145 flat_domain_map = create_flat_map(*(graph.getDomainMap()), block_size);
146
147 // Build range map if necessary
148 if (flat_range_map == Teuchos::null)
149 flat_range_map = create_flat_map(*(graph.getRangeMap()), block_size);
150
151 // Build column map
152 RCP<const Map> flat_col_map =
153 create_flat_map(*(graph.getColMap()), block_size);
154
155 // Build row map if necessary
156 // Check if range_map == row_map, then we can use flat_range_map
157 // as the flattened row map
158 RCP<const Map> flat_row_map;
159 if (graph.getRangeMap() == graph.getRowMap())
160 flat_row_map = flat_range_map;
161 else
162 flat_row_map = create_flat_map(*(graph.getRowMap()), block_size);
163
164 // Build Cijk graph if necessary
165 if (cijk_graph == Teuchos::null)
166 cijk_graph = create_cijk_crs_graph<LocalOrdinal,GlobalOrdinal,Device>(
167 cijk,
168 flat_domain_map->getComm(),
169 matrix_pce_size);
170
171 // Build flattened graph that is the Kronecker product of the given
172 // graph and cijk_graph
173
174 // Loop over outer rows
175 typename Graph::local_inds_host_view_type outer_cols;
176 typename Graph::local_inds_host_view_type inner_cols;
177 size_t max_num_row_entries = graph.getLocalMaxNumRowEntries()*block_size;
178 Array<LocalOrdinal> flat_col_indices;
179 flat_col_indices.reserve(max_num_row_entries);
180 RCP<Graph> flat_graph = rcp(new Graph(flat_row_map, flat_col_map, max_num_row_entries));
181 const LocalOrdinal num_outer_rows = graph.getLocalNumRows();
182 for (LocalOrdinal outer_row=0; outer_row < num_outer_rows; outer_row++) {
183
184 // Get outer columns for this outer row
185 Kokkos::fence();
186 graph.getLocalRowView(outer_row, outer_cols);
187 const LocalOrdinal num_outer_cols = outer_cols.size();
188
189 // Loop over inner rows
190 for (LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
191
192 // Compute flat row index
193 const LocalOrdinal flat_row = outer_row*block_size + inner_row;
194
195 // Get inner columns for this inner row
196 Kokkos::fence();
197 cijk_graph->getLocalRowView(inner_row, inner_cols);
198 const LocalOrdinal num_inner_cols = inner_cols.size();
199
200 // Create space to store all column indices for this flat row
201 flat_col_indices.resize(0);
202
203 // Loop over outer cols
204 for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
205 ++outer_entry) {
206 const LocalOrdinal outer_col = outer_cols[outer_entry];
207
208 // Loop over inner cols
209 for (LocalOrdinal inner_entry=0; inner_entry<num_inner_cols;
210 ++inner_entry) {
211 const LocalOrdinal inner_col = inner_cols[inner_entry];
212
213 // Compute and store flat column index
214 const LocalOrdinal flat_col = outer_col*block_size + inner_col;
215 flat_col_indices.push_back(flat_col);
216 }
217
218 }
219
220 // Insert all indices for this flat row
221 flat_graph->insertLocalIndices(flat_row, flat_col_indices());
222
223 }
224
225 }
226 flat_graph->fillComplete(flat_domain_map, flat_range_map);
227
228 return flat_graph;
229 }
230
231 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
232 // returned vector is a view of the original
233 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
234 typename Device>
235 Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
237 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
239 const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
241 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
242 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
243 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
244 using Teuchos::RCP;
245 using Teuchos::rcp;
246
247 typedef typename Storage::value_type BaseScalar;
248 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
249 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
250 typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
251
252 // Have to do a nasty const-cast because getLocalViewDevice(ReadWrite) is a
253 // non-const method, yet getLocalViewDevice(ReadOnly) returns a const-view
254 // (i.e., with a constant scalar type), and there is no way to make a
255 // MultiVector out of it!
256 typedef Tpetra::MultiVector<Sacado::UQ::PCE<Storage>, LocalOrdinal,GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode<Device> > mv_type;
257 mv_type& vec_nc = const_cast<mv_type&>(vec);
258
259 // Create flattenend view using special reshaping view assignment operator
260 flat_view_type flat_vals = vec_nc.getLocalViewDevice(Tpetra::Access::ReadWrite);
261
262 // Create flat vector
263 RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
264
265 return flat_vec;
266 }
267
268 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
269 // returned vector is a view of the original
270 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
271 typename Device>
272 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
274 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
276 Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
278 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
279 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
280 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
281 using Teuchos::RCP;
282 using Teuchos::rcp;
283
284 typedef typename Storage::value_type BaseScalar;
285 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
286 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
287 typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
288
289 // Create flattenend view using special reshaping view assignment operator
290 flat_view_type flat_vals = vec.getLocalViewDevice(Tpetra::Access::ReadWrite);
291
292 // Create flat vector
293 RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
294
295 return flat_vec;
296 }
297
298 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
299 // returned vector is a view of the original. This version creates the
300 // map if necessary
301 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
302 typename Device>
303 Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
305 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
307 const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
309 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
310 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
311 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
312 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
313 if (flat_map == Teuchos::null) {
314 const LocalOrdinal pce_size =
315 Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
316 flat_map = create_flat_map(*(vec.getMap()), pce_size);
317 }
318 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
319 return create_flat_vector_view(vec, const_flat_map);
320 }
321
322 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
323 // returned vector is a view of the original. This version creates the
324 // map if necessary
325 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
326 typename Device>
327 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
329 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
331 Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
333 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
334 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
335 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
336 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
337 if (flat_map == Teuchos::null) {
338 const LocalOrdinal pce_size =
339 Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
340 flat_map = create_flat_map(*(vec.getMap()), pce_size);
341 }
342 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
343 return create_flat_vector_view(vec, const_flat_map);
344 }
345
346 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
347 // returned vector is a view of the original
348 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
349 typename Device>
350 Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
352 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
354 const Tpetra::Vector<Sacado::UQ::PCE<Storage>,
356 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec_const,
357 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
358 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
359 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
360 const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec_const;
361 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
362 return flat_mv->getVector(0);
363 }
364
365 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
366 // returned vector is a view of the original. This version creates the
367 // map if necessary
368 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
369 typename Device>
370 Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
372 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
374 const Tpetra::Vector<Sacado::UQ::PCE<Storage>,
376 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
377 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
378 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
379 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
380 if (flat_map == Teuchos::null) {
381 const LocalOrdinal pce_size =
382 Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
383 flat_map = create_flat_map(*(vec.getMap()), pce_size);
384 }
385 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
386 return create_flat_vector_view(vec, const_flat_map);
387 }
388
389 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
390 // returned vector is a view of the original
391 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
392 typename Device>
393 Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
395 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
397 Tpetra::Vector<Sacado::UQ::PCE<Storage>,
399 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
400 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
401 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
402 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
403 Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec;
404 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
405 return flat_mv->getVectorNonConst(0);
406 }
407
408 // Create a flattened vector by unrolling the UQ::PCE scalar type. The
409 // returned vector is a view of the original. This version creates the
410 // map if necessary
411 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
412 typename Device>
413 Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
415 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
417 Tpetra::Vector<Sacado::UQ::PCE<Storage>,
419 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
420 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
421 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
422 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
423 if (flat_map == Teuchos::null) {
424 const LocalOrdinal pce_size =
425 Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
426 flat_map = create_flat_map(*(vec.getMap()), pce_size);
427 }
428 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
429 return create_flat_vector_view(vec, const_flat_map);
430 }
431
432 // Create a flattened matrix by unrolling the UQ::PCE scalar type. The
433 // returned matrix is NOT a view of the original (and can't be)
434 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
435 typename Device, typename CijkType>
436 Teuchos::RCP< Tpetra::CrsMatrix<typename Storage::value_type,
438 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
440 const Tpetra::CrsMatrix<Sacado::UQ::PCE<Storage>,
441 LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& mat,
442 const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_graph,
443 const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& cijk_graph,
444 const CijkType& cijk_dev) {
445 using Teuchos::ArrayView;
446 using Teuchos::Array;
447 using Teuchos::RCP;
448 using Teuchos::rcp;
449
450 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
452 typedef typename Storage::value_type BaseScalar;
453 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Matrix;
454 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
455
456 // Code below accesses cijk entries on the host, so make sure it is
457 // accessible there
458 auto cijk = create_mirror_view(cijk_dev);
459 deep_copy(cijk, cijk_dev);
460
461 const LocalOrdinal block_size = cijk.dimension();
462 const LocalOrdinal matrix_pce_size =
463 Kokkos::dimension_scalar(mat.getLocalMatrixDevice().values);
464
465 // Create flat matrix
466 RCP<FlatMatrix> flat_mat = rcp(new FlatMatrix(flat_graph));
467
468 // Fill flat matrix
469 typename Matrix::values_host_view_type outer_values;
470 typename Matrix::local_inds_host_view_type outer_cols;
471 typename Matrix::local_inds_host_view_type inner_cols;
472 typename Matrix::local_inds_host_view_type flat_cols;
473 Array<BaseScalar> flat_values;
474 flat_values.reserve(flat_graph->getLocalMaxNumRowEntries());
475 const LocalOrdinal num_outer_rows = mat.getLocalNumRows();
476 for (LocalOrdinal outer_row=0; outer_row < num_outer_rows; outer_row++) {
477
478 // Get outer columns and values for this outer row
479 mat.getLocalRowView(outer_row, outer_cols, outer_values);
480 const LocalOrdinal num_outer_cols = outer_cols.size();
481
482 // Loop over inner rows
483 for (LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
484
485 // Compute flat row index
486 const LocalOrdinal flat_row = outer_row*block_size + inner_row;
487
488 // Get cijk column indices for this row
489 cijk_graph->getLocalRowView(inner_row, inner_cols);
490 const LocalOrdinal num_inner_cols = inner_cols.size();
491 ArrayView<const LocalOrdinal> inner_cols_av =
492 Kokkos::Compat::getConstArrayView(inner_cols);
493
494 // Create space to store all values for this flat row
495 const LocalOrdinal num_flat_indices = num_outer_cols*num_inner_cols;
496 //flat_values.resize(num_flat_indices);
497 flat_values.assign(num_flat_indices, BaseScalar(0));
498
499 if (matrix_pce_size == 1) {
500 // Mean-based case
501
502 // Loop over outer cols
503 for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
504 ++outer_entry) {
505
506 // Extract mean PCE entry for each outer column
507 flat_values[outer_entry] = outer_values[outer_entry].coeff(0);
508
509 }
510
511 }
512 else {
513
514 // Loop over cijk non-zeros for this inner row
515 const size_t num_entry = cijk.num_entry(inner_row);
516 const size_t entry_beg = cijk.entry_begin(inner_row);
517 const size_t entry_end = entry_beg + num_entry;
518 for (size_t entry = entry_beg; entry < entry_end; ++entry) {
519 const LocalOrdinal j = cijk.coord(entry,0);
520 const LocalOrdinal k = cijk.coord(entry,1);
521 const BaseScalar c = cijk.value(entry);
522
523 // Find column offset for j
524 typedef typename ArrayView<const LocalOrdinal>::iterator iterator;
525 iterator ptr_j =
526 std::find(inner_cols_av.begin(), inner_cols_av.end(), j);
527 iterator ptr_k =
528 std::find(inner_cols_av.begin(), inner_cols_av.end(), k);
529 const LocalOrdinal j_offset = ptr_j - inner_cols_av.begin();
530 const LocalOrdinal k_offset = ptr_k - inner_cols_av.begin();
531
532 // Loop over outer cols
533 for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
534 ++outer_entry) {
535
536 // Add contributions for each outer column
537 flat_values[outer_entry*num_inner_cols + j_offset] +=
538 c*outer_values[outer_entry].coeff(k);
539 flat_values[outer_entry*num_inner_cols + k_offset] +=
540 c*outer_values[outer_entry].coeff(j);
541
542 }
543
544 }
545
546 }
547
548 // Set values in flat matrix
549 flat_graph->getLocalRowView(flat_row, flat_cols);
550 flat_mat->replaceLocalValues(flat_row, Kokkos::Compat::getConstArrayView(flat_cols), flat_values());
551
552 }
553
554 }
555 flat_mat->fillComplete(flat_graph->getDomainMap(),
556 flat_graph->getRangeMap());
557
558 return flat_mat;
559 }
560
561
562} // namespace Stokhos
563
564#endif // STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP
KokkosClassic::DefaultNode::DefaultNodeType Node
Stokhos::StandardStorage< int, double > Storage
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< View< T, P... > >::value, unsigned >::type dimension_scalar(const View< T, P... > &view)
Top-level namespace for Stokhos classes and functions.
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_pce_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &graph, const CijkType &cijk, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_range_map, Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const size_t matrix_pce_size)
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_map)
void deep_copy(const CrsProductTensor< ValueType, DstDevice, DstMemory > &dst, const CrsProductTensor< ValueType, SrcDevice, SrcMemory > &src)
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_cijk_crs_graph(const CijkType &cijk_dev, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const size_t matrix_pce_size)
CrsProductTensor< ValueType, Device, Memory >::HostMirror create_mirror_view(const CrsProductTensor< ValueType, Device, Memory > &src)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const CijkType &cijk_dev)
Teuchos::RCP< Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > create_flat_map(const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map, const LocalOrdinal block_size)