52 #ifndef AMESOS2_UTIL_HPP 53 #define AMESOS2_UTIL_HPP 55 #include "Amesos2_config.h" 57 #include <Teuchos_RCP.hpp> 58 #include <Teuchos_BLAS_types.hpp> 59 #include <Teuchos_ArrayView.hpp> 60 #include <Teuchos_FancyOStream.hpp> 62 #include <Tpetra_Map.hpp> 63 #include <Tpetra_DistObject_decl.hpp> 64 #include <Tpetra_ComputeGatherMap.hpp> 69 #ifdef HAVE_TPETRA_INST_INT_INT 70 #ifdef HAVE_AMESOS2_EPETRA 71 #include <Epetra_Map.h> 87 using Teuchos::ArrayView;
90 using Meta::if_then_else;
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 );
113 template <
typename LO,
typename GO,
typename GS,
typename Node>
114 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
116 GS num_global_elements,
117 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
119 const Teuchos::RCP<
const Tpetra::Map<LO,GO,Node> >& map = Teuchos::null);
122 #ifdef HAVE_TPETRA_INST_INT_INT 123 #ifdef HAVE_AMESOS2_EPETRA 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);
139 template <
typename LO,
typename GO,
typename GS,
typename Node>
141 tpetra_map_to_epetra_map(
const Tpetra::Map<LO,GO,Node>& map);
148 const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
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 164 template <
typename Scalar,
165 typename GlobalOrdinal,
166 typename GlobalSizeT>
168 ArrayView<GlobalOrdinal> indices,
169 ArrayView<GlobalSizeT> ptr,
170 ArrayView<Scalar> trans_vals,
171 ArrayView<GlobalOrdinal> trans_indices,
172 ArrayView<GlobalSizeT> trans_ptr);
187 template <
typename Scalar1,
typename Scalar2>
188 void scale(ArrayView<Scalar1> vals,
size_t l,
189 size_t ld, ArrayView<Scalar2> s);
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);
215 void printLine( Teuchos::FancyOStream &out );
220 template <
class T0,
class T1 >
221 struct getStdCplxType
223 using common_type =
typename std::common_type<T0,T1>::type;
224 using type = common_type;
227 template <
class T0,
class T1 >
228 struct getStdCplxType< T0, T1* >
230 using common_type =
typename std::common_type<T0,T1>::type;
231 using type = common_type;
234 #if defined(HAVE_TEUCHOS_COMPLEX) && defined(HAVE_AMESOS2_KOKKOS) 235 template <
class T0 >
236 struct getStdCplxType< T0, Kokkos::complex<T0>* >
238 using type = std::complex<T0>;
241 template <
class T0 ,
class T1 >
242 struct getStdCplxType< T0, Kokkos::complex<T1>* >
244 using common_type =
typename std::common_type<T0,T1>::type;
245 using type = std::complex<common_type>;
253 #ifndef DOXYGEN_SHOULD_SKIP_THIS 276 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
277 struct same_gs_helper
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,
285 const Tpetra::Map<
typename M::local_ordinal_t,
286 typename M::global_ordinal_t,
287 typename M::node_t> > map,
291 Op::apply(mat, nzvals, indices, pointers, nnz, map, distribution, ordering);
295 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
296 struct diff_gs_helper
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,
304 const Tpetra::Map<
typename M::local_ordinal_t,
305 typename M::global_ordinal_t,
306 typename M::node_t> > map,
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;
315 Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, distribution, ordering);
317 for (i = 0; i < size; ++i){
318 pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
320 nnz = Teuchos::as<GS>(nnz_tmp);
324 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
325 struct same_go_helper
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,
333 const Tpetra::Map<
typename M::local_ordinal_t,
334 typename M::global_ordinal_t,
335 typename M::node_t> > map,
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,
344 distribution, ordering);
348 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
349 struct diff_go_helper
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,
357 const Tpetra::Map<
typename M::local_ordinal_t,
358 typename M::global_ordinal_t,
359 typename M::node_t> > map,
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);
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,
372 distribution, ordering);
374 for (i = 0; i < size; ++i){
375 indices[i] = Teuchos::as<GO>(indices_tmp[i]);
380 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
381 struct same_scalar_helper
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,
389 const Tpetra::Map<
typename M::local_ordinal_t,
390 typename M::global_ordinal_t,
391 typename M::node_t> > map,
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,
400 distribution, ordering);
404 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
405 struct diff_scalar_helper
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,
413 const Tpetra::Map<
typename M::local_ordinal_t,
414 typename M::global_ordinal_t,
415 typename M::node_t> > map,
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);
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,
428 distribution, ordering);
430 for (i = 0; i < size; ++i){
431 nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
435 #endif // DOXYGEN_SHOULD_SKIP_THIS 449 template<
class Matrix,
typename S,
typename GO,
typename GS,
class Op>
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,
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;
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),
471 Op::getMapFromMatrix(mat)
474 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
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,
489 const Teuchos::RCP<
const Tpetra::Map<
typename Matrix::local_ordinal_t,
490 typename Matrix::global_ordinal_t,
491 typename Matrix::node_t> > map
493 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
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,
506 const Tpetra::Map<
typename Matrix::local_ordinal_t,
507 typename Matrix::global_ordinal_t,
508 typename Matrix::node_t> > map,
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,
519 distribution, ordering);
523 #ifndef DOXYGEN_SHOULD_SKIP_THIS 528 template<
class Matrix>
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,
537 const Tpetra::Map<
typename Matrix::local_ordinal_t,
538 typename Matrix::global_ordinal_t,
539 typename Matrix::node_t> > map,
543 mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering, distribution);
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)
553 return mat->getMap();
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)
562 return mat->getColMap();
566 typename Matrix::global_size_t
567 get_dimension(
const Teuchos::Ptr<const Matrix> mat)
569 return mat->getGlobalNumCols();
573 template<
class Matrix>
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,
582 const Tpetra::Map<
typename Matrix::local_ordinal_t,
583 typename Matrix::global_ordinal_t,
584 typename Matrix::node_t> > map,
588 mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering, distribution);
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)
598 return mat->getMap();
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)
607 return mat->getRowMap();
611 typename Matrix::global_size_t
612 get_dimension(
const Teuchos::Ptr<const Matrix> mat)
614 return mat->getGlobalNumRows();
617 #endif // DOXYGEN_SHOULD_SKIP_THIS 656 template<
class Matrix,
typename S,
typename GO,
typename GS>
667 template<
class Matrix,
typename S,
typename GO,
typename GS>
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 )
684 Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
689 template <
typename LO,
typename GO,
typename GS,
typename Node>
690 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
692 GS num_global_elements,
693 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
695 const Teuchos::RCP<
const Tpetra::Map<LO,GO,Node> >& map)
699 switch( distribution ){
702 return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
704 return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
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; }
711 return rcp(
new Tpetra::Map<LO,GO, Node>(num_global_elements,
712 my_num_elems, indexBase, comm));
716 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
717 = getGatherMap<LO,GO,GS,Node>( map );
721 TEUCHOS_TEST_FOR_EXCEPTION(
true,
723 "Control should never reach this point. " 724 "Please contact the Amesos2 developers." );
729 #ifdef HAVE_TPETRA_INST_INT_INT 730 #ifdef HAVE_AMESOS2_EPETRA 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)
742 int num_my_elements = map.NumMyElements();
743 Teuchos::Array<int> my_global_elements(num_my_elements);
744 map.MyGlobalElements(my_global_elements.getRawPtr());
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()))));
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)
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]);
769 RCP<Epetra_Map> emap = rcp(
new Epetra_Map(-1,
771 my_global_elements.getRawPtr(),
772 as<GO>(map.getIndexBase()),
773 *to_epetra_comm(map.getComm())));
776 #endif // HAVE_AMESOS2_EPETRA 777 #endif // HAVE_TPETRA_INST_INT_INT 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)
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." );
810 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
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]);
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);
825 trans_ptr.assign(count);
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]]);
853 template <
typename Scalar1,
typename Scalar2>
855 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
856 size_t ld, Teuchos::ArrayView<Scalar2> s)
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" );
866 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
876 template <
typename Scalar1,
typename Scalar2,
class BinaryOp>
878 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
879 size_t ld, Teuchos::ArrayView<Scalar2> s,
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" );
890 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
896 vals[i] = binary_op(vals[i], s[s_i]);
906 #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:668
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:450
Definition: Amesos2_TypeDecl.hpp:143
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
Definition: Amesos2_TypeDecl.hpp:125
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
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
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition: Amesos2_Util.cpp:123
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:657
Definition: Amesos2_TypeDecl.hpp:126
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
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
Definition: Amesos2_TypeDecl.hpp:124
Enum and other types declarations for Amesos2.
Definition: Amesos2_TypeDecl.hpp:127
EDistribution
Definition: Amesos2_TypeDecl.hpp:123
Definition: Amesos2_TypeDecl.hpp:128