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"
69 #include "Amesos2_Kokkos_View_Copy_Assign.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 
255 
269  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
270  struct same_gs_helper_kokkos_view
271  {
272  static void do_get(const Teuchos::Ptr<const M> mat,
273  KV_S& nzvals,
274  KV_GO& indices,
275  KV_GS& pointers,
276  typename KV_GS::value_type& nnz,
277  const Teuchos::Ptr<
278  const Tpetra::Map<typename M::local_ordinal_t,
279  typename M::global_ordinal_t,
280  typename M::node_t> > map,
281  EDistribution distribution,
282  EStorage_Ordering ordering)
283  {
284  Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
285  indices, pointers, nnz, map, distribution, ordering);
286  }
287  };
288 
289  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
290  struct diff_gs_helper_kokkos_view
291  {
292  static void do_get(const Teuchos::Ptr<const M> mat,
293  KV_S& nzvals,
294  KV_GO& indices,
295  KV_GS& pointers,
296  typename KV_GS::value_type& nnz,
297  const Teuchos::Ptr<
298  const Tpetra::Map<typename M::local_ordinal_t,
299  typename M::global_ordinal_t,
300  typename M::node_t> > map,
301  EDistribution distribution,
302  EStorage_Ordering ordering)
303  {
304  typedef typename M::global_size_t mat_gs_t;
305  typedef typename Kokkos::View<mat_gs_t*, Kokkos::HostSpace> KV_TMP;
306  size_t i, size = pointers.extent(0);
307  KV_TMP pointers_tmp(Kokkos::ViewAllocateWithoutInitializing("pointers_tmp"), size);
308 
309  mat_gs_t nnz_tmp = 0;
310  Op::template apply_kokkos_view<KV_S, KV_GO, KV_TMP>(mat, nzvals,
311  indices, pointers_tmp, nnz_tmp, Teuchos::ptrInArg(*map), distribution, ordering);
312  nnz = Teuchos::as<typename KV_GS::value_type>(nnz_tmp);
313 
314  typedef typename KV_GS::value_type view_gs_t;
315  for (i = 0; i < size; ++i){
316  pointers(i) = Teuchos::as<view_gs_t>(pointers_tmp(i));
317  }
318  nnz = Teuchos::as<view_gs_t>(nnz_tmp);
319  }
320  };
321 
322  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
323  struct same_go_helper_kokkos_view
324  {
325  static void do_get(const Teuchos::Ptr<const M> mat,
326  KV_S& nzvals,
327  KV_GO& indices,
328  KV_GS& pointers,
329  typename KV_GS::value_type& nnz,
330  const Teuchos::Ptr<
331  const Tpetra::Map<typename M::local_ordinal_t,
332  typename M::global_ordinal_t,
333  typename M::node_t> > map,
334  EDistribution distribution,
335  EStorage_Ordering ordering)
336  {
337  typedef typename M::global_size_t mat_gs_t;
338  typedef typename KV_GS::value_type view_gs_t;
339  if_then_else<is_same<view_gs_t,mat_gs_t>::value,
340  same_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
341  diff_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
342  pointers, nnz, map,
343  distribution, ordering);
344  }
345  };
346 
347  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
348  struct diff_go_helper_kokkos_view
349  {
350  static void do_get(const Teuchos::Ptr<const M> mat,
351  KV_S& nzvals,
352  KV_GO& indices,
353  KV_GS& pointers,
354  typename KV_GS::value_type& nnz,
355  const Teuchos::Ptr<
356  const Tpetra::Map<typename M::local_ordinal_t,
357  typename M::global_ordinal_t,
358  typename M::node_t> > map,
359  EDistribution distribution,
360  EStorage_Ordering ordering)
361  {
362  typedef typename M::global_ordinal_t mat_go_t;
363  typedef typename M::global_size_t mat_gs_t;
364  typedef typename Kokkos::View<mat_go_t*, Kokkos::HostSpace> KV_TMP;
365  size_t i, size = indices.extent(0);
366  KV_TMP indices_tmp(Kokkos::ViewAllocateWithoutInitializing("indices_tmp"), size);
367 
368  typedef typename KV_GO::value_type view_go_t;
369  typedef typename KV_GS::value_type view_gs_t;
370  if_then_else<is_same<view_gs_t,mat_gs_t>::value,
371  same_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op>,
372  diff_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op> >::type::do_get(mat, nzvals, indices_tmp,
373  pointers, nnz, map,
374  distribution, ordering);
375  for (i = 0; i < size; ++i){
376  indices(i) = Teuchos::as<view_go_t>(indices_tmp(i));
377  }
378  }
379  };
380 
381  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
382  struct same_scalar_helper_kokkos_view
383  {
384  static void do_get(const Teuchos::Ptr<const M> mat,
385  KV_S& nzvals,
386  KV_GO& indices,
387  KV_GS& pointers,
388  typename KV_GS::value_type& nnz,
389  const Teuchos::Ptr<
390  const Tpetra::Map<typename M::local_ordinal_t,
391  typename M::global_ordinal_t,
392  typename M::node_t> > map,
393  EDistribution distribution,
394  EStorage_Ordering ordering)
395  {
396  typedef typename M::global_ordinal_t mat_go_t;
397  typedef typename KV_GO::value_type view_go_t;
398  if_then_else<is_same<view_go_t,mat_go_t>::value,
399  same_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
400  diff_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
401  pointers, nnz, map,
402  distribution, ordering);
403  }
404  };
405 
406  template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
407  struct diff_scalar_helper_kokkos_view
408  {
409  static void do_get(const Teuchos::Ptr<const M> mat,
410  KV_S& nzvals,
411  KV_GO& indices,
412  KV_GS& pointers,
413  typename KV_GS::value_type& 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::global_ordinal_t mat_go_t;
422  typedef typename Kokkos::ArithTraits<typename M::scalar_t>::val_type mat_scalar_t;
423  typedef typename Kokkos::View<mat_scalar_t*, Kokkos::HostSpace> KV_TMP;
424  size_t i, size = nzvals.extent(0);
425  KV_TMP nzvals_tmp(Kokkos::ViewAllocateWithoutInitializing("nzvals_tmp"), size);
426 
427  typedef typename KV_S::value_type view_scalar_t;
428  typedef typename KV_GO::value_type view_go_t;
429  if_then_else<is_same<view_go_t,mat_go_t>::value,
430  same_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op>,
431  diff_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals_tmp, indices,
432  pointers, nnz, map,
433  distribution, ordering);
434 
435  for (i = 0; i < size; ++i){
436  nzvals(i) = Teuchos::as<view_scalar_t>(nzvals_tmp(i));
437  }
438  }
439  };
440 
441 
442  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS, class Op>
443  struct get_cxs_helper_kokkos_view
444  {
445  static void do_get(const Teuchos::Ptr<const Matrix> mat,
446  KV_S& nzvals,
447  KV_GO& indices,
448  KV_GS& pointers,
449  typename KV_GS::value_type& nnz,
450  EDistribution distribution,
451  EStorage_Ordering ordering=ARBITRARY,
452  typename KV_GO::value_type indexBase = 0)
453  {
454  typedef typename Matrix::local_ordinal_t lo_t;
455  typedef typename Matrix::global_ordinal_t go_t;
456  typedef typename Matrix::global_size_t gs_t;
457  typedef typename Matrix::node_t node_t;
458 
459  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
460  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
461  Op::get_dimension(mat),
462  mat->getComm(),
463  indexBase,
464  Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
465  );
466  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
467  }
468 
473  static void do_get(const Teuchos::Ptr<const Matrix> mat,
474  KV_S& nzvals,
475  KV_GO& indices,
476  KV_GS& pointers,
477  typename KV_GS::value_type& nnz,
478  EDistribution distribution, // Does this one need a distribution argument??
479  EStorage_Ordering ordering=ARBITRARY)
480  {
481  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
482  typename Matrix::global_ordinal_t,
483  typename Matrix::node_t> > map
484  = Op::getMap(mat);
485  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
486  }
487 
492  static void do_get(const Teuchos::Ptr<const Matrix> mat,
493  KV_S& nzvals,
494  KV_GO& indices,
495  KV_GS& pointers,
496  typename KV_GS::value_type& nnz,
497  const Teuchos::Ptr<
498  const Tpetra::Map<typename Matrix::local_ordinal_t,
499  typename Matrix::global_ordinal_t,
500  typename Matrix::node_t> > map,
501  EDistribution distribution,
502  EStorage_Ordering ordering=ARBITRARY)
503  {
504  typedef typename Matrix::scalar_t mat_scalar;
505  typedef typename KV_S::value_type view_scalar_t;
506 
507  if_then_else<is_same<mat_scalar,view_scalar_t>::value,
508  same_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op>,
509  diff_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op> >::type::do_get(mat,
510  nzvals, indices,
511  pointers, nnz,
512  map,
513  distribution, ordering);
514  }
515  };
516 
517 #ifndef DOXYGEN_SHOULD_SKIP_THIS
518  /*
519  * These two function-like classes are meant to be used as the \c
520  * Op template parameter for the \c get_cxs_helper template class.
521  */
522  template<class Matrix>
523  struct get_ccs_func
524  {
525  template<typename KV_S, typename KV_GO, typename KV_GS>
526  static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
527  KV_S& nzvals,
528  KV_GO& rowind,
529  KV_GS& colptr,
530  typename Matrix::global_size_t& nnz,
531  const Teuchos::Ptr<
532  const Tpetra::Map<typename Matrix::local_ordinal_t,
533  typename Matrix::global_ordinal_t,
534  typename Matrix::node_t> > map,
535  EDistribution distribution,
536  EStorage_Ordering ordering)
537  {
538  mat->getCcs_kokkos_view(nzvals, rowind, colptr, nnz, map, ordering, distribution);
539  }
540 
541  static
542  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
543  typename Matrix::global_ordinal_t,
544  typename Matrix::node_t> >
545  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
546  {
547  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
548  }
549 
550  static
551  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
552  typename Matrix::global_ordinal_t,
553  typename Matrix::node_t> >
554  getMap(const Teuchos::Ptr<const Matrix> mat)
555  {
556  return mat->getColMap();
557  }
558 
559  static
560  typename Matrix::global_size_t
561  get_dimension(const Teuchos::Ptr<const Matrix> mat)
562  {
563  return mat->getGlobalNumCols();
564  }
565  };
566 
567  template<class Matrix>
568  struct get_crs_func
569  {
570  template<typename KV_S, typename KV_GO, typename KV_GS>
571  static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
572  KV_S& nzvals,
573  KV_GO& colind,
574  KV_GS& rowptr,
575  typename Matrix::global_size_t& nnz,
576  const Teuchos::Ptr<
577  const Tpetra::Map<typename Matrix::local_ordinal_t,
578  typename Matrix::global_ordinal_t,
579  typename Matrix::node_t> > map,
580  EDistribution distribution,
581  EStorage_Ordering ordering)
582  {
583  mat->getCrs_kokkos_view(nzvals, colind, rowptr, nnz, map, ordering, distribution);
584  }
585 
586  static
587  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
588  typename Matrix::global_ordinal_t,
589  typename Matrix::node_t> >
590  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
591  {
592  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
593  }
594 
595  static
596  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
597  typename Matrix::global_ordinal_t,
598  typename Matrix::node_t> >
599  getMap(const Teuchos::Ptr<const Matrix> mat)
600  {
601  return mat->getRowMap();
602  }
603 
604  static
605  typename Matrix::global_size_t
606  get_dimension(const Teuchos::Ptr<const Matrix> mat)
607  {
608  return mat->getGlobalNumRows();
609  }
610  };
611 #endif // DOXYGEN_SHOULD_SKIP_THIS
612 
650  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
651  struct get_ccs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_ccs_func<Matrix> >
652  {};
653 
661  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
662  struct get_crs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_crs_func<Matrix> >
663  {};
664  /* End Matrix/MultiVector Utilities */
665 
666 
668  // Definitions //
670 
671 
672  template <typename LO, typename GO, typename GS, typename Node>
673  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
674  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
675  {
676  //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
677  Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
678  return gather_map;
679  }
680 
681 
682  template <typename LO, typename GO, typename GS, typename Node>
683  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
684  getDistributionMap(EDistribution distribution,
685  GS num_global_elements,
686  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
687  GO indexBase,
688  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
689  {
690  // TODO: Need to add indexBase to cases other than ROOTED
691  // We do not support these maps in any solver now.
692  switch( distribution ){
693  case DISTRIBUTED:
694  case DISTRIBUTED_NO_OVERLAP:
695  return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
696  case GLOBALLY_REPLICATED:
697  return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
698  case ROOTED:
699  {
700  int rank = Teuchos::rank(*comm);
701  size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
702  if( rank == 0 ) { my_num_elems = num_global_elements; }
703 
704  return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
705  my_num_elems, indexBase, comm));
706  }
707  case CONTIGUOUS_AND_ROOTED:
708  {
709  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
710  = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
711  return gathermap;
712  }
713  default:
714  TEUCHOS_TEST_FOR_EXCEPTION( true,
715  std::logic_error,
716  "Control should never reach this point. "
717  "Please contact the Amesos2 developers." );
718  }
719  }
720 
721 
722 #ifdef HAVE_AMESOS2_EPETRA
723 
724  //#pragma message "include 3"
725  //#include <Epetra_Map.h>
726 
727  template <typename LO, typename GO, typename GS, typename Node>
728  Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
729  epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
730  {
731  using Teuchos::as;
732  using Teuchos::rcp;
733 
734  int num_my_elements = map.NumMyElements();
735  Teuchos::Array<int> my_global_elements(num_my_elements);
736  map.MyGlobalElements(my_global_elements.getRawPtr());
737 
738  Teuchos::Array<GO> my_gbl_inds_buf;
739  Teuchos::ArrayView<GO> my_gbl_inds;
740  if (! std::is_same<int, GO>::value) {
741  my_gbl_inds_buf.resize (num_my_elements);
742  my_gbl_inds = my_gbl_inds_buf ();
743  for (int k = 0; k < num_my_elements; ++k) {
744  my_gbl_inds[k] = static_cast<GO> (my_global_elements[k]);
745  }
746  }
747  else {
748  using Teuchos::av_reinterpret_cast;
749  my_gbl_inds = av_reinterpret_cast<GO> (my_global_elements ());
750  }
751 
752  typedef Tpetra::Map<LO,GO,Node> map_t;
753  RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
754  my_gbl_inds(),
755  as<GO>(map.IndexBase()),
756  to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
757  return tmap;
758  }
759 
760  template <typename LO, typename GO, typename GS, typename Node>
761  Teuchos::RCP<Epetra_Map>
762  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
763  {
764  using Teuchos::as;
765 
766  Teuchos::Array<GO> elements_tmp;
767  elements_tmp = map.getLocalElementList();
768  int num_my_elements = elements_tmp.size();
769  Teuchos::Array<int> my_global_elements(num_my_elements);
770  for (int i = 0; i < num_my_elements; ++i){
771  my_global_elements[i] = as<int>(elements_tmp[i]);
772  }
773 
774  using Teuchos::rcp;
775  RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
776  num_my_elements,
777  my_global_elements.getRawPtr(),
778  as<GO>(map.getIndexBase()),
779  *to_epetra_comm(map.getComm())));
780  return emap;
781  }
782 #endif // HAVE_AMESOS2_EPETRA
783 
784  template <typename Scalar,
785  typename GlobalOrdinal,
786  typename GlobalSizeT>
787  void transpose(Teuchos::ArrayView<Scalar> vals,
788  Teuchos::ArrayView<GlobalOrdinal> indices,
789  Teuchos::ArrayView<GlobalSizeT> ptr,
790  Teuchos::ArrayView<Scalar> trans_vals,
791  Teuchos::ArrayView<GlobalOrdinal> trans_indices,
792  Teuchos::ArrayView<GlobalSizeT> trans_ptr)
793  {
794  /* We have a compressed-row storage format of this matrix. We
795  * transform this into a compressed-column format using a
796  * distribution-counting sort algorithm, which is described by
797  * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
798  */
799 
800 #ifdef HAVE_AMESOS2_DEBUG
801  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
802  ind_begin = indices.begin();
803  ind_end = indices.end();
804  size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
805  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
806  std::invalid_argument,
807  "Transpose pointer size not large enough." );
808  TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
809  std::invalid_argument,
810  "Transpose values array not large enough." );
811  TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
812  std::invalid_argument,
813  "Transpose indices array not large enough." );
814 #else
815  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
816 #endif
817  // Count the number of entries in each column
818  Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
819  ind_end = indices.end();
820  for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
821  ++(count[(*ind_it) + 1]);
822  }
823  // Accumulate
824  typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
825  cnt_end = count.end();
826  for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
827  *cnt_it = *cnt_it + *(cnt_it - 1);
828  }
829  // This becomes the array of column pointers
830  trans_ptr.assign(count);
831 
832  /* Move the nonzero values into their final place in nzval, based on the
833  * counts found previously.
834  *
835  * This sequence deviates from Knuth's algorithm a bit, following more
836  * closely the description presented in Gustavson, Fred G. "Two Fast
837  * Algorithms for Sparse Matrices: Multiplication and Permuted
838  * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
839  * 250--269, http://doi.acm.org/10.1145/355791.355796.
840  *
841  * The output indices end up in sorted order
842  */
843 
844  GlobalSizeT size = ptr.size();
845  for( GlobalSizeT i = 0; i < size - 1; ++i ){
846  GlobalOrdinal u = ptr[i];
847  GlobalOrdinal v = ptr[i + 1];
848  for( GlobalOrdinal j = u; j < v; ++j ){
849  GlobalOrdinal k = count[indices[j]];
850  trans_vals[k] = vals[j];
851  trans_indices[k] = i;
852  ++(count[indices[j]]);
853  }
854  }
855  }
856 
857 
858  template <typename Scalar1, typename Scalar2>
859  void
860  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
861  size_t ld, Teuchos::ArrayView<Scalar2> s)
862  {
863  size_t vals_size = vals.size();
864 #ifdef HAVE_AMESOS2_DEBUG
865  size_t s_size = s.size();
866  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
867  std::invalid_argument,
868  "Scale vector must have length at least that of the vector" );
869 #endif
870  size_t i, s_i;
871  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
872  if( s_i == l ){
873  // bring i to the next multiple of ld
874  i += ld - s_i;
875  s_i = 0;
876  }
877  vals[i] *= s[s_i];
878  }
879  }
880 
881  template <typename Scalar1, typename Scalar2, class BinaryOp>
882  void
883  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
884  size_t ld, Teuchos::ArrayView<Scalar2> s,
885  BinaryOp binary_op)
886  {
887  size_t vals_size = vals.size();
888 #ifdef HAVE_AMESOS2_DEBUG
889  size_t s_size = s.size();
890  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
891  std::invalid_argument,
892  "Scale vector must have length at least that of the vector" );
893 #endif
894  size_t i, s_i;
895  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
896  if( s_i == l ){
897  // bring i to the next multiple of ld
898  i += ld - s_i;
899  s_i = 0;
900  }
901  vals[i] = binary_op(vals[i], s[s_i]);
902  }
903  }
904 
905  template<class row_ptr_view_t, class cols_view_t, class per_view_t>
906  void
907  reorder(row_ptr_view_t & row_ptr, cols_view_t & cols,
908  per_view_t & perm, per_view_t & peri, size_t & nnz,
909  bool permute_matrix)
910  {
911  #ifndef HAVE_AMESOS2_METIS
912  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
913  "Cannot reorder for cuSolver because no METIS is available.");
914  #else
915  typedef typename cols_view_t::value_type ordinal_type;
916  typedef typename row_ptr_view_t::value_type size_type;
917 
918  // begin on host where we'll run metis reorder
919  auto host_row_ptr = Kokkos::create_mirror_view(row_ptr);
920  auto host_cols = Kokkos::create_mirror_view(cols);
921  Kokkos::deep_copy(host_row_ptr, row_ptr);
922  Kokkos::deep_copy(host_cols, cols);
923 
924  // strip out the diagonals - metis will just crash with them included.
925  // make space for the stripped version
926  typedef Kokkos::View<idx_t*, Kokkos::HostSpace> host_metis_array;
927  const ordinal_type size = row_ptr.size() - 1;
928  size_type max_nnz = host_row_ptr(size);
929  host_metis_array host_strip_diag_row_ptr(
930  Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_row_ptr"),
931  size+1);
932  host_metis_array host_strip_diag_cols(
933  Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_cols"),
934  max_nnz);
935 
936  size_type new_nnz = 0;
937  for(ordinal_type i = 0; i < size; ++i) {
938  host_strip_diag_row_ptr(i) = new_nnz;
939  for(size_type j = host_row_ptr(i); j < host_row_ptr(i+1); ++j) {
940  if (i != host_cols(j)) {
941  host_strip_diag_cols(new_nnz++) = host_cols(j);
942  }
943  }
944  }
945  host_strip_diag_row_ptr(size) = new_nnz;
946 
947  // we'll get original permutations on host
948  host_metis_array host_perm(
949  Kokkos::ViewAllocateWithoutInitializing("host_perm"), size);
950  host_metis_array host_peri(
951  Kokkos::ViewAllocateWithoutInitializing("host_peri"), size);
952 
953  // If we want to remove metis.h included in this header we can move this
954  // to the cpp, but we need to decide how to handle the idx_t declaration.
955  idx_t metis_size = size;
956  int err = METIS_NodeND(&metis_size, host_strip_diag_row_ptr.data(), host_strip_diag_cols.data(),
957  NULL, NULL, host_perm.data(), host_peri.data());
958 
959  TEUCHOS_TEST_FOR_EXCEPTION(err != METIS_OK, std::runtime_error,
960  "METIS_NodeND failed to sort matrix.");
961 
962  // put the permutations on our saved device ptrs
963  // these will be used to permute x and b when we solve
964  typedef typename cols_view_t::execution_space exec_space_t;
965  auto device_perm = Kokkos::create_mirror_view(exec_space_t(), host_perm);
966  auto device_peri = Kokkos::create_mirror_view(exec_space_t(), host_peri);
967  deep_copy(device_perm, host_perm);
968  deep_copy(device_peri, host_peri);
969 
970  // also set the permutation which may need to convert the type from
971  // metis to the native ordinal_type
972  deep_copy_or_assign_view(perm, device_perm);
973  deep_copy_or_assign_view(peri, device_peri);
974 
975  if (permute_matrix) {
976  // we'll permute matrix on device to a new set of arrays
977  row_ptr_view_t new_row_ptr(
978  Kokkos::ViewAllocateWithoutInitializing("new_row_ptr"), row_ptr.size());
979  cols_view_t new_cols(
980  Kokkos::ViewAllocateWithoutInitializing("new_cols"), cols.size() - new_nnz/2);
981 
982  // permute row indices
983  Kokkos::RangePolicy<exec_space_t> policy_row(0, row_ptr.size());
984  Kokkos::parallel_scan(policy_row, KOKKOS_LAMBDA(
985  ordinal_type i, size_type & update, const bool &final) {
986  if(final) {
987  new_row_ptr(i) = update;
988  }
989  if(i < size) {
990  ordinal_type count = 0;
991  const ordinal_type row = device_perm(i);
992  for(ordinal_type k = row_ptr(row); k < row_ptr(row + 1); ++k) {
993  const ordinal_type j = device_peri(cols(k));
994  count += (i >= j);
995  }
996  update += count;
997  }
998  });
999 
1000  // permute col indices
1001  Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1002  Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1003  const ordinal_type kbeg = new_row_ptr(i);
1004  const ordinal_type row = device_perm(i);
1005  const ordinal_type col_beg = row_ptr(row);
1006  const ordinal_type col_end = row_ptr(row + 1);
1007  const ordinal_type nk = col_end - col_beg;
1008  for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1009  const ordinal_type tk = kbeg + t;
1010  const ordinal_type sk = col_beg + k;
1011  const ordinal_type j = device_peri(cols(sk));
1012  if(i >= j) {
1013  new_cols(tk) = j;
1014  ++t;
1015  }
1016  }
1017  });
1018 
1019  // finally set the inputs to the new sorted arrays
1020  row_ptr = new_row_ptr;
1021  cols = new_cols;
1022  }
1023 
1024  nnz = new_nnz;
1025  #endif // HAVE_AMESOS2_METIS
1026  }
1027 
1028  template<class values_view_t, class row_ptr_view_t,
1029  class cols_view_t, class per_view_t>
1030  void
1031  reorder_values(values_view_t & values, const row_ptr_view_t & orig_row_ptr,
1032  const row_ptr_view_t & new_row_ptr,
1033  const cols_view_t & orig_cols, const per_view_t & perm, const per_view_t & peri,
1034  size_t nnz)
1035  {
1036  typedef typename cols_view_t::value_type ordinal_type;
1037  typedef typename cols_view_t::execution_space exec_space_t;
1038 
1039  auto device_perm = Kokkos::create_mirror_view(exec_space_t(), perm);
1040  auto device_peri = Kokkos::create_mirror_view(exec_space_t(), peri);
1041  deep_copy(device_perm, perm);
1042  deep_copy(device_peri, peri);
1043 
1044  const ordinal_type size = orig_row_ptr.size() - 1;
1045 
1046  auto host_orig_row_ptr = Kokkos::create_mirror_view(orig_row_ptr);
1047  auto new_nnz = host_orig_row_ptr(size); // TODO: Maybe optimize this by caching
1048 
1049  values_view_t new_values(
1050  Kokkos::ViewAllocateWithoutInitializing("new_values"), values.size() - new_nnz/2);
1051 
1052  // permute col indices
1053  Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1054  Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1055  const ordinal_type kbeg = new_row_ptr(i);
1056  const ordinal_type row = device_perm(i);
1057  const ordinal_type col_beg = orig_row_ptr(row);
1058  const ordinal_type col_end = orig_row_ptr(row + 1);
1059  const ordinal_type nk = col_end - col_beg;
1060  for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1061  const ordinal_type tk = kbeg + t;
1062  const ordinal_type sk = col_beg + k;
1063  const ordinal_type j = device_peri(orig_cols(sk));
1064  if(i >= j) {
1065  new_values(tk) = values(sk);
1066  ++t;
1067  }
1068  }
1069  });
1070 
1071  values = new_values;
1072  }
1073 
1074  template<class array_view_t, class per_view_t>
1075  void
1076  apply_reorder_permutation(const array_view_t & array,
1077  array_view_t & permuted_array, const per_view_t & permutation) {
1078  if(permuted_array.extent(0) != array.extent(0) || permuted_array.extent(1) != array.extent(1)) {
1079  permuted_array = array_view_t(
1080  Kokkos::ViewAllocateWithoutInitializing("permuted_array"),
1081  array.extent(0), array.extent(1));
1082  }
1083  typedef typename array_view_t::execution_space exec_space_t;
1084  Kokkos::RangePolicy<exec_space_t> policy(0, array.extent(0));
1085  Kokkos::parallel_for(policy, KOKKOS_LAMBDA(size_t i) {
1086  for(size_t j = 0; j < array.extent(1); ++j) {
1087  permuted_array(i, j) = array(permutation(i), j);
1088  }
1089  });
1090  }
1091 
1094  } // end namespace Util
1095 
1096 } // end namespace Amesos2
1097 
1098 #endif // #ifndef AMESOS2_UTIL_HPP
Amesos2::Util::get_ccs_helper_kokkos_view
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:651
Amesos2::Util::scale
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
Amesos2::Util::same_gs_helper_kokkos_view
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:270
Amesos2_Meta.hpp
Provides some simple meta-programming utilities for Amesos2.
Amesos2::EStorage_Ordering
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
Amesos2::DISTRIBUTED_NO_OVERLAP
Definition: Amesos2_TypeDecl.hpp:125
Amesos2::Util::transpose
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
Amesos2_Kokkos_View_Copy_Assign.hpp
Copy or assign views based on memory spaces.
Amesos2::Util::to_teuchos_comm
const RCP< const Teuchos::Comm< int > > to_teuchos_comm(RCP< const Epetra_Comm > c)
Transform an Epetra_Comm object into a Teuchos::Comm object.
Amesos2
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
Amesos2::Util::printLine
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition: Amesos2_Util.cpp:119
Amesos2::GLOBALLY_REPLICATED
Definition: Amesos2_TypeDecl.hpp:126
Amesos2::Util::to_epetra_comm
const RCP< const Epetra_Comm > to_epetra_comm(RCP< const Teuchos::Comm< int > > c)
Transfrom a Teuchos::Comm object into an Epetra_Comm object.
Amesos2::Util::getGatherMap
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:674
Amesos2::Util::get_crs_helper_kokkos_view
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:662
Amesos2::Util::tpetra_map_to_epetra_map
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:762
Amesos2::DISTRIBUTED
Definition: Amesos2_TypeDecl.hpp:124
Amesos2_TypeDecl.hpp
Enum and other types declarations for Amesos2.
Amesos2::Util::epetra_map_to_tpetra_map
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:729
Amesos2::ROOTED
Definition: Amesos2_TypeDecl.hpp:127
Amesos2::EDistribution
EDistribution
Definition: Amesos2_TypeDecl.hpp:123
Amesos2::CONTIGUOUS_AND_ROOTED
Definition: Amesos2_TypeDecl.hpp:128
Generated by   doxygen 1.8.14