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  template <class M, typename S, typename GO, typename GS, class Op>
278  struct same_gs_helper
279  {
280  static void do_get(const Teuchos::Ptr<const M> mat,
281  const ArrayView<typename M::scalar_t> nzvals,
282  const ArrayView<typename M::global_ordinal_t> indices,
283  const ArrayView<GS> pointers,
284  GS& nnz,
285  const Teuchos::Ptr<
286  const Tpetra::Map<typename M::local_ordinal_t,
287  typename M::global_ordinal_t,
288  typename M::node_t> > map,
289  EDistribution distribution,
290  EStorage_Ordering ordering)
291  {
292  Op::apply(mat, nzvals, indices, pointers, nnz, map, distribution, ordering);
293  }
294  };
295 
296  template <class M, typename S, typename GO, typename GS, class Op>
297  struct diff_gs_helper
298  {
299  static void do_get(const Teuchos::Ptr<const M> mat,
300  const ArrayView<typename M::scalar_t> nzvals,
301  const ArrayView<typename M::global_ordinal_t> indices,
302  const ArrayView<GS> pointers,
303  GS& nnz,
304  const Teuchos::Ptr<
305  const Tpetra::Map<typename M::local_ordinal_t,
306  typename M::global_ordinal_t,
307  typename M::node_t> > map,
308  EDistribution distribution,
309  EStorage_Ordering ordering)
310  {
311  typedef typename M::global_size_t mat_gs_t;
312  typename ArrayView<GS>::size_type i, size = pointers.size();
313  Teuchos::Array<mat_gs_t> pointers_tmp(size);
314  mat_gs_t nnz_tmp = 0;
315 
316  Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, distribution, ordering);
317 
318  for (i = 0; i < size; ++i){
319  pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
320  }
321  nnz = Teuchos::as<GS>(nnz_tmp);
322  }
323  };
324 
325  template <class M, typename S, typename GO, typename GS, class Op>
326  struct same_go_helper
327  {
328  static void do_get(const Teuchos::Ptr<const M> mat,
329  const ArrayView<typename M::scalar_t> nzvals,
330  const ArrayView<GO> indices,
331  const ArrayView<GS> pointers,
332  GS& nnz,
333  const Teuchos::Ptr<
334  const Tpetra::Map<typename M::local_ordinal_t,
335  typename M::global_ordinal_t,
336  typename M::node_t> > map,
337  EDistribution distribution,
338  EStorage_Ordering ordering)
339  {
340  typedef typename M::global_size_t mat_gs_t;
341  if_then_else<is_same<GS,mat_gs_t>::value,
342  same_gs_helper<M,S,GO,GS,Op>,
343  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
344  pointers, nnz, map,
345  distribution, ordering);
346  }
347  };
348 
349  template <class M, typename S, typename GO, typename GS, class Op>
350  struct diff_go_helper
351  {
352  static void do_get(const Teuchos::Ptr<const M> mat,
353  const ArrayView<typename M::scalar_t> nzvals,
354  const ArrayView<GO> indices,
355  const ArrayView<GS> pointers,
356  GS& nnz,
357  const Teuchos::Ptr<
358  const Tpetra::Map<typename M::local_ordinal_t,
359  typename M::global_ordinal_t,
360  typename M::node_t> > map,
361  EDistribution distribution,
362  EStorage_Ordering ordering)
363  {
364  typedef typename M::global_ordinal_t mat_go_t;
365  typedef typename M::global_size_t mat_gs_t;
366  typename ArrayView<GO>::size_type i, size = indices.size();
367  Teuchos::Array<mat_go_t> indices_tmp(size);
368 
369  if_then_else<is_same<GS,mat_gs_t>::value,
370  same_gs_helper<M,S,GO,GS,Op>,
371  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
372  pointers, nnz, map,
373  distribution, ordering);
374 
375  for (i = 0; i < size; ++i){
376  indices[i] = Teuchos::as<GO>(indices_tmp[i]);
377  }
378  }
379  };
380 
381  template <class M, typename S, typename GO, typename GS, class Op>
382  struct same_scalar_helper
383  {
384  static void do_get(const Teuchos::Ptr<const M> mat,
385  const ArrayView<S> nzvals,
386  const ArrayView<GO> indices,
387  const ArrayView<GS> pointers,
388  GS& 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  if_then_else<is_same<GO,mat_go_t>::value,
398  same_go_helper<M,S,GO,GS,Op>,
399  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
400  pointers, nnz, map,
401  distribution, ordering);
402  }
403  };
404 
405  template <class M, typename S, typename GO, typename GS, class Op>
406  struct diff_scalar_helper
407  {
408  static void do_get(const Teuchos::Ptr<const M> mat,
409  const ArrayView<S> nzvals,
410  const ArrayView<GO> indices,
411  const ArrayView<GS> pointers,
412  GS& nnz,
413  const Teuchos::Ptr<
414  const Tpetra::Map<typename M::local_ordinal_t,
415  typename M::global_ordinal_t,
416  typename M::node_t> > map,
417  EDistribution distribution,
418  EStorage_Ordering ordering)
419  {
420  typedef typename M::scalar_t mat_scalar_t;
421  typedef typename M::global_ordinal_t mat_go_t;
422  typename ArrayView<S>::size_type i, size = nzvals.size();
423  Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
424 
425  if_then_else<is_same<GO,mat_go_t>::value,
426  same_go_helper<M,S,GO,GS,Op>,
427  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
428  pointers, nnz, map,
429  distribution, ordering);
430 
431  for (i = 0; i < size; ++i){
432  nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
433  }
434  }
435  };
436 
437 #endif // DOXYGEN_SHOULD_SKIP_THIS
438 
451  template<class Matrix, typename S, typename GO, typename GS, class Op>
453  {
454  static void do_get(const Teuchos::Ptr<const Matrix> mat,
455  const ArrayView<S> nzvals,
456  const ArrayView<GO> indices,
457  const ArrayView<GS> pointers,
458  GS& nnz,
459  EDistribution distribution,
460  EStorage_Ordering ordering=ARBITRARY,
461  GO indexBase = 0)
462  {
463  typedef typename Matrix::local_ordinal_t lo_t;
464  typedef typename Matrix::global_ordinal_t go_t;
465  typedef typename Matrix::global_size_t gs_t;
466  typedef typename Matrix::node_t node_t;
467 
468  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
469  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
470  Op::get_dimension(mat),
471  mat->getComm(),
472  indexBase,
473  Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
474  );
475 
476  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
477  }
478 
483  static void do_get(const Teuchos::Ptr<const Matrix> mat,
484  const ArrayView<S> nzvals,
485  const ArrayView<GO> indices,
486  const ArrayView<GS> pointers,
487  GS& nnz,
488  EDistribution distribution, // Does this one need a distribution argument??
489  EStorage_Ordering ordering=ARBITRARY)
490  {
491  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
492  typename Matrix::global_ordinal_t,
493  typename Matrix::node_t> > map
494  = Op::getMap(mat);
495  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
496  }
497 
502  static void do_get(const Teuchos::Ptr<const Matrix> mat,
503  const ArrayView<S> nzvals,
504  const ArrayView<GO> indices,
505  const ArrayView<GS> pointers,
506  GS& nnz,
507  const Teuchos::Ptr<
508  const Tpetra::Map<typename Matrix::local_ordinal_t,
509  typename Matrix::global_ordinal_t,
510  typename Matrix::node_t> > map,
511  EDistribution distribution,
512  EStorage_Ordering ordering=ARBITRARY)
513  {
514  typedef typename Matrix::scalar_t mat_scalar;
515  if_then_else<is_same<mat_scalar,S>::value,
516  same_scalar_helper<Matrix,S,GO,GS,Op>,
517  diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
518  nzvals, indices,
519  pointers, nnz,
520  map,
521  distribution, ordering);
522  }
523  };
524 
525  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS, class Op>
526  struct get_cxs_helper_kokkos_view
527  {
528  static void do_get(const Teuchos::Ptr<const Matrix> mat,
529  KV_S& nzvals,
530  KV_GO& indices,
531  KV_GS& pointers,
532  typename KV_GS::value_type& nnz,
533  EDistribution distribution,
534  EStorage_Ordering ordering=ARBITRARY,
535  typename KV_GO::value_type indexBase = 0)
536  {
537  typedef typename Matrix::local_ordinal_t lo_t;
538  typedef typename Matrix::global_ordinal_t go_t;
539  typedef typename Matrix::global_size_t gs_t;
540  typedef typename Matrix::node_t node_t;
541 
542  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
543  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
544  Op::get_dimension(mat),
545  mat->getComm(),
546  indexBase,
547  Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
548  );
549  typename Matrix::global_size_t nnz_temp = 0; // only setting because Cuda gives warning used before unset
550  Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
551  indices, pointers, nnz_temp, Teuchos::ptrInArg(*map), distribution, ordering);
552  nnz = Teuchos::as<typename KV_GS::value_type>(nnz_temp);
553  }
554 
559  static void do_get(const Teuchos::Ptr<const Matrix> mat,
560  KV_S& nzvals,
561  KV_GO& indices,
562  KV_GS& pointers,
563  typename KV_GS::value_type& nnz,
564  EDistribution distribution, // Does this one need a distribution argument??
565  EStorage_Ordering ordering=ARBITRARY)
566  {
567  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
568  typename Matrix::global_ordinal_t,
569  typename Matrix::node_t> > map
570  = Op::getMap(mat);
571  typename Matrix::global_size_t nnz_temp;
572  Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
573  indices, pointers, nnz_temp, map, distribution, ordering);
574  nnz = Teuchos::as<typename KV_GS::value_type>(nnz_temp);
575  }
576 
581  static void do_get(const Teuchos::Ptr<const Matrix> mat,
582  KV_S& nzvals,
583  KV_GO& indices,
584  KV_GS& pointers,
585  typename KV_GS::value_type& nnz,
586  const Teuchos::Ptr<
587  const Tpetra::Map<typename Matrix::local_ordinal_t,
588  typename Matrix::global_ordinal_t,
589  typename Matrix::node_t> > map,
590  EDistribution distribution,
591  EStorage_Ordering ordering=ARBITRARY)
592  {
593  typename Matrix::global_size_t nnz_temp;
594  Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
595  indices, pointers, nnz_temp, map, distribution, ordering);
596  nnz = Teuchos::as<typename KV_GS::value_type>(nnz_temp);
597  }
598  };
599 
600 #ifndef DOXYGEN_SHOULD_SKIP_THIS
601  /*
602  * These two function-like classes are meant to be used as the \c
603  * Op template parameter for the \c get_cxs_helper template class.
604  */
605  template<class Matrix>
606  struct get_ccs_func
607  {
608  static void apply(const Teuchos::Ptr<const Matrix> mat,
609  const ArrayView<typename Matrix::scalar_t> nzvals,
610  const ArrayView<typename Matrix::global_ordinal_t> rowind,
611  const ArrayView<typename Matrix::global_size_t> colptr,
612  typename Matrix::global_size_t& nnz,
613  const Teuchos::Ptr<
614  const Tpetra::Map<typename Matrix::local_ordinal_t,
615  typename Matrix::global_ordinal_t,
616  typename Matrix::node_t> > map,
617  EDistribution distribution,
618  EStorage_Ordering ordering)
619  {
620  mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering, distribution);
621  //mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
622  }
623 
624  template<typename KV_S, typename KV_GO, typename KV_GS>
625  static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
626  KV_S& nzvals,
627  KV_GO& rowind,
628  KV_GS& colptr,
629  typename Matrix::global_size_t& nnz,
630  const Teuchos::Ptr<
631  const Tpetra::Map<typename Matrix::local_ordinal_t,
632  typename Matrix::global_ordinal_t,
633  typename Matrix::node_t> > map,
634  EDistribution distribution,
635  EStorage_Ordering ordering)
636  {
637  mat->getCcs_kokkos_view(nzvals, rowind, colptr, nnz, map, ordering, distribution);
638  }
639 
640  static
641  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
642  typename Matrix::global_ordinal_t,
643  typename Matrix::node_t> >
644  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
645  {
646  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
647  }
648 
649  static
650  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
651  typename Matrix::global_ordinal_t,
652  typename Matrix::node_t> >
653  getMap(const Teuchos::Ptr<const Matrix> mat)
654  {
655  return mat->getColMap();
656  }
657 
658  static
659  typename Matrix::global_size_t
660  get_dimension(const Teuchos::Ptr<const Matrix> mat)
661  {
662  return mat->getGlobalNumCols();
663  }
664  };
665 
666  template<class Matrix>
667  struct get_crs_func
668  {
669  static void apply(const Teuchos::Ptr<const Matrix> mat,
670  const ArrayView<typename Matrix::scalar_t> nzvals,
671  const ArrayView<typename Matrix::global_ordinal_t> colind,
672  const ArrayView<typename Matrix::global_size_t> rowptr,
673  typename Matrix::global_size_t& nnz,
674  const Teuchos::Ptr<
675  const Tpetra::Map<typename Matrix::local_ordinal_t,
676  typename Matrix::global_ordinal_t,
677  typename Matrix::node_t> > map,
678  EDistribution distribution,
679  EStorage_Ordering ordering)
680  {
681  mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering, distribution);
682  }
683 
684  template<typename KV_S, typename KV_GO, typename KV_GS>
685  static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
686  KV_S& nzvals,
687  KV_GO& colind,
688  KV_GS& rowptr,
689  typename Matrix::global_size_t& nnz,
690  const Teuchos::Ptr<
691  const Tpetra::Map<typename Matrix::local_ordinal_t,
692  typename Matrix::global_ordinal_t,
693  typename Matrix::node_t> > map,
694  EDistribution distribution,
695  EStorage_Ordering ordering)
696  {
697  mat->getCrs_kokkos_view(nzvals, colind, rowptr, nnz, map, ordering, distribution);
698  }
699 
700  static
701  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
702  typename Matrix::global_ordinal_t,
703  typename Matrix::node_t> >
704  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
705  {
706  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
707  }
708 
709  static
710  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
711  typename Matrix::global_ordinal_t,
712  typename Matrix::node_t> >
713  getMap(const Teuchos::Ptr<const Matrix> mat)
714  {
715  return mat->getRowMap();
716  }
717 
718  static
719  typename Matrix::global_size_t
720  get_dimension(const Teuchos::Ptr<const Matrix> mat)
721  {
722  return mat->getGlobalNumRows();
723  }
724  };
725 #endif // DOXYGEN_SHOULD_SKIP_THIS
726 
764  template<class Matrix, typename S, typename GO, typename GS>
765  struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
766  {};
767 
768  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
769  struct get_ccs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_ccs_func<Matrix> >
770  {};
771 
779  template<class Matrix, typename S, typename GO, typename GS>
780  struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
781  {};
782 
783  template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
784  struct get_crs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_crs_func<Matrix> >
785  {};
786  /* End Matrix/MultiVector Utilities */
787 
788 
790  // Definitions //
792 
793 
794  template <typename LO, typename GO, typename GS, typename Node>
795  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
796  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
797  {
798  //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
799  Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
800  return gather_map;
801  }
802 
803 
804  template <typename LO, typename GO, typename GS, typename Node>
805  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
806  getDistributionMap(EDistribution distribution,
807  GS num_global_elements,
808  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
809  GO indexBase,
810  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
811  {
812  // TODO: Need to add indexBase to cases other than ROOTED
813  // We do not support these maps in any solver now.
814  switch( distribution ){
815  case DISTRIBUTED:
816  case DISTRIBUTED_NO_OVERLAP:
817  return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
818  case GLOBALLY_REPLICATED:
819  return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
820  case ROOTED:
821  {
822  int rank = Teuchos::rank(*comm);
823  size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
824  if( rank == 0 ) { my_num_elems = num_global_elements; }
825 
826  return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
827  my_num_elems, indexBase, comm));
828  }
829  case CONTIGUOUS_AND_ROOTED:
830  {
831  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
832  = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
833  return gathermap;
834  }
835  default:
836  TEUCHOS_TEST_FOR_EXCEPTION( true,
837  std::logic_error,
838  "Control should never reach this point. "
839  "Please contact the Amesos2 developers." );
840  }
841  }
842 
843 
844 #ifdef HAVE_AMESOS2_EPETRA
845 
846  //#pragma message "include 3"
847  //#include <Epetra_Map.h>
848 
849  template <typename LO, typename GO, typename GS, typename Node>
850  Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
851  epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
852  {
853  using Teuchos::as;
854  using Teuchos::rcp;
855 
856  int num_my_elements = map.NumMyElements();
857  Teuchos::Array<int> my_global_elements(num_my_elements);
858  map.MyGlobalElements(my_global_elements.getRawPtr());
859 
860  Teuchos::Array<GO> my_gbl_inds_buf;
861  Teuchos::ArrayView<GO> my_gbl_inds;
862  if (! std::is_same<int, GO>::value) {
863  my_gbl_inds_buf.resize (num_my_elements);
864  my_gbl_inds = my_gbl_inds_buf ();
865  for (int k = 0; k < num_my_elements; ++k) {
866  my_gbl_inds[k] = static_cast<GO> (my_global_elements[k]);
867  }
868  }
869  else {
870  using Teuchos::av_reinterpret_cast;
871  my_gbl_inds = av_reinterpret_cast<GO> (my_global_elements ());
872  }
873 
874  typedef Tpetra::Map<LO,GO,Node> map_t;
875  RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
876  my_gbl_inds(),
877  as<GO>(map.IndexBase()),
878  to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
879  return tmap;
880  }
881 
882  template <typename LO, typename GO, typename GS, typename Node>
883  Teuchos::RCP<Epetra_Map>
884  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
885  {
886  using Teuchos::as;
887 
888  Teuchos::Array<GO> elements_tmp;
889  elements_tmp = map.getNodeElementList();
890  int num_my_elements = elements_tmp.size();
891  Teuchos::Array<int> my_global_elements(num_my_elements);
892  for (int i = 0; i < num_my_elements; ++i){
893  my_global_elements[i] = as<int>(elements_tmp[i]);
894  }
895 
896  using Teuchos::rcp;
897  RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
898  num_my_elements,
899  my_global_elements.getRawPtr(),
900  as<GO>(map.getIndexBase()),
901  *to_epetra_comm(map.getComm())));
902  return emap;
903  }
904 #endif // HAVE_AMESOS2_EPETRA
905 
906  template <typename Scalar,
907  typename GlobalOrdinal,
908  typename GlobalSizeT>
909  void transpose(Teuchos::ArrayView<Scalar> vals,
910  Teuchos::ArrayView<GlobalOrdinal> indices,
911  Teuchos::ArrayView<GlobalSizeT> ptr,
912  Teuchos::ArrayView<Scalar> trans_vals,
913  Teuchos::ArrayView<GlobalOrdinal> trans_indices,
914  Teuchos::ArrayView<GlobalSizeT> trans_ptr)
915  {
916  /* We have a compressed-row storage format of this matrix. We
917  * transform this into a compressed-column format using a
918  * distribution-counting sort algorithm, which is described by
919  * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
920  */
921 
922 #ifdef HAVE_AMESOS2_DEBUG
923  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
924  ind_begin = indices.begin();
925  ind_end = indices.end();
926  size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
927  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
928  std::invalid_argument,
929  "Transpose pointer size not large enough." );
930  TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
931  std::invalid_argument,
932  "Transpose values array not large enough." );
933  TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
934  std::invalid_argument,
935  "Transpose indices array not large enough." );
936 #else
937  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
938 #endif
939  // Count the number of entries in each column
940  Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
941  ind_end = indices.end();
942  for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
943  ++(count[(*ind_it) + 1]);
944  }
945  // Accumulate
946  typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
947  cnt_end = count.end();
948  for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
949  *cnt_it = *cnt_it + *(cnt_it - 1);
950  }
951  // This becomes the array of column pointers
952  trans_ptr.assign(count);
953 
954  /* Move the nonzero values into their final place in nzval, based on the
955  * counts found previously.
956  *
957  * This sequence deviates from Knuth's algorithm a bit, following more
958  * closely the description presented in Gustavson, Fred G. "Two Fast
959  * Algorithms for Sparse Matrices: Multiplication and Permuted
960  * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
961  * 250--269, http://doi.acm.org/10.1145/355791.355796.
962  *
963  * The output indices end up in sorted order
964  */
965 
966  GlobalSizeT size = ptr.size();
967  for( GlobalSizeT i = 0; i < size - 1; ++i ){
968  GlobalOrdinal u = ptr[i];
969  GlobalOrdinal v = ptr[i + 1];
970  for( GlobalOrdinal j = u; j < v; ++j ){
971  GlobalOrdinal k = count[indices[j]];
972  trans_vals[k] = vals[j];
973  trans_indices[k] = i;
974  ++(count[indices[j]]);
975  }
976  }
977  }
978 
979 
980  template <typename Scalar1, typename Scalar2>
981  void
982  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
983  size_t ld, Teuchos::ArrayView<Scalar2> s)
984  {
985  size_t vals_size = vals.size();
986 #ifdef HAVE_AMESOS2_DEBUG
987  size_t s_size = s.size();
988  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
989  std::invalid_argument,
990  "Scale vector must have length at least that of the vector" );
991 #endif
992  size_t i, s_i;
993  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
994  if( s_i == l ){
995  // bring i to the next multiple of ld
996  i += ld - s_i;
997  s_i = 0;
998  }
999  vals[i] *= s[s_i];
1000  }
1001  }
1002 
1003  template <typename Scalar1, typename Scalar2, class BinaryOp>
1004  void
1005  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
1006  size_t ld, Teuchos::ArrayView<Scalar2> s,
1007  BinaryOp binary_op)
1008  {
1009  size_t vals_size = vals.size();
1010 #ifdef HAVE_AMESOS2_DEBUG
1011  size_t s_size = s.size();
1012  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
1013  std::invalid_argument,
1014  "Scale vector must have length at least that of the vector" );
1015 #endif
1016  size_t i, s_i;
1017  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
1018  if( s_i == l ){
1019  // bring i to the next multiple of ld
1020  i += ld - s_i;
1021  s_i = 0;
1022  }
1023  vals[i] = binary_op(vals[i], s[s_i]);
1024  }
1025  }
1026 
1027  template<class row_ptr_view_t, class cols_view_t, class per_view_t>
1028  void
1029  reorder(row_ptr_view_t & row_ptr, cols_view_t & cols,
1030  per_view_t & perm, per_view_t & peri, size_t & nnz)
1031  {
1032  #ifndef HAVE_AMESOS2_METIS
1033  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
1034  "Cannot reorder for cuSolver because no METIS is available.");
1035  #else
1036  typedef typename cols_view_t::value_type ordinal_type;
1037  typedef typename row_ptr_view_t::value_type size_type;
1038 
1039  // begin on host where we'll run metis reorder
1040  auto host_row_ptr = Kokkos::create_mirror_view(row_ptr);
1041  auto host_cols = Kokkos::create_mirror_view(cols);
1042  Kokkos::deep_copy(host_row_ptr, row_ptr);
1043  Kokkos::deep_copy(host_cols, cols);
1044 
1045  // strip out the diagonals - metis will just crash with them included.
1046  // make space for the stripped version
1047  typedef Kokkos::View<idx_t*, Kokkos::HostSpace> host_metis_array;
1048  const ordinal_type size = row_ptr.size() - 1;
1049  host_metis_array host_strip_diag_row_ptr(
1050  Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_row_ptr"),
1051  row_ptr.size());
1052  host_metis_array host_strip_diag_cols(
1053  Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_cols"),
1054  cols.size() - size); // dropping diagonals.
1055 
1056  size_type new_nnz = 0;
1057  for(ordinal_type i = 0; i < size; ++i) {
1058  host_strip_diag_row_ptr(i) = new_nnz;
1059  for(size_type j = host_row_ptr(i); j < host_row_ptr(i+1); ++j) {
1060  if (i != host_cols(j)) {
1061  host_strip_diag_cols(new_nnz++) = host_cols(j);
1062  }
1063  }
1064  }
1065  host_strip_diag_row_ptr(size) = new_nnz;
1066 
1067  // we'll get original permutations on host
1068  host_metis_array host_perm(
1069  Kokkos::ViewAllocateWithoutInitializing("host_perm"), size);
1070  host_metis_array host_peri(
1071  Kokkos::ViewAllocateWithoutInitializing("host_peri"), size);
1072 
1073  // If we want to remove metis.h included in this header we can move this
1074  // to the cpp, but we need to decide how to handle the idx_t declaration.
1075  idx_t metis_size = size;
1076  int err = METIS_NodeND(&metis_size, host_strip_diag_row_ptr.data(), host_strip_diag_cols.data(),
1077  NULL, NULL, host_perm.data(), host_peri.data());
1078 
1079  TEUCHOS_TEST_FOR_EXCEPTION(err != METIS_OK, std::runtime_error,
1080  "METIS_NodeND failed to sort matrix.");
1081 
1082  // put the permutations on our saved device ptrs
1083  // these will be used to permute x and b when we solve
1084  typedef typename cols_view_t::execution_space exec_space_t;
1085  auto device_perm = Kokkos::create_mirror_view(exec_space_t(), host_perm);
1086  auto device_peri = Kokkos::create_mirror_view(exec_space_t(), host_peri);
1087  deep_copy(device_perm, host_perm);
1088  deep_copy(device_peri, host_peri);
1089 
1090  // also set the permutation which may need to convert the type from
1091  // metis to the native ordinal_type
1092  deep_copy_or_assign_view(perm, device_perm);
1093  deep_copy_or_assign_view(peri, device_peri);
1094 
1095  // we'll permute matrix on device to a new set of arrays
1096  row_ptr_view_t new_row_ptr(
1097  Kokkos::ViewAllocateWithoutInitializing("new_row_ptr"), row_ptr.size());
1098  cols_view_t new_cols(
1099  Kokkos::ViewAllocateWithoutInitializing("new_cols"), cols.size() - new_nnz/2);
1100 
1101  // permute row indices
1102  Kokkos::RangePolicy<exec_space_t> policy_row(0, row_ptr.size());
1103  Kokkos::parallel_scan(policy_row, KOKKOS_LAMBDA(
1104  ordinal_type i, size_type & update, const bool &final) {
1105  if(final) {
1106  new_row_ptr(i) = update;
1107  }
1108  if(i < size) {
1109  ordinal_type count = 0;
1110  const ordinal_type row = device_perm(i);
1111  for(ordinal_type k = row_ptr(row); k < row_ptr(row + 1); ++k) {
1112  const ordinal_type j = device_peri(cols(k));
1113  count += (i >= j);
1114  }
1115  update += count;
1116  }
1117  });
1118 
1119  // permute col indices
1120  Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1121  Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1122  const ordinal_type kbeg = new_row_ptr(i);
1123  const ordinal_type row = device_perm(i);
1124  const ordinal_type col_beg = row_ptr(row);
1125  const ordinal_type col_end = row_ptr(row + 1);
1126  const ordinal_type nk = col_end - col_beg;
1127  for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1128  const ordinal_type tk = kbeg + t;
1129  const ordinal_type sk = col_beg + k;
1130  const ordinal_type j = device_peri(cols(sk));
1131  if(i >= j) {
1132  new_cols(tk) = j;
1133  ++t;
1134  }
1135  }
1136  });
1137 
1138  // finally set the inputs to the new sorted arrays
1139  row_ptr = new_row_ptr;
1140  cols = new_cols;
1141 
1142  nnz = new_nnz;
1143  #endif // HAVE_AMESOS2_METIS
1144  }
1145 
1146  template<class values_view_t, class row_ptr_view_t,
1147  class cols_view_t, class per_view_t>
1148  void
1149  reorder_values(values_view_t & values, const row_ptr_view_t & orig_row_ptr,
1150  const row_ptr_view_t & new_row_ptr,
1151  const cols_view_t & orig_cols, const per_view_t & perm, const per_view_t & peri,
1152  size_t nnz)
1153  {
1154  typedef typename cols_view_t::value_type ordinal_type;
1155  typedef typename cols_view_t::execution_space exec_space_t;
1156 
1157  auto device_perm = Kokkos::create_mirror_view(exec_space_t(), perm);
1158  auto device_peri = Kokkos::create_mirror_view(exec_space_t(), peri);
1159  deep_copy(device_perm, perm);
1160  deep_copy(device_peri, peri);
1161 
1162  const ordinal_type size = orig_row_ptr.size() - 1;
1163 
1164  auto host_orig_row_ptr = Kokkos::create_mirror_view(orig_row_ptr);
1165  auto new_nnz = host_orig_row_ptr(size); // TODO: Maybe optimize this by caching
1166 
1167  values_view_t new_values(
1168  Kokkos::ViewAllocateWithoutInitializing("new_values"), values.size() - new_nnz/2);
1169 
1170  // permute col indices
1171  Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1172  Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1173  const ordinal_type kbeg = new_row_ptr(i);
1174  const ordinal_type row = device_perm(i);
1175  const ordinal_type col_beg = orig_row_ptr(row);
1176  const ordinal_type col_end = orig_row_ptr(row + 1);
1177  const ordinal_type nk = col_end - col_beg;
1178  for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1179  const ordinal_type tk = kbeg + t;
1180  const ordinal_type sk = col_beg + k;
1181  const ordinal_type j = device_peri(orig_cols(sk));
1182  if(i >= j) {
1183  new_values(tk) = values(sk);
1184  ++t;
1185  }
1186  }
1187  });
1188 
1189  values = new_values;
1190  }
1191 
1192  template<class array_view_t, class per_view_t>
1193  void
1194  apply_reorder_permutation(const array_view_t & array,
1195  array_view_t & permuted_array, const per_view_t & permutation) {
1196  if(permuted_array.extent(0) != array.extent(0) || permuted_array.extent(1) != array.extent(1)) {
1197  permuted_array = array_view_t(
1198  Kokkos::ViewAllocateWithoutInitializing("permuted_array"),
1199  array.extent(0), array.extent(1));
1200  }
1201  typedef typename array_view_t::execution_space exec_space_t;
1202  Kokkos::RangePolicy<exec_space_t> policy(0, array.extent(0));
1203  Kokkos::parallel_for(policy, KOKKOS_LAMBDA(size_t i) {
1204  for(size_t j = 0; j < array.extent(1); ++j) {
1205  permuted_array(i, j) = array(permutation(i), j);
1206  }
1207  });
1208  }
1209 
1212  } // end namespace Util
1213 
1214 } // end namespace Amesos2
1215 
1216 #endif // #ifndef AMESOS2_UTIL_HPP
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:780
Provides some simple meta-programming utilities for Amesos2.
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:452
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
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.
static void do_get(const Teuchos::Ptr< const Matrix > mat, const ArrayView< S > nzvals, const ArrayView< GO > indices, const ArrayView< GS > pointers, GS &nnz, const Teuchos::Ptr< const Tpetra::Map< typename Matrix::local_ordinal_t, typename Matrix::global_ordinal_t, typename Matrix::node_t > > map, EDistribution distribution, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:502
const RCP< const Teuchos::Comm< int > > to_teuchos_comm(RCP< const Epetra_Comm > c)
Transform an Epetra_Comm object into a Teuchos::Comm object.
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 &quot;-&quot;s on std::cout.
Definition: Amesos2_Util.cpp:119
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:765
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:796
static void do_get(const Teuchos::Ptr< const Matrix > mat, const ArrayView< S > nzvals, const ArrayView< GO > indices, const ArrayView< GS > pointers, GS &nnz, EDistribution distribution, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:483
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:884
Enum and other types declarations for Amesos2.
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:851
EDistribution
Definition: Amesos2_TypeDecl.hpp:123