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);
616 return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
620 int rank = Teuchos::rank(*comm);
621 size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
622 if( rank == 0 ) my_num_elems = num_global_elements;
623 return rcp(
new Tpetra::Map<LO,GO, Node>(num_global_elements,
624 my_num_elems, indexBase, comm));
628 TEUCHOS_TEST_FOR_EXCEPTION(
true,
630 "Control should never reach this point. " 631 "Please contact the Amesos2 developers." );
637 #ifdef HAVE_TPETRA_INST_INT_INT 638 #ifdef HAVE_AMESOS2_EPETRA 643 template <
typename LO,
typename GO,
typename GS,
typename Node>
644 Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
645 epetra_map_to_tpetra_map(
const Epetra_BlockMap& map)
650 int num_my_elements = map.NumMyElements();
651 Teuchos::Array<int> my_global_elements(num_my_elements);
652 map.MyGlobalElements(my_global_elements.getRawPtr());
654 typedef Tpetra::Map<LO,GO,Node> map_t;
655 RCP<map_t> tmap = rcp(
new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
656 my_global_elements(),
657 as<GO>(map.IndexBase()),
658 to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
662 template <
typename LO,
typename GO,
typename GS,
typename Node>
663 Teuchos::RCP<Epetra_Map>
664 tpetra_map_to_epetra_map(
const Tpetra::Map<LO,GO,Node>& map)
668 Teuchos::Array<GO> elements_tmp;
669 elements_tmp = map.getNodeElementList();
670 int num_my_elements = elements_tmp.size();
671 Teuchos::Array<int> my_global_elements(num_my_elements);
672 for (
int i = 0; i < num_my_elements; ++i){
673 my_global_elements[i] = as<int>(elements_tmp[i]);
677 RCP<Epetra_Map> emap = rcp(
new Epetra_Map(-1,
679 my_global_elements.getRawPtr(),
680 as<GO>(map.getIndexBase()),
681 *to_epetra_comm(map.getComm())));
684 #endif // HAVE_AMESOS2_EPETRA 685 #endif // HAVE_TPETRA_INST_INT_INT 687 template <
typename Scalar,
688 typename GlobalOrdinal,
689 typename GlobalSizeT>
690 void transpose(Teuchos::ArrayView<Scalar> vals,
691 Teuchos::ArrayView<GlobalOrdinal> indices,
692 Teuchos::ArrayView<GlobalSizeT> ptr,
693 Teuchos::ArrayView<Scalar> trans_vals,
694 Teuchos::ArrayView<GlobalOrdinal> trans_indices,
695 Teuchos::ArrayView<GlobalSizeT> trans_ptr)
703 #ifdef HAVE_AMESOS2_DEBUG 704 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
705 ind_begin = indices.begin();
706 ind_end = indices.end();
707 size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
708 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
709 std::invalid_argument,
710 "Transpose pointer size not large enough." );
711 TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
712 std::invalid_argument,
713 "Transpose values array not large enough." );
714 TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
715 std::invalid_argument,
716 "Transpose indices array not large enough." );
718 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
722 Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
723 ind_end = indices.end();
724 for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
725 ++(count[(*ind_it) + 1]);
729 typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
730 cnt_end = count.end();
731 for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
732 *cnt_it = *cnt_it + *(cnt_it - 1);
735 trans_ptr.assign(count);
749 GlobalSizeT size = ptr.size();
750 for( GlobalSizeT i = 0; i < size - 1; ++i ){
751 GlobalOrdinal u = ptr[i];
752 GlobalOrdinal v = ptr[i + 1];
753 for( GlobalOrdinal j = u; j < v; ++j ){
754 GlobalOrdinal k = count[indices[j]];
755 trans_vals[k] = vals[j];
756 trans_indices[k] = i;
757 ++(count[indices[j]]);
763 template <
typename Scalar1,
typename Scalar2>
765 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
766 size_t ld, Teuchos::ArrayView<Scalar2> s)
768 size_t vals_size = vals.size();
769 #ifdef HAVE_AMESOS2_DEBUG 770 size_t s_size = s.size();
771 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
772 std::invalid_argument,
773 "Scale vector must have length at least that of the vector" );
776 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
786 template <
typename Scalar1,
typename Scalar2,
class BinaryOp>
788 scale(Teuchos::ArrayView<Scalar1> vals,
size_t l,
789 size_t ld, Teuchos::ArrayView<Scalar2> s,
792 size_t vals_size = vals.size();
793 #ifdef HAVE_AMESOS2_DEBUG 794 size_t s_size = s.size();
795 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
796 std::invalid_argument,
797 "Scale vector must have length at least that of the vector" );
800 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
806 vals[i] = binary_op(vals[i], s[s_i]);
816 #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