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