42 #ifndef STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP 43 #define STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP 46 #include "Tpetra_Map.hpp" 47 #include "Tpetra_MultiVector.hpp" 48 #include "Tpetra_CrsGraph.hpp" 49 #include "Tpetra_CrsMatrix.hpp" 55 template <
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
56 Teuchos::RCP< Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
59 using Tpetra::global_size_t;
60 using Teuchos::ArrayView;
65 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
68 const global_size_t num_global_entries = map.getGlobalNumElements();
69 const size_t num_local_entries = map.getNodeNumElements();
71 ArrayView<const GlobalOrdinal> element_list =
72 map.getNodeElementList();
75 const global_size_t flat_num_global_entries = num_global_entries*block_size;
76 const size_t flat_num_local_entries = num_local_entries * block_size;
78 Array<GlobalOrdinal> flat_element_list(flat_num_local_entries);
79 for (
size_t i=0; i<num_local_entries; ++i)
81 flat_element_list[i*block_size+
j] = element_list[i]*block_size+
j;
85 rcp(
new Map(flat_num_global_entries, flat_element_list(),
86 flat_index_base, map.getComm(), map.getNode()));
95 template <
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
96 Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >
98 const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>& graph,
99 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_domain_map,
100 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_range_map,
102 using Teuchos::ArrayRCP;
106 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
107 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
110 if (flat_domain_map == Teuchos::null)
111 flat_domain_map =
create_flat_map(*(graph.getDomainMap()), block_size);
114 if (flat_range_map == Teuchos::null)
118 RCP<const Map> flat_col_map =
124 RCP<const Map> flat_row_map;
125 if (graph.getRangeMap() == graph.getRowMap())
126 flat_row_map = flat_range_map;
131 ArrayRCP<const size_t> row_offsets = graph.getNodeRowPtrs();
132 ArrayRCP<const LocalOrdinal> col_indices = graph.getNodePackedIndices();
133 const size_t num_row = graph.getNodeNumRows();
134 const size_t num_col_indices = col_indices.size();
135 ArrayRCP<size_t> flat_row_offsets(num_row*block_size+1);
136 ArrayRCP<LocalOrdinal> flat_col_indices(num_col_indices * block_size);
137 for (
size_t row=0; row<num_row; ++row) {
138 const size_t row_beg = row_offsets[row];
139 const size_t row_end = row_offsets[row+1];
140 const size_t num_col = row_end - row_beg;
142 const size_t flat_row = row*block_size +
j;
143 const size_t flat_row_beg = row_beg*block_size +
j*num_col;
144 flat_row_offsets[flat_row] = flat_row_beg;
145 for (
size_t entry=0; entry<num_col; ++entry) {
148 flat_col_indices[flat_row_beg+entry] = flat_col;
152 flat_row_offsets[num_row*block_size] = num_col_indices*block_size;
155 RCP<Graph> flat_graph =
156 rcp(
new Graph(flat_row_map, flat_col_map,
157 flat_row_offsets, flat_col_indices));
158 flat_graph->fillComplete(flat_domain_map, flat_range_map);
163 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR) 169 template <
typename LocalOrdinal,
typename GlobalOrdinal,
typename Device>
171 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
174 Teuchos::RCP<
const Tpetra::Map<
LocalOrdinal,
GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_domain_map,
175 Teuchos::RCP<
const Tpetra::Map<
LocalOrdinal,
GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_range_map,
177 using Teuchos::ArrayRCP;
181 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
182 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
183 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
184 typedef typename Graph::local_graph_type::row_map_type::non_const_type RowPtrs;
185 typedef typename Graph::local_graph_type::entries_type::non_const_type LocalIndices;
188 if (flat_domain_map == Teuchos::null)
189 flat_domain_map =
create_flat_map(*(graph.getDomainMap()), block_size);
192 if (flat_range_map == Teuchos::null)
196 RCP<const Map> flat_col_map =
202 RCP<const Map> flat_row_map;
203 if (graph.getRangeMap() == graph.getRowMap())
204 flat_row_map = flat_range_map;
209 ArrayRCP<const size_t> row_offsets = graph.getNodeRowPtrs();
210 ArrayRCP<const LocalOrdinal> col_indices = graph.getNodePackedIndices();
211 const size_t num_row = graph.getNodeNumRows();
212 const size_t num_col_indices = col_indices.size();
213 RowPtrs flat_row_offsets(
"row_ptrs", num_row*block_size+1);
214 LocalIndices flat_col_indices(
"col_indices", num_col_indices * block_size);
215 for (
size_t row=0; row<num_row; ++row) {
216 const size_t row_beg = row_offsets[row];
217 const size_t row_end = row_offsets[row+1];
218 const size_t num_col = row_end - row_beg;
220 const size_t flat_row = row*block_size +
j;
221 const size_t flat_row_beg = row_beg*block_size +
j*num_col;
222 flat_row_offsets[flat_row] = flat_row_beg;
223 for (
size_t entry=0; entry<num_col; ++entry) {
226 flat_col_indices[flat_row_beg+entry] = flat_col;
230 flat_row_offsets[num_row*block_size] = num_col_indices*block_size;
233 RCP<Graph> flat_graph =
234 rcp(
new Graph(flat_row_map, flat_col_map,
235 flat_row_offsets, flat_col_indices));
236 flat_graph->fillComplete(flat_domain_map, flat_range_map);
252 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
253 using Teuchos::ArrayRCP;
259 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Vector;
260 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
267 Vector& vec =
const_cast<Vector&
>(vec_const);
270 ArrayRCP<Scalar> vec_vals = vec.get1dViewNonConst();
271 const size_t vec_size = vec_vals.size();
274 BaseScalar *flat_vec_ptr =
275 reinterpret_cast<BaseScalar*
>(vec_vals.getRawPtr());
276 ArrayRCP<BaseScalar> flat_vec_vals =
277 Teuchos::arcp(flat_vec_ptr, 0, vec_size*mp_size,
false);
280 const size_t stride = vec.getStride();
281 const size_t flat_stride = stride * mp_size;
282 const size_t num_vecs = vec.getNumVectors();
283 RCP<FlatVector> flat_vec =
284 rcp(
new FlatVector(flat_map, flat_vec_vals, flat_stride, num_vecs));
299 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
300 if (flat_map == Teuchos::null) {
304 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
317 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
318 using Teuchos::ArrayRCP;
324 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
330 ArrayRCP<Scalar> vec_vals = vec.get1dViewNonConst();
331 const size_t vec_size = vec_vals.size();
334 BaseScalar *flat_vec_ptr =
335 reinterpret_cast<BaseScalar*
>(vec_vals.getRawPtr());
336 ArrayRCP<BaseScalar> flat_vec_vals =
337 Teuchos::arcp(flat_vec_ptr, 0, vec_size*mp_size,
false);
340 const size_t stride = vec.getStride();
341 const size_t flat_stride = stride * mp_size;
342 const size_t num_vecs = vec.getNumVectors();
343 RCP<FlatVector> flat_vec =
344 rcp(
new FlatVector(flat_map, flat_vec_vals, flat_stride, num_vecs));
359 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
360 if (flat_map == Teuchos::null) {
364 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
377 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
379 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv =
create_flat_vector_view(mv, flat_map);
380 return flat_mv->getVector(0);
393 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
394 if (flat_map == Teuchos::null) {
398 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
411 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
413 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv =
create_flat_vector_view(mv, flat_map);
414 return flat_mv->getVectorNonConst(0);
427 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
428 if (flat_map == Teuchos::null) {
432 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
436 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR) 444 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
448 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
450 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
455 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
456 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
457 typedef typename FlatVector::dual_view_type flat_view_type;
460 flat_view_type flat_vals = vec.getDualView();
463 RCP<FlatVector> flat_vec = rcp(
new FlatVector(flat_map, flat_vals));
474 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
478 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
480 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
485 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
486 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
487 typedef typename FlatVector::dual_view_type flat_view_type;
490 flat_view_type flat_vals = vec.getDualView();
493 RCP<FlatVector> flat_vec = rcp(
new FlatVector(flat_map, flat_vals));
509 const Teuchos::RCP<
const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& flat_graph,
511 using Teuchos::ArrayView;
512 using Teuchos::Array;
518 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
521 RCP<FlatMatrix> flat_mat = rcp(
new FlatMatrix(flat_graph));
524 const size_t num_rows = mat.getNodeNumRows();
525 const size_t max_cols = mat.getNodeMaxNumRowEntries();
526 ArrayView<const LocalOrdinal> indices, flat_indices;
527 ArrayView<const Scalar> values;
528 Array<BaseScalar> flat_values(max_cols);
529 for (
size_t row=0; row<num_rows; ++row) {
530 mat.getLocalRowView(row, indices, values);
531 const size_t num_col = mat.getNumEntriesInLocalRow(row);
534 for (
size_t j=0;
j<num_col; ++
j)
535 flat_values[
j] = values[
j].coeff(i);
536 flat_graph->getLocalRowView(flat_row, flat_indices);
537 flat_mat->replaceLocalValues(flat_row, flat_indices,
538 flat_values(0, num_col));
541 flat_mat->fillComplete(flat_graph->getDomainMap(),
542 flat_graph->getRangeMap());
549 #endif // STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP Stokhos::StandardStorage< int, double > Storage
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)
Top-level namespace for Stokhos classes and functions.
KokkosClassic::DefaultNode::DefaultNodeType Node
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_flat_mp_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_range_map, const LocalOrdinal block_size)
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)
Teuchos::RCP< Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > create_flat_map(const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map, const LocalOrdinal block_size)