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_Array.hpp"
60 #include "Teuchos_ArrayView.hpp"
61 #include "Teuchos_FancyOStream.hpp"
63 #include <Tpetra_Map.hpp>
64 #include <Tpetra_DistObject_decl.hpp>
65 #include <Tpetra_ComputeGatherMap.hpp>
71 #ifdef HAVE_AMESOS2_EPETRA
72 #include <Epetra_Map.h>
75 #ifdef HAVE_AMESOS2_METIS
90 using Teuchos::ArrayView;
93 using Meta::if_then_else;
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 );
116 template <
typename LO,
typename GO,
typename GS,
typename Node>
117 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
119 GS num_global_elements,
120 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
122 const Teuchos::RCP<
const Tpetra::Map<LO,GO,Node> >& map = Teuchos::null);
125 #ifdef HAVE_AMESOS2_EPETRA
132 template <
typename LO,
typename GO,
typename GS,
typename Node>
133 RCP<Tpetra::Map<LO,GO,Node> >
141 template <
typename LO,
typename GO,
typename GS,
typename Node>
150 const RCP<const Teuchos::Comm<int> >
to_teuchos_comm(RCP<const Epetra_Comm> c);
157 const RCP<const Epetra_Comm>
to_epetra_comm(RCP<
const Teuchos::Comm<int> > c);
158 #endif // HAVE_AMESOS2_EPETRA
165 template <
typename Scalar,
166 typename GlobalOrdinal,
167 typename GlobalSizeT>
169 ArrayView<GlobalOrdinal> indices,
170 ArrayView<GlobalSizeT> ptr,
171 ArrayView<Scalar> trans_vals,
172 ArrayView<GlobalOrdinal> trans_indices,
173 ArrayView<GlobalSizeT> trans_ptr);
188 template <
typename Scalar1,
typename Scalar2>
189 void scale(ArrayView<Scalar1> vals,
size_t l,
190 size_t ld, ArrayView<Scalar2> s);
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);
216 void printLine( Teuchos::FancyOStream &out );
221 template <
class T0,
class T1 >
222 struct getStdCplxType
224 using common_type =
typename std::common_type<T0,T1>::type;
225 using type = common_type;
228 template <
class T0,
class T1 >
229 struct getStdCplxType< T0, T1* >
231 using common_type =
typename std::common_type<T0,T1>::type;
232 using type = common_type;
235 #if defined(HAVE_TEUCHOS_COMPLEX) && defined(HAVE_AMESOS2_KOKKOS)
236 template <
class T0 >
237 struct getStdCplxType< T0, Kokkos::complex<T0>* >
239 using type = std::complex<T0>;
242 template <
class T0 ,
class T1 >
243 struct getStdCplxType< T0, Kokkos::complex<T1>* >
245 using common_type =
typename std::common_type<T0,T1>::type;
246 using type = std::complex<common_type>;
254 #ifndef DOXYGEN_SHOULD_SKIP_THIS
277 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
278 struct same_gs_helper
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,
286 const Tpetra::Map<
typename M::local_ordinal_t,
287 typename M::global_ordinal_t,
288 typename M::node_t> > map,
292 Op::apply(mat, nzvals, indices, pointers, nnz, map, distribution, ordering);
296 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
297 struct diff_gs_helper
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,
305 const Tpetra::Map<
typename M::local_ordinal_t,
306 typename M::global_ordinal_t,
307 typename M::node_t> > map,
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;
316 Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, distribution, ordering);
318 for (i = 0; i < size; ++i){
319 pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
321 nnz = Teuchos::as<GS>(nnz_tmp);
325 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
326 struct same_go_helper
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,
334 const Tpetra::Map<
typename M::local_ordinal_t,
335 typename M::global_ordinal_t,
336 typename M::node_t> > map,
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,
345 distribution, ordering);
349 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
350 struct diff_go_helper
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,
358 const Tpetra::Map<
typename M::local_ordinal_t,
359 typename M::global_ordinal_t,
360 typename M::node_t> > map,
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);
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,
373 distribution, ordering);
375 for (i = 0; i < size; ++i){
376 indices[i] = Teuchos::as<GO>(indices_tmp[i]);
381 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
382 struct same_scalar_helper
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,
390 const Tpetra::Map<
typename M::local_ordinal_t,
391 typename M::global_ordinal_t,
392 typename M::node_t> > map,
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,
401 distribution, ordering);
405 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
406 struct diff_scalar_helper
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,
414 const Tpetra::Map<
typename M::local_ordinal_t,
415 typename M::global_ordinal_t,
416 typename M::node_t> > map,
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);
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,
429 distribution, ordering);
431 for (i = 0; i < size; ++i){
432 nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
437 #endif // DOXYGEN_SHOULD_SKIP_THIS
451 template<
class Matrix,
typename S,
typename GO,
typename GS,
class Op>
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,
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;
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),
473 Op::getMapFromMatrix(mat)
476 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
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,
491 const Teuchos::RCP<
const Tpetra::Map<
typename Matrix::local_ordinal_t,
492 typename Matrix::global_ordinal_t,
493 typename Matrix::node_t> > map
495 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
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,
508 const Tpetra::Map<
typename Matrix::local_ordinal_t,
509 typename Matrix::global_ordinal_t,
510 typename Matrix::node_t> > map,
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,
521 distribution, ordering);
525 template<
class Matrix,
typename KV_S,
typename KV_GO,
typename KV_GS,
class Op>
526 struct get_cxs_helper_kokkos_view
528 static void do_get(
const Teuchos::Ptr<const Matrix> mat,
532 typename KV_GS::value_type& nnz,
535 typename KV_GO::value_type indexBase = 0)
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;
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),
547 Op::getMapFromMatrix(mat)
549 typename Matrix::global_size_t nnz_temp = 0;
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);
559 static void do_get(
const Teuchos::Ptr<const Matrix> mat,
563 typename KV_GS::value_type& nnz,
567 const Teuchos::RCP<
const Tpetra::Map<
typename Matrix::local_ordinal_t,
568 typename Matrix::global_ordinal_t,
569 typename Matrix::node_t> > map
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);
581 static void do_get(
const Teuchos::Ptr<const Matrix> mat,
585 typename KV_GS::value_type& nnz,
587 const Tpetra::Map<
typename Matrix::local_ordinal_t,
588 typename Matrix::global_ordinal_t,
589 typename Matrix::node_t> > map,
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);
600 #ifndef DOXYGEN_SHOULD_SKIP_THIS
605 template<
class Matrix>
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,
614 const Tpetra::Map<
typename Matrix::local_ordinal_t,
615 typename Matrix::global_ordinal_t,
616 typename Matrix::node_t> > map,
620 mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering, distribution);
624 template<
typename KV_S,
typename KV_GO,
typename KV_GS>
625 static void apply_kokkos_view(
const Teuchos::Ptr<const Matrix> mat,
629 typename Matrix::global_size_t& nnz,
631 const Tpetra::Map<
typename Matrix::local_ordinal_t,
632 typename Matrix::global_ordinal_t,
633 typename Matrix::node_t> > map,
637 mat->getCcs_kokkos_view(nzvals, rowind, colptr, nnz, map, ordering, distribution);
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)
646 return mat->getMap();
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)
655 return mat->getColMap();
659 typename Matrix::global_size_t
660 get_dimension(
const Teuchos::Ptr<const Matrix> mat)
662 return mat->getGlobalNumCols();
666 template<
class Matrix>
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,
675 const Tpetra::Map<
typename Matrix::local_ordinal_t,
676 typename Matrix::global_ordinal_t,
677 typename Matrix::node_t> > map,
681 mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering, distribution);
684 template<
typename KV_S,
typename KV_GO,
typename KV_GS>
685 static void apply_kokkos_view(
const Teuchos::Ptr<const Matrix> mat,
689 typename Matrix::global_size_t& nnz,
691 const Tpetra::Map<
typename Matrix::local_ordinal_t,
692 typename Matrix::global_ordinal_t,
693 typename Matrix::node_t> > map,
697 mat->getCrs_kokkos_view(nzvals, colind, rowptr, nnz, map, ordering, distribution);
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)
706 return mat->getMap();
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)
715 return mat->getRowMap();
719 typename Matrix::global_size_t
720 get_dimension(
const Teuchos::Ptr<const Matrix> mat)
722 return mat->getGlobalNumRows();
725 #endif // DOXYGEN_SHOULD_SKIP_THIS
764 template<
class Matrix,
typename S,
typename GO,
typename GS>
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> >
779 template<
class Matrix,
typename S,
typename GO,
typename GS>
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> >
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 )
799 Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
804 template <
typename LO,
typename GO,
typename GS,
typename Node>
805 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
807 GS num_global_elements,
808 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
810 const Teuchos::RCP<
const Tpetra::Map<LO,GO,Node> >& map)
814 switch( distribution ){
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);
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; }
826 return rcp(
new Tpetra::Map<LO,GO, Node>(num_global_elements,
827 my_num_elems, indexBase, comm));
829 case CONTIGUOUS_AND_ROOTED:
831 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
832 = getGatherMap<LO,GO,GS,Node>( map );
836 TEUCHOS_TEST_FOR_EXCEPTION(
true,
838 "Control should never reach this point. "
839 "Please contact the Amesos2 developers." );
844 #ifdef HAVE_AMESOS2_EPETRA
849 template <
typename LO,
typename GO,
typename GS,
typename Node>
850 Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
856 int num_my_elements = map.NumMyElements();
857 Teuchos::Array<int> my_global_elements(num_my_elements);
858 map.MyGlobalElements(my_global_elements.getRawPtr());
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]);
870 using Teuchos::av_reinterpret_cast;
871 my_gbl_inds = av_reinterpret_cast<GO> (my_global_elements ());
874 typedef Tpetra::Map<LO,GO,Node> map_t;
875 RCP<map_t> tmap = rcp(
new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
877 as<GO>(map.IndexBase()),
882 template <
typename LO,
typename GO,
typename GS,
typename Node>
883 Teuchos::RCP<Epetra_Map>
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]);
897 RCP<Epetra_Map> emap = rcp(
new Epetra_Map(-1,
899 my_global_elements.getRawPtr(),
900 as<GO>(map.getIndexBase()),
904 #endif // HAVE_AMESOS2_EPETRA
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)
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." );
937 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
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]);
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);
952 trans_ptr.assign(count);
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]]);
980 template <
typename Scalar1,
typename Scalar2>
982 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
983 size_t ld, Teuchos::ArrayView<Scalar2> s)
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" );
993 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
1003 template <
typename Scalar1,
typename Scalar2,
class BinaryOp>
1005 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
1006 size_t ld, Teuchos::ArrayView<Scalar2> s,
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" );
1017 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
1023 vals[i] = binary_op(vals[i], s[s_i]);
1027 template<
class row_ptr_view_t,
class cols_view_t,
class per_view_t>
1029 reorder(row_ptr_view_t & row_ptr, cols_view_t & cols,
1030 per_view_t & perm, per_view_t & peri,
size_t & nnz)
1032 #ifndef HAVE_AMESOS2_METIS
1033 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
1034 "Cannot reorder for cuSolver because no METIS is available.");
1036 typedef typename cols_view_t::value_type ordinal_type;
1037 typedef typename row_ptr_view_t::value_type size_type;
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);
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"),
1052 host_metis_array host_strip_diag_cols(
1053 Kokkos::ViewAllocateWithoutInitializing(
"host_strip_diag_cols"),
1054 cols.size() - size);
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);
1065 host_strip_diag_row_ptr(size) = new_nnz;
1068 host_metis_array host_perm(
1069 Kokkos::ViewAllocateWithoutInitializing(
"host_perm"), size);
1070 host_metis_array host_peri(
1071 Kokkos::ViewAllocateWithoutInitializing(
"host_peri"), size);
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());
1079 TEUCHOS_TEST_FOR_EXCEPTION(err != METIS_OK, std::runtime_error,
1080 "METIS_NodeND failed to sort matrix.");
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);
1092 deep_copy_or_assign_view(perm, device_perm);
1093 deep_copy_or_assign_view(peri, device_peri);
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);
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) {
1106 new_row_ptr(i) = update;
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));
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));
1139 row_ptr = new_row_ptr;
1143 #endif // HAVE_AMESOS2_METIS
1146 template<
class values_view_t,
class row_ptr_view_t,
1147 class cols_view_t,
class per_view_t>
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,
1154 typedef typename cols_view_t::value_type ordinal_type;
1155 typedef typename cols_view_t::execution_space exec_space_t;
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);
1162 const ordinal_type size = orig_row_ptr.size() - 1;
1164 auto host_orig_row_ptr = Kokkos::create_mirror_view(orig_row_ptr);
1165 auto new_nnz = host_orig_row_ptr(size);
1167 values_view_t new_values(
1168 Kokkos::ViewAllocateWithoutInitializing(
"new_values"), values.size() - new_nnz/2);
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));
1183 new_values(tk) = values(sk);
1189 values = new_values;
1192 template<
class array_view_t,
class per_view_t>
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));
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);
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
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 "-"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