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> 67 #ifdef HAVE_TPETRA_INST_INT_INT 68 #ifdef HAVE_AMESOS2_EPETRA 69 #include <Epetra_Map.h> 85 using Teuchos::ArrayView;
88 using Meta::if_then_else;
105 template <
typename LO,
typename GO,
typename GS,
typename Node>
106 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
108 GS num_global_elements,
109 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
113 #ifdef HAVE_TPETRA_INST_INT_INT 114 #ifdef HAVE_AMESOS2_EPETRA 121 template <
typename LO,
typename GO,
typename GS,
typename Node>
122 RCP<Tpetra::Map<LO,GO,Node> >
123 epetra_map_to_tpetra_map(
const Epetra_BlockMap& map);
130 template <
typename LO,
typename GO,
typename GS,
typename Node>
132 tpetra_map_to_epetra_map(
const Tpetra::Map<LO,GO,Node>& map);
139 const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
146 const RCP<const Epetra_Comm> to_epetra_comm(RCP<
const Teuchos::Comm<int> > c);
147 #endif // HAVE_AMESOS2_EPETRA 148 #endif // HAVE_TPETRA_INST_INT_INT 155 template <
typename Scalar,
156 typename GlobalOrdinal,
157 typename GlobalSizeT>
159 ArrayView<GlobalOrdinal> indices,
160 ArrayView<GlobalSizeT> ptr,
161 ArrayView<Scalar> trans_vals,
162 ArrayView<GlobalOrdinal> trans_indices,
163 ArrayView<GlobalSizeT> trans_ptr);
178 template <
typename Scalar1,
typename Scalar2>
179 void scale(ArrayView<Scalar1> vals,
size_t l,
180 size_t ld, ArrayView<Scalar2> s);
200 template <
typename Scalar1,
typename Scalar2,
class BinaryOp>
201 void scale(ArrayView<Scalar1> vals,
size_t l,
202 size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
206 void printLine( Teuchos::FancyOStream &out );
213 #ifndef DOXYGEN_SHOULD_SKIP_THIS 236 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
237 struct same_gs_helper
239 static void do_get(
const Teuchos::Ptr<const M> mat,
240 const ArrayView<typename M::scalar_t> nzvals,
241 const ArrayView<typename M::global_ordinal_t> indices,
242 const ArrayView<GS> pointers,
245 const Tpetra::Map<
typename M::local_ordinal_t,
246 typename M::global_ordinal_t,
247 typename M::node_t> > map,
250 Op::apply(mat, nzvals, indices, pointers, nnz, map, ordering);
254 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
255 struct diff_gs_helper
257 static void do_get(
const Teuchos::Ptr<const M> mat,
258 const ArrayView<typename M::scalar_t> nzvals,
259 const ArrayView<typename M::global_ordinal_t> indices,
260 const ArrayView<GS> pointers,
263 const Tpetra::Map<
typename M::local_ordinal_t,
264 typename M::global_ordinal_t,
265 typename M::node_t> > map,
268 typedef typename M::global_size_t mat_gs_t;
269 typename ArrayView<GS>::size_type i, size = pointers.size();
270 Teuchos::Array<mat_gs_t> pointers_tmp(size);
271 mat_gs_t nnz_tmp = 0;
273 Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, ordering);
275 for (i = 0; i < size; ++i){
276 pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
278 nnz = Teuchos::as<GS>(nnz_tmp);
282 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
283 struct same_go_helper
285 static void do_get(
const Teuchos::Ptr<const M> mat,
286 const ArrayView<typename M::scalar_t> nzvals,
287 const ArrayView<GO> indices,
288 const ArrayView<GS> pointers,
291 const Tpetra::Map<
typename M::local_ordinal_t,
292 typename M::global_ordinal_t,
293 typename M::node_t> > map,
296 typedef typename M::global_size_t mat_gs_t;
297 if_then_else<is_same<GS,mat_gs_t>::value,
298 same_gs_helper<M,S,GO,GS,Op>,
299 diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
305 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
306 struct diff_go_helper
308 static void do_get(
const Teuchos::Ptr<const M> mat,
309 const ArrayView<typename M::scalar_t> nzvals,
310 const ArrayView<GO> indices,
311 const ArrayView<GS> pointers,
314 const Tpetra::Map<
typename M::local_ordinal_t,
315 typename M::global_ordinal_t,
316 typename M::node_t> > map,
319 typedef typename M::global_ordinal_t mat_go_t;
320 typedef typename M::global_size_t mat_gs_t;
321 typename ArrayView<GO>::size_type i, size = indices.size();
322 Teuchos::Array<mat_go_t> indices_tmp(size);
324 if_then_else<is_same<GS,mat_gs_t>::value,
325 same_gs_helper<M,S,GO,GS,Op>,
326 diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
330 for (i = 0; i < size; ++i){
331 indices[i] = Teuchos::as<GO>(indices_tmp[i]);
336 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
337 struct same_scalar_helper
339 static void do_get(
const Teuchos::Ptr<const M> mat,
340 const ArrayView<S> nzvals,
341 const ArrayView<GO> indices,
342 const ArrayView<GS> pointers,
345 const Tpetra::Map<
typename M::local_ordinal_t,
346 typename M::global_ordinal_t,
347 typename M::node_t> > map,
350 typedef typename M::global_ordinal_t mat_go_t;
351 if_then_else<is_same<GO,mat_go_t>::value,
352 same_go_helper<M,S,GO,GS,Op>,
353 diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
359 template <
class M,
typename S,
typename GO,
typename GS,
class Op>
360 struct diff_scalar_helper
362 static void do_get(
const Teuchos::Ptr<const M> mat,
363 const ArrayView<S> nzvals,
364 const ArrayView<GO> indices,
365 const ArrayView<GS> pointers,
368 const Tpetra::Map<
typename M::local_ordinal_t,
369 typename M::global_ordinal_t,
370 typename M::node_t> > map,
373 typedef typename M::scalar_t mat_scalar_t;
374 typedef typename M::global_ordinal_t mat_go_t;
375 typename ArrayView<S>::size_type i, size = nzvals.size();
376 Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
378 if_then_else<is_same<GO,mat_go_t>::value,
379 same_go_helper<M,S,GO,GS,Op>,
380 diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
384 for (i = 0; i < size; ++i){
385 nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
389 #endif // DOXYGEN_SHOULD_SKIP_THIS 403 template<
class Matrix,
typename S,
typename GO,
typename GS,
class Op>
406 static void do_get(
const Teuchos::Ptr<const Matrix> mat,
407 const ArrayView<S> nzvals,
408 const ArrayView<GO> indices,
409 const ArrayView<GS> pointers,
414 typedef typename Matrix::local_ordinal_t lo_t;
415 typedef typename Matrix::global_ordinal_t go_t;
416 typedef typename Matrix::global_size_t gs_t;
417 typedef typename Matrix::node_t node_t;
419 const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
420 = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
421 Op::get_dimension(mat),
422 mat->getComm(), indexBase);
423 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
430 static void do_get(
const Teuchos::Ptr<const Matrix> mat,
431 const ArrayView<S> nzvals,
432 const ArrayView<GO> indices,
433 const ArrayView<GS> pointers,
436 const Teuchos::RCP<
const Tpetra::Map<
typename Matrix::local_ordinal_t,
437 typename Matrix::global_ordinal_t,
438 typename Matrix::node_t> > map
440 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
447 static void do_get(
const Teuchos::Ptr<const Matrix> mat,
448 const ArrayView<S> nzvals,
449 const ArrayView<GO> indices,
450 const ArrayView<GS> pointers,
453 const Tpetra::Map<
typename Matrix::local_ordinal_t,
454 typename Matrix::global_ordinal_t,
455 typename Matrix::node_t> > map,
458 typedef typename Matrix::scalar_t mat_scalar;
459 if_then_else<is_same<mat_scalar,S>::value,
460 same_scalar_helper<Matrix,S,GO,GS,Op>,
461 diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
468 #ifndef DOXYGEN_SHOULD_SKIP_THIS 473 template<
class Matrix>
476 static void apply(
const Teuchos::Ptr<const Matrix> mat,
477 const ArrayView<typename Matrix::scalar_t> nzvals,
478 const ArrayView<typename Matrix::global_ordinal_t> rowind,
479 const ArrayView<typename Matrix::global_size_t> colptr,
480 typename Matrix::global_size_t& nnz,
482 const Tpetra::Map<
typename Matrix::local_ordinal_t,
483 typename Matrix::global_ordinal_t,
484 typename Matrix::node_t> > map,
487 mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
491 const Teuchos::RCP<
const Tpetra::Map<
typename Matrix::local_ordinal_t,
492 typename Matrix::global_ordinal_t,
493 typename Matrix::node_t> >
494 getMap(
const Teuchos::Ptr<const Matrix> mat)
496 return mat->getColMap();
500 typename Matrix::global_size_t
501 get_dimension(
const Teuchos::Ptr<const Matrix> mat)
503 return mat->getGlobalNumCols();
507 template<
class Matrix>
510 static void apply(
const Teuchos::Ptr<const Matrix> mat,
511 const ArrayView<typename Matrix::scalar_t> nzvals,
512 const ArrayView<typename Matrix::global_ordinal_t> colind,
513 const ArrayView<typename Matrix::global_size_t> rowptr,
514 typename Matrix::global_size_t& nnz,
516 const Tpetra::Map<
typename Matrix::local_ordinal_t,
517 typename Matrix::global_ordinal_t,
518 typename Matrix::node_t> > map,
521 mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering);
525 const Teuchos::RCP<
const Tpetra::Map<
typename Matrix::local_ordinal_t,
526 typename Matrix::global_ordinal_t,
527 typename Matrix::node_t> >
528 getMap(
const Teuchos::Ptr<const Matrix> mat)
530 return mat->getRowMap();
534 typename Matrix::global_size_t
535 get_dimension(
const Teuchos::Ptr<const Matrix> mat)
537 return mat->getGlobalNumRows();
540 #endif // DOXYGEN_SHOULD_SKIP_THIS 579 template<
class Matrix,
typename S,
typename GO,
typename GS>
590 template<
class Matrix,
typename S,
typename GO,
typename GS>
601 template <
typename LO,
typename GO,
typename GS,
typename Node>
602 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
604 GS num_global_elements,
605 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
610 switch( distribution ){
613 return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
615 return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
618 int rank = Teuchos::rank(*comm);
619 size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
620 if( rank == 0 ) my_num_elems = num_global_elements;
621 return rcp(
new Tpetra::Map<LO,GO, Node>(num_global_elements,
622 my_num_elems, indexBase, comm));
625 TEUCHOS_TEST_FOR_EXCEPTION(
true,
627 "Control should never reach this point. " 628 "Please contact the Amesos2 developers." );
633 #ifdef HAVE_TPETRA_INST_INT_INT 634 #ifdef HAVE_AMESOS2_EPETRA 639 template <
typename LO,
typename GO,
typename GS,
typename Node>
640 Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
641 epetra_map_to_tpetra_map(
const Epetra_BlockMap& map)
646 int num_my_elements = map.NumMyElements();
647 Teuchos::Array<int> my_global_elements(num_my_elements);
648 map.MyGlobalElements(my_global_elements.getRawPtr());
650 typedef Tpetra::Map<LO,GO,Node> map_t;
651 RCP<map_t> tmap = rcp(
new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
652 my_global_elements(),
653 as<GO>(map.IndexBase()),
654 to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
658 template <
typename LO,
typename GO,
typename GS,
typename Node>
659 Teuchos::RCP<Epetra_Map>
660 tpetra_map_to_epetra_map(
const Tpetra::Map<LO,GO,Node>& map)
664 Teuchos::Array<GO> elements_tmp;
665 elements_tmp = map.getNodeElementList();
666 int num_my_elements = elements_tmp.size();
667 Teuchos::Array<int> my_global_elements(num_my_elements);
668 for (
int i = 0; i < num_my_elements; ++i){
669 my_global_elements[i] = as<int>(elements_tmp[i]);
673 RCP<Epetra_Map> emap = rcp(
new Epetra_Map(-1,
675 my_global_elements.getRawPtr(),
676 as<GO>(map.getIndexBase()),
677 *to_epetra_comm(map.getComm())));
680 #endif // HAVE_AMESOS2_EPETRA 681 #endif // HAVE_TPETRA_INST_INT_INT 683 template <
typename Scalar,
684 typename GlobalOrdinal,
685 typename GlobalSizeT>
686 void transpose(Teuchos::ArrayView<Scalar> vals,
687 Teuchos::ArrayView<GlobalOrdinal> indices,
688 Teuchos::ArrayView<GlobalSizeT> ptr,
689 Teuchos::ArrayView<Scalar> trans_vals,
690 Teuchos::ArrayView<GlobalOrdinal> trans_indices,
691 Teuchos::ArrayView<GlobalSizeT> trans_ptr)
699 #ifdef HAVE_AMESOS2_DEBUG 700 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
701 ind_begin = indices.begin();
702 ind_end = indices.end();
703 size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
704 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
705 std::invalid_argument,
706 "Transpose pointer size not large enough." );
707 TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
708 std::invalid_argument,
709 "Transpose values array not large enough." );
710 TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
711 std::invalid_argument,
712 "Transpose indices array not large enough." );
714 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
718 Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
719 ind_end = indices.end();
720 for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
721 ++(count[(*ind_it) + 1]);
725 typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
726 cnt_end = count.end();
727 for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
728 *cnt_it = *cnt_it + *(cnt_it - 1);
731 trans_ptr.assign(count);
745 GlobalSizeT size = ptr.size();
746 for( GlobalSizeT i = 0; i < size - 1; ++i ){
747 GlobalOrdinal u = ptr[i];
748 GlobalOrdinal v = ptr[i + 1];
749 for( GlobalOrdinal j = u; j < v; ++j ){
750 GlobalOrdinal k = count[indices[j]];
751 trans_vals[k] = vals[j];
752 trans_indices[k] = i;
753 ++(count[indices[j]]);
759 template <
typename Scalar1,
typename Scalar2>
761 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
762 size_t ld, Teuchos::ArrayView<Scalar2> s)
764 size_t vals_size = vals.size();
765 #ifdef HAVE_AMESOS2_DEBUG 766 size_t s_size = s.size();
767 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
768 std::invalid_argument,
769 "Scale vector must have length at least that of the vector" );
772 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
782 template <
typename Scalar1,
typename Scalar2,
class BinaryOp>
784 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
785 size_t ld, Teuchos::ArrayView<Scalar2> s,
788 size_t vals_size = vals.size();
789 #ifdef HAVE_AMESOS2_DEBUG 790 size_t s_size = s.size();
791 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
792 std::invalid_argument,
793 "Scale vector must have length at least that of the vector" );
796 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
802 vals[i] = binary_op(vals[i], s[s_i]);
812 #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.
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, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:447
const Teuchos::RCP< const Tpetra::Map< LO, GO, Node > > getDistributionMap(EDistribution distribution, GS num_global_elements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, GO indexBase=0)
Gets a Tpetra::Map described by the EDistribution.
Definition: Amesos2_Util.hpp:603
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:591
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:404
Definition: Amesos2_TypeDecl.hpp:142
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:140
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)
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:580
Definition: Amesos2_TypeDecl.hpp:126
static void do_get(const Teuchos::Ptr< const Matrix > mat, const ArrayView< S > nzvals, const ArrayView< GO > indices, const ArrayView< GS > pointers, GS &nnz, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:430
Definition: Amesos2_TypeDecl.hpp:124
Enum and other types declarations for Amesos2.
Definition: Amesos2_TypeDecl.hpp:127
EDistribution
Definition: Amesos2_TypeDecl.hpp:123