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