Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Util.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Amesos2: Templated Direct Sparse Solver Package
6 // Copyright 2011 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //
42 // @HEADER
43 
52 #ifndef AMESOS2_UTIL_HPP
53 #define AMESOS2_UTIL_HPP
54 
55 #include "Amesos2_config.h"
56 
57 #include "Teuchos_RCP.hpp"
58 #include "Teuchos_BLAS_types.hpp"
59 #include "Teuchos_Array.hpp"
60 #include "Teuchos_ArrayView.hpp"
61 #include "Teuchos_FancyOStream.hpp"
62 
63 #include <Tpetra_Map.hpp>
64 #include <Tpetra_DistObject_decl.hpp>
65 #include <Tpetra_ComputeGatherMap.hpp> // added for gather map... where is the best place??
66 
67 #include "Amesos2_TypeDecl.hpp"
68 #include "Amesos2_Meta.hpp"
70 
71 #ifdef HAVE_AMESOS2_EPETRA
72 #include <Epetra_Map.h>
73 #endif
74 
75 #ifdef HAVE_AMESOS2_METIS
76 #include "metis.h" // to discuss, remove from header?
77 #endif
78 
79 namespace Amesos2 {
80 
81  namespace Util {
82 
89  using Teuchos::RCP;
90  using Teuchos::ArrayView;
91 
92  using Meta::is_same;
93  using Meta::if_then_else;
94 
111  template <typename LO, typename GO, typename GS, typename Node>
112  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
113  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map );
114 
115 
116  template <typename LO, typename GO, typename GS, typename Node>
117  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
118  getDistributionMap(EDistribution distribution,
119  GS num_global_elements,
120  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
121  GO indexBase = 0,
122  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map = Teuchos::null);
123 
124 
125 #ifdef HAVE_AMESOS2_EPETRA
126 
132  template <typename LO, typename GO, typename GS, typename Node>
133  RCP<Tpetra::Map<LO,GO,Node> >
134  epetra_map_to_tpetra_map(const Epetra_BlockMap& map);
135 
141  template <typename LO, typename GO, typename GS, typename Node>
142  RCP<Epetra_Map>
143  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map);
144 
150  const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
151 
157  const RCP<const Epetra_Comm> to_epetra_comm(RCP<const Teuchos::Comm<int> > c);
158 #endif // HAVE_AMESOS2_EPETRA
159 
165  template <typename Scalar,
166  typename GlobalOrdinal,
167  typename GlobalSizeT>
168  void transpose(ArrayView<Scalar> vals,
169  ArrayView<GlobalOrdinal> indices,
170  ArrayView<GlobalSizeT> ptr,
171  ArrayView<Scalar> trans_vals,
172  ArrayView<GlobalOrdinal> trans_indices,
173  ArrayView<GlobalSizeT> trans_ptr);
174 
188  template <typename Scalar1, typename Scalar2>
189  void scale(ArrayView<Scalar1> vals, size_t l,
190  size_t ld, ArrayView<Scalar2> s);
191 
210  template <typename Scalar1, typename Scalar2, class BinaryOp>
211  void scale(ArrayView<Scalar1> vals, size_t l,
212  size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
213 
214 
216  void printLine( Teuchos::FancyOStream &out );
217 
218  // Helper function used to convert Kokkos::complex pointer
219  // to std::complex pointer; needed for optimized code path
220  // when retrieving the CRS raw pointers
221  template < class T0, class T1 >
222  struct getStdCplxType
223  {
224  using common_type = typename std::common_type<T0,T1>::type;
225  using type = common_type;
226  };
227 
228  template < class T0, class T1 >
229  struct getStdCplxType< T0, T1* >
230  {
231  using common_type = typename std::common_type<T0,T1>::type;
232  using type = common_type;
233  };
234 
235 #if defined(HAVE_TEUCHOS_COMPLEX) && defined(HAVE_AMESOS2_KOKKOS)
236  template < class T0 >
237  struct getStdCplxType< T0, Kokkos::complex<T0>* >
238  {
239  using type = std::complex<T0>;
240  };
241 
242  template < class T0 , class T1 >
243  struct getStdCplxType< T0, Kokkos::complex<T1>* >
244  {
245  using common_type = typename std::common_type<T0,T1>::type;
246  using type = std::complex<common_type>;
247  };
248 #endif
249 
251  // Matrix/MultiVector Utilities //
253 
254 #ifndef DOXYGEN_SHOULD_SKIP_THIS
255  /*
256  * The following represents a general way of getting a CRS or CCS
257  * representation of a matrix with implicit type conversions. The
258  * \c get_crs_helper and \c get_ccs_helper classes are templated
259  * on 4 types:
260  *
261  * - A Matrix type (conforming to the Amesos2 MatrixAdapter interface)
262  * - A scalar type
263  * - A global ordinal type
264  * - A global size type
265  *
266  * The last three template types correspond to the input argument
267  * types. For example, if the scalar type is \c double , then we
268  * require that the \c nzvals argument is a \c
269  * Teuchos::ArrayView<double> type.
270  *
271  * These helpers perform any type conversions that must be
272  * performed to go between the Matrix's types and the input types.
273  * If no conversions are necessary the static functions can be
274  * effectively inlined.
275  */
276 
277  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
278  template <class M, typename S, typename GO, typename GS, class Op>
279  struct same_gs_helper
280  {
281  static void do_get(const Teuchos::Ptr<const M> mat,
282  const ArrayView<typename M::scalar_t> nzvals,
283  const ArrayView<typename M::global_ordinal_t> indices,
284  const ArrayView<GS> pointers,
285  GS& nnz,
286  const Teuchos::Ptr<
287  const Tpetra::Map<typename M::local_ordinal_t,
288  typename M::global_ordinal_t,
289  typename M::node_t> > map,
290  EDistribution distribution,
291  EStorage_Ordering ordering)
292  {
293  Op::apply(mat, nzvals, indices, pointers, nnz, map, distribution, ordering);
294  }
295  };
296 
297  template <class M, typename S, typename GO, typename GS, class Op>
298  struct diff_gs_helper
299  {
300  static void do_get(const Teuchos::Ptr<const M> mat,
301  const ArrayView<typename M::scalar_t> nzvals,
302  const ArrayView<typename M::global_ordinal_t> indices,
303  const ArrayView<GS> pointers,
304  GS& nnz,
305  const Teuchos::Ptr<
306  const Tpetra::Map<typename M::local_ordinal_t,
307  typename M::global_ordinal_t,
308  typename M::node_t> > map,
309  EDistribution distribution,
310  EStorage_Ordering ordering)
311  {
312  typedef typename M::global_size_t mat_gs_t;
313  typename ArrayView<GS>::size_type i, size = pointers.size();
314  Teuchos::Array<mat_gs_t> pointers_tmp(size);
315  mat_gs_t nnz_tmp = 0;
316 
317  Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, distribution, ordering);
318 
319  for (i = 0; i < size; ++i){
320  pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
321  }
322  nnz = Teuchos::as<GS>(nnz_tmp);
323  }
324  };
325 
326  template <class M, typename S, typename GO, typename GS, class Op>
327  struct same_go_helper
328  {
329  static void do_get(const Teuchos::Ptr<const M> mat,
330  const ArrayView<typename M::scalar_t> nzvals,
331  const ArrayView<GO> indices,
332  const ArrayView<GS> pointers,
333  GS& nnz,
334  const Teuchos::Ptr<
335  const Tpetra::Map<typename M::local_ordinal_t,
336  typename M::global_ordinal_t,
337  typename M::node_t> > map,
338  EDistribution distribution,
339  EStorage_Ordering ordering)
340  {
341  typedef typename M::global_size_t mat_gs_t;
342  if_then_else<is_same<GS,mat_gs_t>::value,
343  same_gs_helper<M,S,GO,GS,Op>,
344  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
345  pointers, nnz, map,
346  distribution, ordering);
347  }
348  };
349 
350  template <class M, typename S, typename GO, typename GS, class Op>
351  struct diff_go_helper
352  {
353  static void do_get(const Teuchos::Ptr<const M> mat,
354  const ArrayView<typename M::scalar_t> nzvals,
355  const ArrayView<GO> indices,
356  const ArrayView<GS> pointers,
357  GS& nnz,
358  const Teuchos::Ptr<
359  const Tpetra::Map<typename M::local_ordinal_t,
360  typename M::global_ordinal_t,
361  typename M::node_t> > map,
362  EDistribution distribution,
363  EStorage_Ordering ordering)
364  {
365  typedef typename M::global_ordinal_t mat_go_t;
366  typedef typename M::global_size_t mat_gs_t;
367  typename ArrayView<GO>::size_type i, size = indices.size();
368  Teuchos::Array<mat_go_t> indices_tmp(size);
369 
370  if_then_else<is_same<GS,mat_gs_t>::value,
371  same_gs_helper<M,S,GO,GS,Op>,
372  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
373  pointers, nnz, map,
374  distribution, ordering);
375 
376  for (i = 0; i < size; ++i){
377  indices[i] = Teuchos::as<GO>(indices_tmp[i]);
378  }
379  }
380  };
381 
382  template <class M, typename S, typename GO, typename GS, class Op>
383  struct same_scalar_helper
384  {
385  static void do_get(const Teuchos::Ptr<const M> mat,
386  const ArrayView<S> nzvals,
387  const ArrayView<GO> indices,
388  const ArrayView<GS> pointers,
389  GS& nnz,
390  const Teuchos::Ptr<
391  const Tpetra::Map<typename M::local_ordinal_t,
392  typename M::global_ordinal_t,
393  typename M::node_t> > map,
394  EDistribution distribution,
395  EStorage_Ordering ordering)
396  {
397  typedef typename M::global_ordinal_t mat_go_t;
398  if_then_else<is_same<GO,mat_go_t>::value,
399  same_go_helper<M,S,GO,GS,Op>,
400  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
401  pointers, nnz, map,
402  distribution, ordering);
403  }
404  };
405 
406  template <class M, typename S, typename GO, typename GS, class Op>
407  struct diff_scalar_helper
408  {
409  static void do_get(const Teuchos::Ptr<const M> mat,
410  const ArrayView<S> nzvals,
411  const ArrayView<GO> indices,
412  const ArrayView<GS> pointers,
413  GS& nnz,
414  const Teuchos::Ptr<
415  const Tpetra::Map<typename M::local_ordinal_t,
416  typename M::global_ordinal_t,
417  typename M::node_t> > map,
418  EDistribution distribution,
419  EStorage_Ordering ordering)
420  {
421  typedef typename M::scalar_t mat_scalar_t;
422  typedef typename M::global_ordinal_t mat_go_t;
423  typename ArrayView<S>::size_type i, size = nzvals.size();
424  Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
425 
426  if_then_else<is_same<GO,mat_go_t>::value,
427  same_go_helper<M,S,GO,GS,Op>,
428  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
429  pointers, nnz, map,
430  distribution, ordering);
431 
432  for (i = 0; i < size; ++i){
433  nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
434  }
435  }
436  };
437  #endif
438 
439 #endif // DOXYGEN_SHOULD_SKIP_THIS
440 
453  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
454  template<class Matrix, typename S, typename GO, typename GS, class Op>
455  struct get_cxs_helper
456  {
457  static void do_get(const Teuchos::Ptr<const Matrix> mat,
458  const ArrayView<S> nzvals,
459  const ArrayView<GO> indices,
460  const ArrayView<GS> pointers,
461  GS& nnz,
462  EDistribution distribution,
463  EStorage_Ordering ordering=ARBITRARY,
464  GO indexBase = 0)
465  {
466  typedef typename Matrix::local_ordinal_t lo_t;
467  typedef typename Matrix::global_ordinal_t go_t;
468  typedef typename Matrix::global_size_t gs_t;
469  typedef typename Matrix::node_t node_t;
470 
471  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
472  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
473  Op::get_dimension(mat),
474  mat->getComm(),
475  indexBase,
476  Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
477  );
478 
479  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
480  }
481 
486  static void do_get(const Teuchos::Ptr<const Matrix> mat,
487  const ArrayView<S> nzvals,
488  const ArrayView<GO> indices,
489  const ArrayView<GS> pointers,
490  GS& nnz,
491  EDistribution distribution, // Does this one need a distribution argument??
492  EStorage_Ordering ordering=ARBITRARY)
493  {
494  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
495  typename Matrix::global_ordinal_t,
496  typename Matrix::node_t> > map
497  = Op::getMap(mat);
498  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
499  }
500 
505  static void do_get(const Teuchos::Ptr<const Matrix> mat,
506  const ArrayView<S> nzvals,
507  const ArrayView<GO> indices,
508  const ArrayView<GS> pointers,
509  GS& nnz,
510  const Teuchos::Ptr<
511  const Tpetra::Map<typename Matrix::local_ordinal_t,
512  typename Matrix::global_ordinal_t,
513  typename Matrix::node_t> > map,
514  EDistribution distribution,
515  EStorage_Ordering ordering=ARBITRARY)
516  {
517  typedef typename Matrix::scalar_t mat_scalar;
518  if_then_else<is_same<mat_scalar,S>::value,
519  same_scalar_helper<Matrix,S,GO,GS,Op>,
520  diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
521  nzvals, indices,
522  pointers, nnz,
523  map,
524  distribution, ordering);
525  }
526  };
527  #endif
528 
529 
530  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
532  {
533  static void do_get(const Teuchos::Ptr<const M> mat,
534  KV_S& nzvals,
535  KV_GO& indices,
536  KV_GS& pointers,
537  typename KV_GS::value_type& nnz,
538  const Teuchos::Ptr<
539  const Tpetra::Map<typename M::local_ordinal_t,
540  typename M::global_ordinal_t,
541  typename M::node_t> > map,
542  EDistribution distribution,
543  EStorage_Ordering ordering)
544  {
545  Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
546  indices, pointers, nnz, map, distribution, ordering);
547  }
548  };
549 
550  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
551  struct diff_gs_helper_kokkos_view
552  {
553  static void do_get(const Teuchos::Ptr<const M> mat,
554  KV_S& nzvals,
555  KV_GO& indices,
556  KV_GS& pointers,
557  typename KV_GS::value_type& nnz,
558  const Teuchos::Ptr<
559  const Tpetra::Map<typename M::local_ordinal_t,
560  typename M::global_ordinal_t,
561  typename M::node_t> > map,
562  EDistribution distribution,
563  EStorage_Ordering ordering)
564  {
565  typedef typename M::global_size_t mat_gs_t;
566  typedef typename Kokkos::View<mat_gs_t*, Kokkos::HostSpace> KV_TMP;
567  size_t i, size = pointers.extent(0);
568  KV_TMP pointers_tmp(Kokkos::ViewAllocateWithoutInitializing("pointers_tmp"), size);
569 
570  mat_gs_t nnz_tmp = 0;
571  Op::template apply_kokkos_view<KV_S, KV_GO, KV_TMP>(mat, nzvals,
572  indices, pointers_tmp, nnz_tmp, Teuchos::ptrInArg(*map), distribution, ordering);
573  nnz = Teuchos::as<typename KV_GS::value_type>(nnz_tmp);
574 
575  typedef typename KV_GS::value_type view_gs_t;
576  for (i = 0; i < size; ++i){
577  pointers(i) = Teuchos::as<view_gs_t>(pointers_tmp(i));
578  }
579  nnz = Teuchos::as<view_gs_t>(nnz_tmp);
580  }
581  };
582 
583  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
584  struct same_go_helper_kokkos_view
585  {
586  static void do_get(const Teuchos::Ptr<const M> mat,
587  KV_S& nzvals,
588  KV_GO& indices,
589  KV_GS& pointers,
590  typename KV_GS::value_type& nnz,
591  const Teuchos::Ptr<
592  const Tpetra::Map<typename M::local_ordinal_t,
593  typename M::global_ordinal_t,
594  typename M::node_t> > map,
595  EDistribution distribution,
596  EStorage_Ordering ordering)
597  {
598  typedef typename M::global_size_t mat_gs_t;
599  typedef typename KV_GS::value_type view_gs_t;
600  if_then_else<is_same<view_gs_t,mat_gs_t>::value,
601  same_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
602  diff_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
603  pointers, nnz, map,
604  distribution, ordering);
605  }
606  };
607 
608  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
609  struct diff_go_helper_kokkos_view
610  {
611  static void do_get(const Teuchos::Ptr<const M> mat,
612  KV_S& nzvals,
613  KV_GO& indices,
614  KV_GS& pointers,
615  typename KV_GS::value_type& nnz,
616  const Teuchos::Ptr<
617  const Tpetra::Map<typename M::local_ordinal_t,
618  typename M::global_ordinal_t,
619  typename M::node_t> > map,
620  EDistribution distribution,
621  EStorage_Ordering ordering)
622  {
623  typedef typename M::global_ordinal_t mat_go_t;
624  typedef typename M::global_size_t mat_gs_t;
625  typedef typename Kokkos::View<mat_go_t*, Kokkos::HostSpace> KV_TMP;
626  size_t i, size = indices.extent(0);
627  KV_TMP indices_tmp(Kokkos::ViewAllocateWithoutInitializing("indices_tmp"), size);
628 
629  typedef typename KV_GO::value_type view_go_t;
630  typedef typename KV_GS::value_type view_gs_t;
631  if_then_else<is_same<view_gs_t,mat_gs_t>::value,
632  same_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op>,
633  diff_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op> >::type::do_get(mat, nzvals, indices_tmp,
634  pointers, nnz, map,
635  distribution, ordering);
636  for (i = 0; i < size; ++i){
637  indices(i) = Teuchos::as<view_go_t>(indices_tmp(i));
638  }
639  }
640  };
641 
642  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
643  struct same_scalar_helper_kokkos_view
644  {
645  static void do_get(const Teuchos::Ptr<const M> mat,
646  KV_S& nzvals,
647  KV_GO& indices,
648  KV_GS& pointers,
649  typename KV_GS::value_type& nnz,
650  const Teuchos::Ptr<
651  const Tpetra::Map<typename M::local_ordinal_t,
652  typename M::global_ordinal_t,
653  typename M::node_t> > map,
654  EDistribution distribution,
655  EStorage_Ordering ordering)
656  {
657  typedef typename M::global_ordinal_t mat_go_t;
658  typedef typename KV_GO::value_type view_go_t;
659  if_then_else<is_same<view_go_t,mat_go_t>::value,
660  same_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
661  diff_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
662  pointers, nnz, map,
663  distribution, ordering);
664  }
665  };
666 
667  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
668  struct diff_scalar_helper_kokkos_view
669  {
670  static void do_get(const Teuchos::Ptr<const M> mat,
671  KV_S& nzvals,
672  KV_GO& indices,
673  KV_GS& pointers,
674  typename KV_GS::value_type& nnz,
675  const Teuchos::Ptr<
676  const Tpetra::Map<typename M::local_ordinal_t,
677  typename M::global_ordinal_t,
678  typename M::node_t> > map,
679  EDistribution distribution,
680  EStorage_Ordering ordering)
681  {
682  typedef typename M::global_ordinal_t mat_go_t;
683  typedef typename Kokkos::ArithTraits<typename M::scalar_t>::val_type mat_scalar_t;
684  typedef typename Kokkos::View<mat_scalar_t*, Kokkos::HostSpace> KV_TMP;
685  size_t i, size = nzvals.extent(0);
686  KV_TMP nzvals_tmp(Kokkos::ViewAllocateWithoutInitializing("nzvals_tmp"), size);
687 
688  typedef typename KV_S::value_type view_scalar_t;
689  typedef typename KV_GO::value_type view_go_t;
690  if_then_else<is_same<view_go_t,mat_go_t>::value,
691  same_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op>,
692  diff_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals_tmp, indices,
693  pointers, nnz, map,
694  distribution, ordering);
695 
696  for (i = 0; i < size; ++i){
697  nzvals(i) = Teuchos::as<view_scalar_t>(nzvals_tmp(i));
698  }
699  }
700  };
701 
702 
703  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS, class Op>
704  struct get_cxs_helper_kokkos_view
705  {
706  static void do_get(const Teuchos::Ptr<const Matrix> mat,
707  KV_S& nzvals,
708  KV_GO& indices,
709  KV_GS& pointers,
710  typename KV_GS::value_type& nnz,
711  EDistribution distribution,
712  EStorage_Ordering ordering=ARBITRARY,
713  typename KV_GO::value_type indexBase = 0)
714  {
715  typedef typename Matrix::local_ordinal_t lo_t;
716  typedef typename Matrix::global_ordinal_t go_t;
717  typedef typename Matrix::global_size_t gs_t;
718  typedef typename Matrix::node_t node_t;
719 
720  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
721  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
722  Op::get_dimension(mat),
723  mat->getComm(),
724  indexBase,
725  Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
726  );
727  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
728  }
729 
734  static void do_get(const Teuchos::Ptr<const Matrix> mat,
735  KV_S& nzvals,
736  KV_GO& indices,
737  KV_GS& pointers,
738  typename KV_GS::value_type& nnz,
739  EDistribution distribution, // Does this one need a distribution argument??
740  EStorage_Ordering ordering=ARBITRARY)
741  {
742  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
743  typename Matrix::global_ordinal_t,
744  typename Matrix::node_t> > map
745  = Op::getMap(mat);
746  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
747  }
748 
753  static void do_get(const Teuchos::Ptr<const Matrix> mat,
754  KV_S& nzvals,
755  KV_GO& indices,
756  KV_GS& pointers,
757  typename KV_GS::value_type& nnz,
758  const Teuchos::Ptr<
759  const Tpetra::Map<typename Matrix::local_ordinal_t,
760  typename Matrix::global_ordinal_t,
761  typename Matrix::node_t> > map,
762  EDistribution distribution,
763  EStorage_Ordering ordering=ARBITRARY)
764  {
765  typedef typename Matrix::scalar_t mat_scalar;
766  typedef typename KV_S::value_type view_scalar_t;
767 
768  if_then_else<is_same<mat_scalar,view_scalar_t>::value,
769  same_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op>,
770  diff_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op> >::type::do_get(mat,
771  nzvals, indices,
772  pointers, nnz,
773  map,
774  distribution, ordering);
775  }
776  };
777 
778 #ifndef DOXYGEN_SHOULD_SKIP_THIS
779  /*
780  * These two function-like classes are meant to be used as the \c
781  * Op template parameter for the \c get_cxs_helper template class.
782  */
783  template<class Matrix>
784  struct get_ccs_func
785  {
786  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
787  static void apply(const Teuchos::Ptr<const Matrix> mat,
788  const ArrayView<typename Matrix::scalar_t> nzvals,
789  const ArrayView<typename Matrix::global_ordinal_t> rowind,
790  const ArrayView<typename Matrix::global_size_t> colptr,
791  typename Matrix::global_size_t& nnz,
792  const Teuchos::Ptr<
793  const Tpetra::Map<typename Matrix::local_ordinal_t,
794  typename Matrix::global_ordinal_t,
795  typename Matrix::node_t> > map,
796  EDistribution distribution,
797  EStorage_Ordering ordering)
798  {
799  mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering, distribution);
800  //mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
801  }
802  #endif
803 
804  template<typename KV_S, typename KV_GO, typename KV_GS>
805  static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
806  KV_S& nzvals,
807  KV_GO& rowind,
808  KV_GS& colptr,
809  typename Matrix::global_size_t& nnz,
810  const Teuchos::Ptr<
811  const Tpetra::Map<typename Matrix::local_ordinal_t,
812  typename Matrix::global_ordinal_t,
813  typename Matrix::node_t> > map,
814  EDistribution distribution,
815  EStorage_Ordering ordering)
816  {
817  mat->getCcs_kokkos_view(nzvals, rowind, colptr, nnz, map, ordering, distribution);
818  }
819 
820  static
821  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
822  typename Matrix::global_ordinal_t,
823  typename Matrix::node_t> >
824  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
825  {
826  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
827  }
828 
829  static
830  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
831  typename Matrix::global_ordinal_t,
832  typename Matrix::node_t> >
833  getMap(const Teuchos::Ptr<const Matrix> mat)
834  {
835  return mat->getColMap();
836  }
837 
838  static
839  typename Matrix::global_size_t
840  get_dimension(const Teuchos::Ptr<const Matrix> mat)
841  {
842  return mat->getGlobalNumCols();
843  }
844  };
845 
846  template<class Matrix>
847  struct get_crs_func
848  {
849  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
850  static void apply(const Teuchos::Ptr<const Matrix> mat,
851  const ArrayView<typename Matrix::scalar_t> nzvals,
852  const ArrayView<typename Matrix::global_ordinal_t> colind,
853  const ArrayView<typename Matrix::global_size_t> rowptr,
854  typename Matrix::global_size_t& nnz,
855  const Teuchos::Ptr<
856  const Tpetra::Map<typename Matrix::local_ordinal_t,
857  typename Matrix::global_ordinal_t,
858  typename Matrix::node_t> > map,
859  EDistribution distribution,
860  EStorage_Ordering ordering)
861  {
862  mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering, distribution);
863  }
864  #endif
865 
866  template<typename KV_S, typename KV_GO, typename KV_GS>
867  static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
868  KV_S& nzvals,
869  KV_GO& colind,
870  KV_GS& rowptr,
871  typename Matrix::global_size_t& nnz,
872  const Teuchos::Ptr<
873  const Tpetra::Map<typename Matrix::local_ordinal_t,
874  typename Matrix::global_ordinal_t,
875  typename Matrix::node_t> > map,
876  EDistribution distribution,
877  EStorage_Ordering ordering)
878  {
879  mat->getCrs_kokkos_view(nzvals, colind, rowptr, nnz, map, ordering, distribution);
880  }
881 
882  static
883  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
884  typename Matrix::global_ordinal_t,
885  typename Matrix::node_t> >
886  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
887  {
888  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
889  }
890 
891  static
892  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
893  typename Matrix::global_ordinal_t,
894  typename Matrix::node_t> >
895  getMap(const Teuchos::Ptr<const Matrix> mat)
896  {
897  return mat->getRowMap();
898  }
899 
900  static
901  typename Matrix::global_size_t
902  get_dimension(const Teuchos::Ptr<const Matrix> mat)
903  {
904  return mat->getGlobalNumRows();
905  }
906  };
907 #endif // DOXYGEN_SHOULD_SKIP_THIS
908 
946  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
947  template<class Matrix, typename S, typename GO, typename GS>
948  struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
949  {};
950  #endif
951 
952  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
953  struct get_ccs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_ccs_func<Matrix> >
954  {};
955 
963  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
964  template<class Matrix, typename S, typename GO, typename GS>
965  struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
966  {};
967  #endif
968 
969  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
970  struct get_crs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_crs_func<Matrix> >
971  {};
972  /* End Matrix/MultiVector Utilities */
973 
974 
976  // Definitions //
978 
979 
980  template <typename LO, typename GO, typename GS, typename Node>
981  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
982  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
983  {
984  //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
985  Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
986  return gather_map;
987  }
988 
989 
990  template <typename LO, typename GO, typename GS, typename Node>
991  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
992  getDistributionMap(EDistribution distribution,
993  GS num_global_elements,
994  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
995  GO indexBase,
996  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
997  {
998  // TODO: Need to add indexBase to cases other than ROOTED
999  // We do not support these maps in any solver now.
1000  switch( distribution ){
1001  case DISTRIBUTED:
1003  return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
1004  case GLOBALLY_REPLICATED:
1005  return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
1006  case ROOTED:
1007  {
1008  int rank = Teuchos::rank(*comm);
1009  size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
1010  if( rank == 0 ) { my_num_elems = num_global_elements; }
1011 
1012  return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
1013  my_num_elems, indexBase, comm));
1014  }
1015  case CONTIGUOUS_AND_ROOTED:
1016  {
1017  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
1018  = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
1019  return gathermap;
1020  }
1021  default:
1022  TEUCHOS_TEST_FOR_EXCEPTION( true,
1023  std::logic_error,
1024  "Control should never reach this point. "
1025  "Please contact the Amesos2 developers." );
1026  }
1027  }
1028 
1029 
1030 #ifdef HAVE_AMESOS2_EPETRA
1031 
1032  //#pragma message "include 3"
1033  //#include <Epetra_Map.h>
1034 
1035  template <typename LO, typename GO, typename GS, typename Node>
1036  Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
1037  epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
1038  {
1039  using Teuchos::as;
1040  using Teuchos::rcp;
1041 
1042  int num_my_elements = map.NumMyElements();
1043  Teuchos::Array<int> my_global_elements(num_my_elements);
1044  map.MyGlobalElements(my_global_elements.getRawPtr());
1045 
1046  Teuchos::Array<GO> my_gbl_inds_buf;
1047  Teuchos::ArrayView<GO> my_gbl_inds;
1048  if (! std::is_same<int, GO>::value) {
1049  my_gbl_inds_buf.resize (num_my_elements);
1050  my_gbl_inds = my_gbl_inds_buf ();
1051  for (int k = 0; k < num_my_elements; ++k) {
1052  my_gbl_inds[k] = static_cast<GO> (my_global_elements[k]);
1053  }
1054  }
1055  else {
1056  using Teuchos::av_reinterpret_cast;
1057  my_gbl_inds = av_reinterpret_cast<GO> (my_global_elements ());
1058  }
1059 
1060  typedef Tpetra::Map<LO,GO,Node> map_t;
1061  RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
1062  my_gbl_inds(),
1063  as<GO>(map.IndexBase()),
1064  to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
1065  return tmap;
1066  }
1067 
1068  template <typename LO, typename GO, typename GS, typename Node>
1069  Teuchos::RCP<Epetra_Map>
1070  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
1071  {
1072  using Teuchos::as;
1073 
1074  Teuchos::Array<GO> elements_tmp;
1075  elements_tmp = map.getNodeElementList();
1076  int num_my_elements = elements_tmp.size();
1077  Teuchos::Array<int> my_global_elements(num_my_elements);
1078  for (int i = 0; i < num_my_elements; ++i){
1079  my_global_elements[i] = as<int>(elements_tmp[i]);
1080  }
1081 
1082  using Teuchos::rcp;
1083  RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
1084  num_my_elements,
1085  my_global_elements.getRawPtr(),
1086  as<GO>(map.getIndexBase()),
1087  *to_epetra_comm(map.getComm())));
1088  return emap;
1089  }
1090 #endif // HAVE_AMESOS2_EPETRA
1091 
1092  template <typename Scalar,
1093  typename GlobalOrdinal,
1094  typename GlobalSizeT>
1095  void transpose(Teuchos::ArrayView<Scalar> vals,
1096  Teuchos::ArrayView<GlobalOrdinal> indices,
1097  Teuchos::ArrayView<GlobalSizeT> ptr,
1098  Teuchos::ArrayView<Scalar> trans_vals,
1099  Teuchos::ArrayView<GlobalOrdinal> trans_indices,
1100  Teuchos::ArrayView<GlobalSizeT> trans_ptr)
1101  {
1102  /* We have a compressed-row storage format of this matrix. We
1103  * transform this into a compressed-column format using a
1104  * distribution-counting sort algorithm, which is described by
1105  * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
1106  */
1107 
1108 #ifdef HAVE_AMESOS2_DEBUG
1109  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
1110  ind_begin = indices.begin();
1111  ind_end = indices.end();
1112  size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
1113  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
1114  std::invalid_argument,
1115  "Transpose pointer size not large enough." );
1116  TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
1117  std::invalid_argument,
1118  "Transpose values array not large enough." );
1119  TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
1120  std::invalid_argument,
1121  "Transpose indices array not large enough." );
1122 #else
1123  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
1124 #endif
1125  // Count the number of entries in each column
1126  Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
1127  ind_end = indices.end();
1128  for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
1129  ++(count[(*ind_it) + 1]);
1130  }
1131  // Accumulate
1132  typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
1133  cnt_end = count.end();
1134  for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
1135  *cnt_it = *cnt_it + *(cnt_it - 1);
1136  }
1137  // This becomes the array of column pointers
1138  trans_ptr.assign(count);
1139 
1140  /* Move the nonzero values into their final place in nzval, based on the
1141  * counts found previously.
1142  *
1143  * This sequence deviates from Knuth's algorithm a bit, following more
1144  * closely the description presented in Gustavson, Fred G. "Two Fast
1145  * Algorithms for Sparse Matrices: Multiplication and Permuted
1146  * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
1147  * 250--269, http://doi.acm.org/10.1145/355791.355796.
1148  *
1149  * The output indices end up in sorted order
1150  */
1151 
1152  GlobalSizeT size = ptr.size();
1153  for( GlobalSizeT i = 0; i < size - 1; ++i ){
1154  GlobalOrdinal u = ptr[i];
1155  GlobalOrdinal v = ptr[i + 1];
1156  for( GlobalOrdinal j = u; j < v; ++j ){
1157  GlobalOrdinal k = count[indices[j]];
1158  trans_vals[k] = vals[j];
1159  trans_indices[k] = i;
1160  ++(count[indices[j]]);
1161  }
1162  }
1163  }
1164 
1165 
1166  template <typename Scalar1, typename Scalar2>
1167  void
1168  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
1169  size_t ld, Teuchos::ArrayView<Scalar2> s)
1170  {
1171  size_t vals_size = vals.size();
1172 #ifdef HAVE_AMESOS2_DEBUG
1173  size_t s_size = s.size();
1174  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
1175  std::invalid_argument,
1176  "Scale vector must have length at least that of the vector" );
1177 #endif
1178  size_t i, s_i;
1179  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
1180  if( s_i == l ){
1181  // bring i to the next multiple of ld
1182  i += ld - s_i;
1183  s_i = 0;
1184  }
1185  vals[i] *= s[s_i];
1186  }
1187  }
1188 
1189  template <typename Scalar1, typename Scalar2, class BinaryOp>
1190  void
1191  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
1192  size_t ld, Teuchos::ArrayView<Scalar2> s,
1193  BinaryOp binary_op)
1194  {
1195  size_t vals_size = vals.size();
1196 #ifdef HAVE_AMESOS2_DEBUG
1197  size_t s_size = s.size();
1198  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
1199  std::invalid_argument,
1200  "Scale vector must have length at least that of the vector" );
1201 #endif
1202  size_t i, s_i;
1203  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
1204  if( s_i == l ){
1205  // bring i to the next multiple of ld
1206  i += ld - s_i;
1207  s_i = 0;
1208  }
1209  vals[i] = binary_op(vals[i], s[s_i]);
1210  }
1211  }
1212 
1213  template<class row_ptr_view_t, class cols_view_t, class per_view_t>
1214  void
1215  reorder(row_ptr_view_t & row_ptr, cols_view_t & cols,
1216  per_view_t & perm, per_view_t & peri, size_t & nnz,
1217  bool permute_matrix)
1218  {
1219  #ifndef HAVE_AMESOS2_METIS
1220  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
1221  "Cannot reorder for cuSolver because no METIS is available.");
1222  #else
1223  typedef typename cols_view_t::value_type ordinal_type;
1224  typedef typename row_ptr_view_t::value_type size_type;
1225 
1226  // begin on host where we'll run metis reorder
1227  auto host_row_ptr = Kokkos::create_mirror_view(row_ptr);
1228  auto host_cols = Kokkos::create_mirror_view(cols);
1229  Kokkos::deep_copy(host_row_ptr, row_ptr);
1230  Kokkos::deep_copy(host_cols, cols);
1231 
1232  // strip out the diagonals - metis will just crash with them included.
1233  // make space for the stripped version
1234  typedef Kokkos::View<idx_t*, Kokkos::HostSpace> host_metis_array;
1235  const ordinal_type size = row_ptr.size() - 1;
1236  size_type max_nnz = host_row_ptr(size);
1237  host_metis_array host_strip_diag_row_ptr(
1238  Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_row_ptr"),
1239  size+1);
1240  host_metis_array host_strip_diag_cols(
1241  Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_cols"),
1242  max_nnz);
1243 
1244  size_type new_nnz = 0;
1245  for(ordinal_type i = 0; i < size; ++i) {
1246  host_strip_diag_row_ptr(i) = new_nnz;
1247  for(size_type j = host_row_ptr(i); j < host_row_ptr(i+1); ++j) {
1248  if (i != host_cols(j)) {
1249  host_strip_diag_cols(new_nnz++) = host_cols(j);
1250  }
1251  }
1252  }
1253  host_strip_diag_row_ptr(size) = new_nnz;
1254 
1255  // we'll get original permutations on host
1256  host_metis_array host_perm(
1257  Kokkos::ViewAllocateWithoutInitializing("host_perm"), size);
1258  host_metis_array host_peri(
1259  Kokkos::ViewAllocateWithoutInitializing("host_peri"), size);
1260 
1261  // If we want to remove metis.h included in this header we can move this
1262  // to the cpp, but we need to decide how to handle the idx_t declaration.
1263  idx_t metis_size = size;
1264  int err = METIS_NodeND(&metis_size, host_strip_diag_row_ptr.data(), host_strip_diag_cols.data(),
1265  NULL, NULL, host_perm.data(), host_peri.data());
1266 
1267  TEUCHOS_TEST_FOR_EXCEPTION(err != METIS_OK, std::runtime_error,
1268  "METIS_NodeND failed to sort matrix.");
1269 
1270  // put the permutations on our saved device ptrs
1271  // these will be used to permute x and b when we solve
1272  typedef typename cols_view_t::execution_space exec_space_t;
1273  auto device_perm = Kokkos::create_mirror_view(exec_space_t(), host_perm);
1274  auto device_peri = Kokkos::create_mirror_view(exec_space_t(), host_peri);
1275  deep_copy(device_perm, host_perm);
1276  deep_copy(device_peri, host_peri);
1277 
1278  // also set the permutation which may need to convert the type from
1279  // metis to the native ordinal_type
1280  deep_copy_or_assign_view(perm, device_perm);
1281  deep_copy_or_assign_view(peri, device_peri);
1282 
1283  if (permute_matrix) {
1284  // we'll permute matrix on device to a new set of arrays
1285  row_ptr_view_t new_row_ptr(
1286  Kokkos::ViewAllocateWithoutInitializing("new_row_ptr"), row_ptr.size());
1287  cols_view_t new_cols(
1288  Kokkos::ViewAllocateWithoutInitializing("new_cols"), cols.size() - new_nnz/2);
1289 
1290  // permute row indices
1291  Kokkos::RangePolicy<exec_space_t> policy_row(0, row_ptr.size());
1292  Kokkos::parallel_scan(policy_row, KOKKOS_LAMBDA(
1293  ordinal_type i, size_type & update, const bool &final) {
1294  if(final) {
1295  new_row_ptr(i) = update;
1296  }
1297  if(i < size) {
1298  ordinal_type count = 0;
1299  const ordinal_type row = device_perm(i);
1300  for(ordinal_type k = row_ptr(row); k < row_ptr(row + 1); ++k) {
1301  const ordinal_type j = device_peri(cols(k));
1302  count += (i >= j);
1303  }
1304  update += count;
1305  }
1306  });
1307 
1308  // permute col indices
1309  Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1310  Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1311  const ordinal_type kbeg = new_row_ptr(i);
1312  const ordinal_type row = device_perm(i);
1313  const ordinal_type col_beg = row_ptr(row);
1314  const ordinal_type col_end = row_ptr(row + 1);
1315  const ordinal_type nk = col_end - col_beg;
1316  for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1317  const ordinal_type tk = kbeg + t;
1318  const ordinal_type sk = col_beg + k;
1319  const ordinal_type j = device_peri(cols(sk));
1320  if(i >= j) {
1321  new_cols(tk) = j;
1322  ++t;
1323  }
1324  }
1325  });
1326 
1327  // finally set the inputs to the new sorted arrays
1328  row_ptr = new_row_ptr;
1329  cols = new_cols;
1330  }
1331 
1332  nnz = new_nnz;
1333  #endif // HAVE_AMESOS2_METIS
1334  }
1335 
1336  template<class values_view_t, class row_ptr_view_t,
1337  class cols_view_t, class per_view_t>
1338  void
1339  reorder_values(values_view_t & values, const row_ptr_view_t & orig_row_ptr,
1340  const row_ptr_view_t & new_row_ptr,
1341  const cols_view_t & orig_cols, const per_view_t & perm, const per_view_t & peri,
1342  size_t nnz)
1343  {
1344  typedef typename cols_view_t::value_type ordinal_type;
1345  typedef typename cols_view_t::execution_space exec_space_t;
1346 
1347  auto device_perm = Kokkos::create_mirror_view(exec_space_t(), perm);
1348  auto device_peri = Kokkos::create_mirror_view(exec_space_t(), peri);
1349  deep_copy(device_perm, perm);
1350  deep_copy(device_peri, peri);
1351 
1352  const ordinal_type size = orig_row_ptr.size() - 1;
1353 
1354  auto host_orig_row_ptr = Kokkos::create_mirror_view(orig_row_ptr);
1355  auto new_nnz = host_orig_row_ptr(size); // TODO: Maybe optimize this by caching
1356 
1357  values_view_t new_values(
1358  Kokkos::ViewAllocateWithoutInitializing("new_values"), values.size() - new_nnz/2);
1359 
1360  // permute col indices
1361  Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1362  Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1363  const ordinal_type kbeg = new_row_ptr(i);
1364  const ordinal_type row = device_perm(i);
1365  const ordinal_type col_beg = orig_row_ptr(row);
1366  const ordinal_type col_end = orig_row_ptr(row + 1);
1367  const ordinal_type nk = col_end - col_beg;
1368  for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1369  const ordinal_type tk = kbeg + t;
1370  const ordinal_type sk = col_beg + k;
1371  const ordinal_type j = device_peri(orig_cols(sk));
1372  if(i >= j) {
1373  new_values(tk) = values(sk);
1374  ++t;
1375  }
1376  }
1377  });
1378 
1379  values = new_values;
1380  }
1381 
1382  template<class array_view_t, class per_view_t>
1383  void
1384  apply_reorder_permutation(const array_view_t & array,
1385  array_view_t & permuted_array, const per_view_t & permutation) {
1386  if(permuted_array.extent(0) != array.extent(0) || permuted_array.extent(1) != array.extent(1)) {
1387  permuted_array = array_view_t(
1388  Kokkos::ViewAllocateWithoutInitializing("permuted_array"),
1389  array.extent(0), array.extent(1));
1390  }
1391  typedef typename array_view_t::execution_space exec_space_t;
1392  Kokkos::RangePolicy<exec_space_t> policy(0, array.extent(0));
1393  Kokkos::parallel_for(policy, KOKKOS_LAMBDA(size_t i) {
1394  for(size_t j = 0; j < array.extent(1); ++j) {
1395  permuted_array(i, j) = array(permutation(i), j);
1396  }
1397  });
1398  }
1399 
1402  } // end namespace Util
1403 
1404 } // end namespace Amesos2
1405 
1406 #endif // #ifndef AMESOS2_UTIL_HPP
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:953
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:531
Provides some simple meta-programming utilities for Amesos2.
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
Definition: Amesos2_TypeDecl.hpp:125
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
Copy or assign views based on memory spaces.
const RCP< const Teuchos::Comm< int > > to_teuchos_comm(RCP< const Epetra_Comm > c)
Transform an Epetra_Comm object into a Teuchos::Comm object.
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition: Amesos2_Util.cpp:119
Definition: Amesos2_TypeDecl.hpp:126
const RCP< const Epetra_Comm > to_epetra_comm(RCP< const Teuchos::Comm< int > > c)
Transfrom a Teuchos::Comm object into an Epetra_Comm object.
const Teuchos::RCP< const Tpetra::Map< LO, GO, Node > > getGatherMap(const Teuchos::RCP< const Tpetra::Map< LO, GO, Node > > &map)
Gets a Tpetra::Map described by the EDistribution.
Definition: Amesos2_Util.hpp:982
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:970
Teuchos::RCP< Epetra_Map > tpetra_map_to_epetra_map(const Tpetra::Map< LO, GO, Node > &map)
Transform a Tpetra::Map object into an Epetra_Map.
Definition: Amesos2_Util.hpp:1070
Definition: Amesos2_TypeDecl.hpp:124
Enum and other types declarations for Amesos2.
Teuchos::RCP< Tpetra::Map< LO, GO, Node > > epetra_map_to_tpetra_map(const Epetra_BlockMap &map)
Transform an Epetra_Map object into a Tpetra::Map.
Definition: Amesos2_Util.hpp:1037
Definition: Amesos2_TypeDecl.hpp:127
EDistribution
Definition: Amesos2_TypeDecl.hpp:123
Definition: Amesos2_TypeDecl.hpp:128