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_ArrayView.hpp>
60 #include <Teuchos_FancyOStream.hpp>
61 
62 #include <Tpetra_Map.hpp>
63 #include <Tpetra_DistObject_decl.hpp>
64 #include <Tpetra_ComputeGatherMap.hpp> // added for gather map... where is the best place??
65 
66 #include "Amesos2_TypeDecl.hpp"
67 #include "Amesos2_Meta.hpp"
68 
69 #ifdef HAVE_TPETRA_INST_INT_INT
70 #ifdef HAVE_AMESOS2_EPETRA
71 #include <Epetra_Map.h>
72 #endif
73 #endif
74 
75 
76 namespace Amesos2 {
77 
78  namespace Util {
79 
86  using Teuchos::RCP;
87  using Teuchos::ArrayView;
88 
89  using Meta::is_same;
90  using Meta::if_then_else;
91 
108  template <typename LO, typename GO, typename GS, typename Node>
109  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
110  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map );
111 
112 
113  template <typename LO, typename GO, typename GS, typename Node>
114  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
115  getDistributionMap(EDistribution distribution,
116  GS num_global_elements,
117  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
118  GO indexBase = 0,
119  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map = Teuchos::null);
120 
121 
122 #ifdef HAVE_TPETRA_INST_INT_INT
123 #ifdef HAVE_AMESOS2_EPETRA
124 
130  template <typename LO, typename GO, typename GS, typename Node>
131  RCP<Tpetra::Map<LO,GO,Node> >
132  epetra_map_to_tpetra_map(const Epetra_BlockMap& map);
133 
139  template <typename LO, typename GO, typename GS, typename Node>
140  RCP<Epetra_Map>
141  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map);
142 
148  const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
149 
155  const RCP<const Epetra_Comm> to_epetra_comm(RCP<const Teuchos::Comm<int> > c);
156 #endif // HAVE_AMESOS2_EPETRA
157 #endif // HAVE_TPETRA_INST_INT_INT
158 
164  template <typename Scalar,
165  typename GlobalOrdinal,
166  typename GlobalSizeT>
167  void transpose(ArrayView<Scalar> vals,
168  ArrayView<GlobalOrdinal> indices,
169  ArrayView<GlobalSizeT> ptr,
170  ArrayView<Scalar> trans_vals,
171  ArrayView<GlobalOrdinal> trans_indices,
172  ArrayView<GlobalSizeT> trans_ptr);
173 
187  template <typename Scalar1, typename Scalar2>
188  void scale(ArrayView<Scalar1> vals, size_t l,
189  size_t ld, ArrayView<Scalar2> s);
190 
209  template <typename Scalar1, typename Scalar2, class BinaryOp>
210  void scale(ArrayView<Scalar1> vals, size_t l,
211  size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
212 
213 
215  void printLine( Teuchos::FancyOStream &out );
216 
217  // Helper function used to convert Kokkos::complex pointer
218  // to std::complex pointer; needed for optimized code path
219  // when retrieving the CRS raw pointers
220  template < class T0, class T1 >
221  struct getStdCplxType
222  {
223  using common_type = typename std::common_type<T0,T1>::type;
224  using type = common_type;
225  };
226 
227  template < class T0, class T1 >
228  struct getStdCplxType< T0, T1* >
229  {
230  using common_type = typename std::common_type<T0,T1>::type;
231  using type = common_type;
232  };
233 
234 #if defined(HAVE_TEUCHOS_COMPLEX) && defined(HAVE_AMESOS2_KOKKOS)
235  template < class T0 >
236  struct getStdCplxType< T0, Kokkos::complex<T0>* >
237  {
238  using type = std::complex<T0>;
239  };
240 
241  template < class T0 , class T1 >
242  struct getStdCplxType< T0, Kokkos::complex<T1>* >
243  {
244  using common_type = typename std::common_type<T0,T1>::type;
245  using type = std::complex<common_type>;
246  };
247 #endif
248 
250  // Matrix/MultiVector Utilities //
252 
253 #ifndef DOXYGEN_SHOULD_SKIP_THIS
254  /*
255  * The following represents a general way of getting a CRS or CCS
256  * representation of a matrix with implicit type conversions. The
257  * \c get_crs_helper and \c get_ccs_helper classes are templated
258  * on 4 types:
259  *
260  * - A Matrix type (conforming to the Amesos2 MatrixAdapter interface)
261  * - A scalar type
262  * - A global ordinal type
263  * - A global size type
264  *
265  * The last three template types correspond to the input argument
266  * types. For example, if the scalar type is \c double , then we
267  * require that the \c nzvals argument is a \c
268  * Teuchos::ArrayView<double> type.
269  *
270  * These helpers perform any type conversions that must be
271  * performed to go between the Matrix's types and the input types.
272  * If no conversions are necessary the static functions can be
273  * effectively inlined.
274  */
275 
276  template <class M, typename S, typename GO, typename GS, class Op>
277  struct same_gs_helper
278  {
279  static void do_get(const Teuchos::Ptr<const M> mat,
280  const ArrayView<typename M::scalar_t> nzvals,
281  const ArrayView<typename M::global_ordinal_t> indices,
282  const ArrayView<GS> pointers,
283  GS& nnz,
284  const Teuchos::Ptr<
285  const Tpetra::Map<typename M::local_ordinal_t,
286  typename M::global_ordinal_t,
287  typename M::node_t> > map,
288  EDistribution distribution,
289  EStorage_Ordering ordering)
290  {
291  Op::apply(mat, nzvals, indices, pointers, nnz, map, distribution, ordering);
292  }
293  };
294 
295  template <class M, typename S, typename GO, typename GS, class Op>
296  struct diff_gs_helper
297  {
298  static void do_get(const Teuchos::Ptr<const M> mat,
299  const ArrayView<typename M::scalar_t> nzvals,
300  const ArrayView<typename M::global_ordinal_t> indices,
301  const ArrayView<GS> pointers,
302  GS& nnz,
303  const Teuchos::Ptr<
304  const Tpetra::Map<typename M::local_ordinal_t,
305  typename M::global_ordinal_t,
306  typename M::node_t> > map,
307  EDistribution distribution,
308  EStorage_Ordering ordering)
309  {
310  typedef typename M::global_size_t mat_gs_t;
311  typename ArrayView<GS>::size_type i, size = pointers.size();
312  Teuchos::Array<mat_gs_t> pointers_tmp(size);
313  mat_gs_t nnz_tmp = 0;
314 
315  Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, distribution, ordering);
316 
317  for (i = 0; i < size; ++i){
318  pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
319  }
320  nnz = Teuchos::as<GS>(nnz_tmp);
321  }
322  };
323 
324  template <class M, typename S, typename GO, typename GS, class Op>
325  struct same_go_helper
326  {
327  static void do_get(const Teuchos::Ptr<const M> mat,
328  const ArrayView<typename M::scalar_t> nzvals,
329  const ArrayView<GO> indices,
330  const ArrayView<GS> pointers,
331  GS& nnz,
332  const Teuchos::Ptr<
333  const Tpetra::Map<typename M::local_ordinal_t,
334  typename M::global_ordinal_t,
335  typename M::node_t> > map,
336  EDistribution distribution,
337  EStorage_Ordering ordering)
338  {
339  typedef typename M::global_size_t mat_gs_t;
340  if_then_else<is_same<GS,mat_gs_t>::value,
341  same_gs_helper<M,S,GO,GS,Op>,
342  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
343  pointers, nnz, map,
344  distribution, ordering);
345  }
346  };
347 
348  template <class M, typename S, typename GO, typename GS, class Op>
349  struct diff_go_helper
350  {
351  static void do_get(const Teuchos::Ptr<const M> mat,
352  const ArrayView<typename M::scalar_t> nzvals,
353  const ArrayView<GO> indices,
354  const ArrayView<GS> pointers,
355  GS& nnz,
356  const Teuchos::Ptr<
357  const Tpetra::Map<typename M::local_ordinal_t,
358  typename M::global_ordinal_t,
359  typename M::node_t> > map,
360  EDistribution distribution,
361  EStorage_Ordering ordering)
362  {
363  typedef typename M::global_ordinal_t mat_go_t;
364  typedef typename M::global_size_t mat_gs_t;
365  typename ArrayView<GO>::size_type i, size = indices.size();
366  Teuchos::Array<mat_go_t> indices_tmp(size);
367 
368  if_then_else<is_same<GS,mat_gs_t>::value,
369  same_gs_helper<M,S,GO,GS,Op>,
370  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
371  pointers, nnz, map,
372  distribution, ordering);
373 
374  for (i = 0; i < size; ++i){
375  indices[i] = Teuchos::as<GO>(indices_tmp[i]);
376  }
377  }
378  };
379 
380  template <class M, typename S, typename GO, typename GS, class Op>
381  struct same_scalar_helper
382  {
383  static void do_get(const Teuchos::Ptr<const M> mat,
384  const ArrayView<S> nzvals,
385  const ArrayView<GO> indices,
386  const ArrayView<GS> pointers,
387  GS& nnz,
388  const Teuchos::Ptr<
389  const Tpetra::Map<typename M::local_ordinal_t,
390  typename M::global_ordinal_t,
391  typename M::node_t> > map,
392  EDistribution distribution,
393  EStorage_Ordering ordering)
394  {
395  typedef typename M::global_ordinal_t mat_go_t;
396  if_then_else<is_same<GO,mat_go_t>::value,
397  same_go_helper<M,S,GO,GS,Op>,
398  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
399  pointers, nnz, map,
400  distribution, ordering);
401  }
402  };
403 
404  template <class M, typename S, typename GO, typename GS, class Op>
405  struct diff_scalar_helper
406  {
407  static void do_get(const Teuchos::Ptr<const M> mat,
408  const ArrayView<S> nzvals,
409  const ArrayView<GO> indices,
410  const ArrayView<GS> pointers,
411  GS& nnz,
412  const Teuchos::Ptr<
413  const Tpetra::Map<typename M::local_ordinal_t,
414  typename M::global_ordinal_t,
415  typename M::node_t> > map,
416  EDistribution distribution,
417  EStorage_Ordering ordering)
418  {
419  typedef typename M::scalar_t mat_scalar_t;
420  typedef typename M::global_ordinal_t mat_go_t;
421  typename ArrayView<S>::size_type i, size = nzvals.size();
422  Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
423 
424  if_then_else<is_same<GO,mat_go_t>::value,
425  same_go_helper<M,S,GO,GS,Op>,
426  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
427  pointers, nnz, map,
428  distribution, ordering);
429 
430  for (i = 0; i < size; ++i){
431  nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
432  }
433  }
434  };
435 #endif // DOXYGEN_SHOULD_SKIP_THIS
436 
449  template<class Matrix, typename S, typename GO, typename GS, class Op>
450  struct get_cxs_helper
451  {
452  static void do_get(const Teuchos::Ptr<const Matrix> mat,
453  const ArrayView<S> nzvals,
454  const ArrayView<GO> indices,
455  const ArrayView<GS> pointers,
456  GS& nnz,
457  EDistribution distribution,
458  EStorage_Ordering ordering=ARBITRARY,
459  GO indexBase = 0)
460  {
461  typedef typename Matrix::local_ordinal_t lo_t;
462  typedef typename Matrix::global_ordinal_t go_t;
463  typedef typename Matrix::global_size_t gs_t;
464  typedef typename Matrix::node_t node_t;
465 
466  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
467  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
468  Op::get_dimension(mat),
469  mat->getComm(),
470  indexBase,
471  Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
472  );
473 
474  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
475  }
476 
481  static void do_get(const Teuchos::Ptr<const Matrix> mat,
482  const ArrayView<S> nzvals,
483  const ArrayView<GO> indices,
484  const ArrayView<GS> pointers,
485  GS& nnz,
486  EDistribution distribution, // Does this one need a distribution argument??
487  EStorage_Ordering ordering=ARBITRARY)
488  {
489  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
490  typename Matrix::global_ordinal_t,
491  typename Matrix::node_t> > map
492  = Op::getMap(mat);
493  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
494  }
495 
500  static void do_get(const Teuchos::Ptr<const Matrix> mat,
501  const ArrayView<S> nzvals,
502  const ArrayView<GO> indices,
503  const ArrayView<GS> pointers,
504  GS& nnz,
505  const Teuchos::Ptr<
506  const Tpetra::Map<typename Matrix::local_ordinal_t,
507  typename Matrix::global_ordinal_t,
508  typename Matrix::node_t> > map,
509  EDistribution distribution,
510  EStorage_Ordering ordering=ARBITRARY)
511  {
512  typedef typename Matrix::scalar_t mat_scalar;
513  if_then_else<is_same<mat_scalar,S>::value,
514  same_scalar_helper<Matrix,S,GO,GS,Op>,
515  diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
516  nzvals, indices,
517  pointers, nnz,
518  map,
519  distribution, ordering);
520  }
521  };
522 
523 #ifndef DOXYGEN_SHOULD_SKIP_THIS
524  /*
525  * These two function-like classes are meant to be used as the \c
526  * Op template parameter for the \c get_cxs_helper template class.
527  */
528  template<class Matrix>
529  struct get_ccs_func
530  {
531  static void apply(const Teuchos::Ptr<const Matrix> mat,
532  const ArrayView<typename Matrix::scalar_t> nzvals,
533  const ArrayView<typename Matrix::global_ordinal_t> rowind,
534  const ArrayView<typename Matrix::global_size_t> colptr,
535  typename Matrix::global_size_t& nnz,
536  const Teuchos::Ptr<
537  const Tpetra::Map<typename Matrix::local_ordinal_t,
538  typename Matrix::global_ordinal_t,
539  typename Matrix::node_t> > map,
540  EDistribution distribution,
541  EStorage_Ordering ordering)
542  {
543  mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering, distribution);
544  //mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
545  }
546 
547  static
548  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
549  typename Matrix::global_ordinal_t,
550  typename Matrix::node_t> >
551  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
552  {
553  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
554  }
555 
556  static
557  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
558  typename Matrix::global_ordinal_t,
559  typename Matrix::node_t> >
560  getMap(const Teuchos::Ptr<const Matrix> mat)
561  {
562  return mat->getColMap();
563  }
564 
565  static
566  typename Matrix::global_size_t
567  get_dimension(const Teuchos::Ptr<const Matrix> mat)
568  {
569  return mat->getGlobalNumCols();
570  }
571  };
572 
573  template<class Matrix>
574  struct get_crs_func
575  {
576  static void apply(const Teuchos::Ptr<const Matrix> mat,
577  const ArrayView<typename Matrix::scalar_t> nzvals,
578  const ArrayView<typename Matrix::global_ordinal_t> colind,
579  const ArrayView<typename Matrix::global_size_t> rowptr,
580  typename Matrix::global_size_t& nnz,
581  const Teuchos::Ptr<
582  const Tpetra::Map<typename Matrix::local_ordinal_t,
583  typename Matrix::global_ordinal_t,
584  typename Matrix::node_t> > map,
585  EDistribution distribution,
586  EStorage_Ordering ordering)
587  {
588  mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering, distribution);
589  //mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering);
590  }
591 
592  static
593  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
594  typename Matrix::global_ordinal_t,
595  typename Matrix::node_t> >
596  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
597  {
598  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
599  }
600 
601  static
602  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
603  typename Matrix::global_ordinal_t,
604  typename Matrix::node_t> >
605  getMap(const Teuchos::Ptr<const Matrix> mat)
606  {
607  return mat->getRowMap();
608  }
609 
610  static
611  typename Matrix::global_size_t
612  get_dimension(const Teuchos::Ptr<const Matrix> mat)
613  {
614  return mat->getGlobalNumRows();
615  }
616  };
617 #endif // DOXYGEN_SHOULD_SKIP_THIS
618 
656  template<class Matrix, typename S, typename GO, typename GS>
657  struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
658  {};
659 
667  template<class Matrix, typename S, typename GO, typename GS>
668  struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
669  {};
670 
671  /* End Matrix/MultiVector Utilities */
672 
673 
675  // Definitions //
677 
678 
679  template <typename LO, typename GO, typename GS, typename Node>
680  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
681  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
682  {
683  //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
684  Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
685  return gather_map;
686  }
687 
688 
689  template <typename LO, typename GO, typename GS, typename Node>
690  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
691  getDistributionMap(EDistribution distribution,
692  GS num_global_elements,
693  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
694  GO indexBase,
695  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
696  {
697  // TODO: Need to add indexBase to cases other than ROOTED
698  // We do not support these maps in any solver now.
699  switch( distribution ){
700  case DISTRIBUTED:
701  case DISTRIBUTED_NO_OVERLAP:
702  return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
703  case GLOBALLY_REPLICATED:
704  return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
705  case ROOTED:
706  {
707  int rank = Teuchos::rank(*comm);
708  size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
709  if( rank == 0 ) { my_num_elems = num_global_elements; }
710 
711  return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
712  my_num_elems, indexBase, comm));
713  }
714  case CONTIGUOUS_AND_ROOTED:
715  {
716  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
717  = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
718  return gathermap;
719  }
720  default:
721  TEUCHOS_TEST_FOR_EXCEPTION( true,
722  std::logic_error,
723  "Control should never reach this point. "
724  "Please contact the Amesos2 developers." );
725  }
726  }
727 
728 
729 #ifdef HAVE_TPETRA_INST_INT_INT
730 #ifdef HAVE_AMESOS2_EPETRA
731 
732  //#pragma message "include 3"
733  //#include <Epetra_Map.h>
734 
735  template <typename LO, typename GO, typename GS, typename Node>
736  Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
737  epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
738  {
739  using Teuchos::as;
740  using Teuchos::rcp;
741 
742  int num_my_elements = map.NumMyElements();
743  Teuchos::Array<int> my_global_elements(num_my_elements);
744  map.MyGlobalElements(my_global_elements.getRawPtr());
745 
746  typedef Tpetra::Map<LO,GO,Node> map_t;
747  RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
748  my_global_elements(),
749  as<GO>(map.IndexBase()),
750  to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
751  return tmap;
752  }
753 
754  template <typename LO, typename GO, typename GS, typename Node>
755  Teuchos::RCP<Epetra_Map>
756  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
757  {
758  using Teuchos::as;
759 
760  Teuchos::Array<GO> elements_tmp;
761  elements_tmp = map.getNodeElementList();
762  int num_my_elements = elements_tmp.size();
763  Teuchos::Array<int> my_global_elements(num_my_elements);
764  for (int i = 0; i < num_my_elements; ++i){
765  my_global_elements[i] = as<int>(elements_tmp[i]);
766  }
767 
768  using Teuchos::rcp;
769  RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
770  num_my_elements,
771  my_global_elements.getRawPtr(),
772  as<GO>(map.getIndexBase()),
773  *to_epetra_comm(map.getComm())));
774  return emap;
775  }
776 #endif // HAVE_AMESOS2_EPETRA
777 #endif // HAVE_TPETRA_INST_INT_INT
778 
779  template <typename Scalar,
780  typename GlobalOrdinal,
781  typename GlobalSizeT>
782  void transpose(Teuchos::ArrayView<Scalar> vals,
783  Teuchos::ArrayView<GlobalOrdinal> indices,
784  Teuchos::ArrayView<GlobalSizeT> ptr,
785  Teuchos::ArrayView<Scalar> trans_vals,
786  Teuchos::ArrayView<GlobalOrdinal> trans_indices,
787  Teuchos::ArrayView<GlobalSizeT> trans_ptr)
788  {
789  /* We have a compressed-row storage format of this matrix. We
790  * transform this into a compressed-column format using a
791  * distribution-counting sort algorithm, which is described by
792  * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
793  */
794 
795 #ifdef HAVE_AMESOS2_DEBUG
796  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
797  ind_begin = indices.begin();
798  ind_end = indices.end();
799  size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
800  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
801  std::invalid_argument,
802  "Transpose pointer size not large enough." );
803  TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
804  std::invalid_argument,
805  "Transpose values array not large enough." );
806  TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
807  std::invalid_argument,
808  "Transpose indices array not large enough." );
809 #else
810  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
811 #endif
812  // Count the number of entries in each column
813  Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
814  ind_end = indices.end();
815  for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
816  ++(count[(*ind_it) + 1]);
817  }
818  // Accumulate
819  typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
820  cnt_end = count.end();
821  for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
822  *cnt_it = *cnt_it + *(cnt_it - 1);
823  }
824  // This becomes the array of column pointers
825  trans_ptr.assign(count);
826 
827  /* Move the nonzero values into their final place in nzval, based on the
828  * counts found previously.
829  *
830  * This sequence deviates from Knuth's algorithm a bit, following more
831  * closely the description presented in Gustavson, Fred G. "Two Fast
832  * Algorithms for Sparse Matrices: Multiplication and Permuted
833  * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
834  * 250--269, http://doi.acm.org/10.1145/355791.355796.
835  *
836  * The output indices end up in sorted order
837  */
838 
839  GlobalSizeT size = ptr.size();
840  for( GlobalSizeT i = 0; i < size - 1; ++i ){
841  GlobalOrdinal u = ptr[i];
842  GlobalOrdinal v = ptr[i + 1];
843  for( GlobalOrdinal j = u; j < v; ++j ){
844  GlobalOrdinal k = count[indices[j]];
845  trans_vals[k] = vals[j];
846  trans_indices[k] = i;
847  ++(count[indices[j]]);
848  }
849  }
850  }
851 
852 
853  template <typename Scalar1, typename Scalar2>
854  void
855  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
856  size_t ld, Teuchos::ArrayView<Scalar2> s)
857  {
858  size_t vals_size = vals.size();
859 #ifdef HAVE_AMESOS2_DEBUG
860  size_t s_size = s.size();
861  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
862  std::invalid_argument,
863  "Scale vector must have length at least that of the vector" );
864 #endif
865  size_t i, s_i;
866  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
867  if( s_i == l ){
868  // bring i to the next multiple of ld
869  i += ld - s_i;
870  s_i = 0;
871  }
872  vals[i] *= s[s_i];
873  }
874  }
875 
876  template <typename Scalar1, typename Scalar2, class BinaryOp>
877  void
878  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
879  size_t ld, Teuchos::ArrayView<Scalar2> s,
880  BinaryOp binary_op)
881  {
882  size_t vals_size = vals.size();
883 #ifdef HAVE_AMESOS2_DEBUG
884  size_t s_size = s.size();
885  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
886  std::invalid_argument,
887  "Scale vector must have length at least that of the vector" );
888 #endif
889  size_t i, s_i;
890  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
891  if( s_i == l ){
892  // bring i to the next multiple of ld
893  i += ld - s_i;
894  s_i = 0;
895  }
896  vals[i] = binary_op(vals[i], s[s_i]);
897  }
898  }
899 
902  } // end namespace Util
903 
904 } // end namespace Amesos2
905 
906 #endif // #ifndef AMESOS2_UTIL_HPP
Amesos2::Util::get_cxs_helper
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:450
Amesos2::EDistribution
EDistribution
Definition: Amesos2_TypeDecl.hpp:123
Amesos2::Util::get_cxs_helper::do_get
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:481
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:681
Amesos2::EStorage_Ordering
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
Amesos2_Util.hpp
Utility functions for Amesos2.
Amesos2::Util::get_cxs_helper::do_get
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:500
Amesos2_Meta.hpp
Provides some simple meta-programming utilities for Amesos2.
Amesos2::Util::get_ccs_helper
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:657
Amesos2::DISTRIBUTED_NO_OVERLAP
Definition: Amesos2_TypeDecl.hpp:125
Amesos2_TypeDecl.hpp
Enum and other types declarations for Amesos2.
Amesos2::ARBITRARY
Definition: Amesos2_TypeDecl.hpp:143
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::CONTIGUOUS_AND_ROOTED
Definition: Amesos2_TypeDecl.hpp:128
Amesos2::ROOTED
Definition: Amesos2_TypeDecl.hpp:127
Amesos2::GLOBALLY_REPLICATED
Definition: Amesos2_TypeDecl.hpp:126
Amesos2::DISTRIBUTED
Definition: Amesos2_TypeDecl.hpp:124
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::Util::printLine
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition: Amesos2_Util.cpp:123
Amesos2::Util::get_crs_helper
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:668
Generated by   doxygen 1.8.16