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"
70
71#ifdef HAVE_AMESOS2_EPETRA
72#include <Epetra_Map.h>
73#endif
74
75#ifdef HAVE_AMESOS2_METIS
76#include "metis.h" // to discuss, remove from header?
77#endif
78
79namespace Amesos2 {
80
81 namespace Util {
82
89 using Teuchos::RCP;
90 using Teuchos::ArrayView;
91
92 using Meta::is_same;
93 using Meta::if_then_else;
94
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 );
114
115
116 template <typename LO, typename GO, typename GS, typename Node>
117 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
118 getDistributionMap(EDistribution distribution,
119 GS num_global_elements,
120 const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
121 GO indexBase = 0,
122 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map = Teuchos::null);
123
124
125#ifdef HAVE_AMESOS2_EPETRA
126
132 template <typename LO, typename GO, typename GS, typename Node>
133 RCP<Tpetra::Map<LO,GO,Node> >
134 epetra_map_to_tpetra_map(const Epetra_BlockMap& map);
135
141 template <typename LO, typename GO, typename GS, typename Node>
142 RCP<Epetra_Map>
143 tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map);
144
150 const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
151
157 const RCP<const Epetra_Comm> to_epetra_comm(RCP<const Teuchos::Comm<int> > c);
158#endif // HAVE_AMESOS2_EPETRA
159
165 template <typename Scalar,
166 typename GlobalOrdinal,
167 typename GlobalSizeT>
168 void transpose(ArrayView<Scalar> vals,
169 ArrayView<GlobalOrdinal> indices,
170 ArrayView<GlobalSizeT> ptr,
171 ArrayView<Scalar> trans_vals,
172 ArrayView<GlobalOrdinal> trans_indices,
173 ArrayView<GlobalSizeT> trans_ptr);
174
188 template <typename Scalar1, typename Scalar2>
189 void scale(ArrayView<Scalar1> vals, size_t l,
190 size_t ld, ArrayView<Scalar2> s);
191
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);
213
214
216 void printLine( Teuchos::FancyOStream &out );
217
218 // Helper function used to convert Kokkos::complex pointer
219 // to std::complex pointer; needed for optimized code path
220 // when retrieving the CRS raw pointers
221 template < class T0, class T1 >
222 struct getStdCplxType
223 {
224 using common_type = typename std::common_type<T0,T1>::type;
225 using type = common_type;
226 };
227
228 template < class T0, class T1 >
229 struct getStdCplxType< T0, T1* >
230 {
231 using common_type = typename std::common_type<T0,T1>::type;
232 using type = common_type;
233 };
234
235#if defined(HAVE_TEUCHOS_COMPLEX) && defined(HAVE_AMESOS2_KOKKOS)
236 template < class T0 >
237 struct getStdCplxType< T0, Kokkos::complex<T0>* >
238 {
239 using type = std::complex<T0>;
240 };
241
242 template < class T0 , class T1 >
243 struct getStdCplxType< T0, Kokkos::complex<T1>* >
244 {
245 using common_type = typename std::common_type<T0,T1>::type;
246 using type = std::complex<common_type>;
247 };
248#endif
249
251 // Matrix/MultiVector Utilities //
253
254
255
269 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
271 {
272 static void do_get(const Teuchos::Ptr<const M> mat,
273 KV_S& nzvals,
274 KV_GO& indices,
275 KV_GS& pointers,
276 typename KV_GS::value_type& nnz,
277 const Teuchos::Ptr<
278 const Tpetra::Map<typename M::local_ordinal_t,
279 typename M::global_ordinal_t,
280 typename M::node_t> > map,
281 EDistribution distribution,
282 EStorage_Ordering ordering)
283 {
284 Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
285 indices, pointers, nnz, map, distribution, ordering);
286 }
287 };
288
289 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
290 struct diff_gs_helper_kokkos_view
291 {
292 static void do_get(const Teuchos::Ptr<const M> mat,
293 KV_S& nzvals,
294 KV_GO& indices,
295 KV_GS& pointers,
296 typename KV_GS::value_type& nnz,
297 const Teuchos::Ptr<
298 const Tpetra::Map<typename M::local_ordinal_t,
299 typename M::global_ordinal_t,
300 typename M::node_t> > map,
301 EDistribution distribution,
302 EStorage_Ordering ordering)
303 {
304 typedef typename M::global_size_t mat_gs_t;
305 typedef typename Kokkos::View<mat_gs_t*, Kokkos::HostSpace> KV_TMP;
306 size_t i, size = pointers.extent(0);
307 KV_TMP pointers_tmp(Kokkos::ViewAllocateWithoutInitializing("pointers_tmp"), size);
308
309 mat_gs_t nnz_tmp = 0;
310 Op::template apply_kokkos_view<KV_S, KV_GO, KV_TMP>(mat, nzvals,
311 indices, pointers_tmp, nnz_tmp, Teuchos::ptrInArg(*map), distribution, ordering);
312 nnz = Teuchos::as<typename KV_GS::value_type>(nnz_tmp);
313
314 typedef typename KV_GS::value_type view_gs_t;
315 for (i = 0; i < size; ++i){
316 pointers(i) = Teuchos::as<view_gs_t>(pointers_tmp(i));
317 }
318 nnz = Teuchos::as<view_gs_t>(nnz_tmp);
319 }
320 };
321
322 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
323 struct same_go_helper_kokkos_view
324 {
325 static void do_get(const Teuchos::Ptr<const M> mat,
326 KV_S& nzvals,
327 KV_GO& indices,
328 KV_GS& pointers,
329 typename KV_GS::value_type& nnz,
330 const Teuchos::Ptr<
331 const Tpetra::Map<typename M::local_ordinal_t,
332 typename M::global_ordinal_t,
333 typename M::node_t> > map,
334 EDistribution distribution,
335 EStorage_Ordering ordering)
336 {
337 typedef typename M::global_size_t mat_gs_t;
338 typedef typename KV_GS::value_type view_gs_t;
339 if_then_else<is_same<view_gs_t,mat_gs_t>::value,
340 same_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
341 diff_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
342 pointers, nnz, map,
343 distribution, ordering);
344 }
345 };
346
347 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
348 struct diff_go_helper_kokkos_view
349 {
350 static void do_get(const Teuchos::Ptr<const M> mat,
351 KV_S& nzvals,
352 KV_GO& indices,
353 KV_GS& pointers,
354 typename KV_GS::value_type& nnz,
355 const Teuchos::Ptr<
356 const Tpetra::Map<typename M::local_ordinal_t,
357 typename M::global_ordinal_t,
358 typename M::node_t> > map,
359 EDistribution distribution,
360 EStorage_Ordering ordering)
361 {
362 typedef typename M::global_ordinal_t mat_go_t;
363 typedef typename M::global_size_t mat_gs_t;
364 typedef typename Kokkos::View<mat_go_t*, Kokkos::HostSpace> KV_TMP;
365 size_t i, size = indices.extent(0);
366 KV_TMP indices_tmp(Kokkos::ViewAllocateWithoutInitializing("indices_tmp"), size);
367
368 typedef typename KV_GO::value_type view_go_t;
369 typedef typename KV_GS::value_type view_gs_t;
370 if_then_else<is_same<view_gs_t,mat_gs_t>::value,
371 same_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op>,
372 diff_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op> >::type::do_get(mat, nzvals, indices_tmp,
373 pointers, nnz, map,
374 distribution, ordering);
375 for (i = 0; i < size; ++i){
376 indices(i) = Teuchos::as<view_go_t>(indices_tmp(i));
377 }
378 }
379 };
380
381 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
382 struct same_scalar_helper_kokkos_view
383 {
384 static void do_get(const Teuchos::Ptr<const M> mat,
385 KV_S& nzvals,
386 KV_GO& indices,
387 KV_GS& pointers,
388 typename KV_GS::value_type& nnz,
389 const Teuchos::Ptr<
390 const Tpetra::Map<typename M::local_ordinal_t,
391 typename M::global_ordinal_t,
392 typename M::node_t> > map,
393 EDistribution distribution,
394 EStorage_Ordering ordering)
395 {
396 typedef typename M::global_ordinal_t mat_go_t;
397 typedef typename KV_GO::value_type view_go_t;
398 if_then_else<is_same<view_go_t,mat_go_t>::value,
399 same_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
400 diff_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
401 pointers, nnz, map,
402 distribution, ordering);
403 }
404 };
405
406 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
407 struct diff_scalar_helper_kokkos_view
408 {
409 static void do_get(const Teuchos::Ptr<const M> mat,
410 KV_S& nzvals,
411 KV_GO& indices,
412 KV_GS& pointers,
413 typename KV_GS::value_type& nnz,
414 const Teuchos::Ptr<
415 const Tpetra::Map<typename M::local_ordinal_t,
416 typename M::global_ordinal_t,
417 typename M::node_t> > map,
418 EDistribution distribution,
419 EStorage_Ordering ordering)
420 {
421 typedef typename M::global_ordinal_t mat_go_t;
422 typedef typename Kokkos::ArithTraits<typename M::scalar_t>::val_type mat_scalar_t;
423 typedef typename Kokkos::View<mat_scalar_t*, Kokkos::HostSpace> KV_TMP;
424 size_t i, size = nzvals.extent(0);
425 KV_TMP nzvals_tmp(Kokkos::ViewAllocateWithoutInitializing("nzvals_tmp"), size);
426
427 typedef typename KV_S::value_type view_scalar_t;
428 typedef typename KV_GO::value_type view_go_t;
429 if_then_else<is_same<view_go_t,mat_go_t>::value,
430 same_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op>,
431 diff_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals_tmp, indices,
432 pointers, nnz, map,
433 distribution, ordering);
434
435 for (i = 0; i < size; ++i){
436 nzvals(i) = Teuchos::as<view_scalar_t>(nzvals_tmp(i));
437 }
438 }
439 };
440
441
442 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS, class Op>
443 struct get_cxs_helper_kokkos_view
444 {
445 static void do_get(const Teuchos::Ptr<const Matrix> mat,
446 KV_S& nzvals,
447 KV_GO& indices,
448 KV_GS& pointers,
449 typename KV_GS::value_type& nnz,
450 EDistribution distribution,
451 EStorage_Ordering ordering=ARBITRARY,
452 typename KV_GO::value_type indexBase = 0)
453 {
454 typedef typename Matrix::local_ordinal_t lo_t;
455 typedef typename Matrix::global_ordinal_t go_t;
456 typedef typename Matrix::global_size_t gs_t;
457 typedef typename Matrix::node_t node_t;
458
459 const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
460 = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
461 Op::get_dimension(mat),
462 mat->getComm(),
463 indexBase,
464 Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
465 );
466 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
467 }
468
473 static void do_get(const Teuchos::Ptr<const Matrix> mat,
474 KV_S& nzvals,
475 KV_GO& indices,
476 KV_GS& pointers,
477 typename KV_GS::value_type& nnz,
478 EDistribution distribution, // Does this one need a distribution argument??
479 EStorage_Ordering ordering=ARBITRARY)
480 {
481 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
482 typename Matrix::global_ordinal_t,
483 typename Matrix::node_t> > map
484 = Op::getMap(mat);
485 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
486 }
487
492 static void do_get(const Teuchos::Ptr<const Matrix> mat,
493 KV_S& nzvals,
494 KV_GO& indices,
495 KV_GS& pointers,
496 typename KV_GS::value_type& nnz,
497 const Teuchos::Ptr<
498 const Tpetra::Map<typename Matrix::local_ordinal_t,
499 typename Matrix::global_ordinal_t,
500 typename Matrix::node_t> > map,
501 EDistribution distribution,
502 EStorage_Ordering ordering=ARBITRARY)
503 {
504 typedef typename Matrix::scalar_t mat_scalar;
505 typedef typename KV_S::value_type view_scalar_t;
506
507 if_then_else<is_same<mat_scalar,view_scalar_t>::value,
508 same_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op>,
509 diff_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op> >::type::do_get(mat,
510 nzvals, indices,
511 pointers, nnz,
512 map,
513 distribution, ordering);
514 }
515 };
516
517#ifndef DOXYGEN_SHOULD_SKIP_THIS
518 /*
519 * These two function-like classes are meant to be used as the \c
520 * Op template parameter for the \c get_cxs_helper template class.
521 */
522 template<class Matrix>
523 struct get_ccs_func
524 {
525 template<typename KV_S, typename KV_GO, typename KV_GS>
526 static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
527 KV_S& nzvals,
528 KV_GO& rowind,
529 KV_GS& colptr,
530 typename Matrix::global_size_t& nnz,
531 const Teuchos::Ptr<
532 const Tpetra::Map<typename Matrix::local_ordinal_t,
533 typename Matrix::global_ordinal_t,
534 typename Matrix::node_t> > map,
535 EDistribution distribution,
536 EStorage_Ordering ordering)
537 {
538 mat->getCcs_kokkos_view(nzvals, rowind, colptr, nnz, map, ordering, distribution);
539 }
540
541 static
542 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
543 typename Matrix::global_ordinal_t,
544 typename Matrix::node_t> >
545 getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
546 {
547 return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
548 }
549
550 static
551 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
552 typename Matrix::global_ordinal_t,
553 typename Matrix::node_t> >
554 getMap(const Teuchos::Ptr<const Matrix> mat)
555 {
556 return mat->getColMap();
557 }
558
559 static
560 typename Matrix::global_size_t
561 get_dimension(const Teuchos::Ptr<const Matrix> mat)
562 {
563 return mat->getGlobalNumCols();
564 }
565 };
566
567 template<class Matrix>
568 struct get_crs_func
569 {
570 template<typename KV_S, typename KV_GO, typename KV_GS>
571 static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
572 KV_S& nzvals,
573 KV_GO& colind,
574 KV_GS& rowptr,
575 typename Matrix::global_size_t& nnz,
576 const Teuchos::Ptr<
577 const Tpetra::Map<typename Matrix::local_ordinal_t,
578 typename Matrix::global_ordinal_t,
579 typename Matrix::node_t> > map,
580 EDistribution distribution,
581 EStorage_Ordering ordering)
582 {
583 mat->getCrs_kokkos_view(nzvals, colind, rowptr, nnz, map, ordering, distribution);
584 }
585
586 static
587 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
588 typename Matrix::global_ordinal_t,
589 typename Matrix::node_t> >
590 getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
591 {
592 return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
593 }
594
595 static
596 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
597 typename Matrix::global_ordinal_t,
598 typename Matrix::node_t> >
599 getMap(const Teuchos::Ptr<const Matrix> mat)
600 {
601 return mat->getRowMap();
602 }
603
604 static
605 typename Matrix::global_size_t
606 get_dimension(const Teuchos::Ptr<const Matrix> mat)
607 {
608 return mat->getGlobalNumRows();
609 }
610 };
611#endif // DOXYGEN_SHOULD_SKIP_THIS
612
650 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
651 struct get_ccs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_ccs_func<Matrix> >
652 {};
653
661 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
662 struct get_crs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_crs_func<Matrix> >
663 {};
664 /* End Matrix/MultiVector Utilities */
665
666
668 // Definitions //
670
671
672 template <typename LO, typename GO, typename GS, typename Node>
673 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
674 getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
675 {
676 //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
677 Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
678 return gather_map;
679 }
680
681
682 template <typename LO, typename GO, typename GS, typename Node>
683 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
684 getDistributionMap(EDistribution distribution,
685 GS num_global_elements,
686 const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
687 GO indexBase,
688 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
689 {
690 // TODO: Need to add indexBase to cases other than ROOTED
691 // We do not support these maps in any solver now.
692 switch( distribution ){
693 case DISTRIBUTED:
695 return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
697 return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
698 case ROOTED:
699 {
700 int rank = Teuchos::rank(*comm);
701 size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
702 if( rank == 0 ) { my_num_elems = num_global_elements; }
703
704 return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
705 my_num_elems, indexBase, comm));
706 }
708 {
709 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
710 = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
711 return gathermap;
712 }
713 default:
714 TEUCHOS_TEST_FOR_EXCEPTION( true,
715 std::logic_error,
716 "Control should never reach this point. "
717 "Please contact the Amesos2 developers." );
718 }
719 }
720
721
722#ifdef HAVE_AMESOS2_EPETRA
723
724 //#pragma message "include 3"
725 //#include <Epetra_Map.h>
726
727 template <typename LO, typename GO, typename GS, typename Node>
728 Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
729 epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
730 {
731 using Teuchos::as;
732 using Teuchos::rcp;
733
734 int num_my_elements = map.NumMyElements();
735 Teuchos::Array<int> my_global_elements(num_my_elements);
736 map.MyGlobalElements(my_global_elements.getRawPtr());
737
738 Teuchos::Array<GO> my_gbl_inds_buf;
739 Teuchos::ArrayView<GO> my_gbl_inds;
740 if (! std::is_same<int, GO>::value) {
741 my_gbl_inds_buf.resize (num_my_elements);
742 my_gbl_inds = my_gbl_inds_buf ();
743 for (int k = 0; k < num_my_elements; ++k) {
744 my_gbl_inds[k] = static_cast<GO> (my_global_elements[k]);
745 }
746 }
747 else {
748 using Teuchos::av_reinterpret_cast;
749 my_gbl_inds = av_reinterpret_cast<GO> (my_global_elements ());
750 }
751
752 typedef Tpetra::Map<LO,GO,Node> map_t;
753 RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
754 my_gbl_inds(),
755 as<GO>(map.IndexBase()),
756 to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
757 return tmap;
758 }
759
760 template <typename LO, typename GO, typename GS, typename Node>
761 Teuchos::RCP<Epetra_Map>
762 tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
763 {
764 using Teuchos::as;
765
766 Teuchos::Array<GO> elements_tmp;
767 elements_tmp = map.getLocalElementList();
768 int num_my_elements = elements_tmp.size();
769 Teuchos::Array<int> my_global_elements(num_my_elements);
770 for (int i = 0; i < num_my_elements; ++i){
771 my_global_elements[i] = as<int>(elements_tmp[i]);
772 }
773
774 using Teuchos::rcp;
775 RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
776 num_my_elements,
777 my_global_elements.getRawPtr(),
778 as<GO>(map.getIndexBase()),
779 *to_epetra_comm(map.getComm())));
780 return emap;
781 }
782#endif // HAVE_AMESOS2_EPETRA
783
784 template <typename Scalar,
785 typename GlobalOrdinal,
786 typename GlobalSizeT>
787 void transpose(Teuchos::ArrayView<Scalar> vals,
788 Teuchos::ArrayView<GlobalOrdinal> indices,
789 Teuchos::ArrayView<GlobalSizeT> ptr,
790 Teuchos::ArrayView<Scalar> trans_vals,
791 Teuchos::ArrayView<GlobalOrdinal> trans_indices,
792 Teuchos::ArrayView<GlobalSizeT> trans_ptr)
793 {
794 /* We have a compressed-row storage format of this matrix. We
795 * transform this into a compressed-column format using a
796 * distribution-counting sort algorithm, which is described by
797 * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
798 */
799
800#ifdef HAVE_AMESOS2_DEBUG
801 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
802 ind_begin = indices.begin();
803 ind_end = indices.end();
804 size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
805 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
806 std::invalid_argument,
807 "Transpose pointer size not large enough." );
808 TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
809 std::invalid_argument,
810 "Transpose values array not large enough." );
811 TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
812 std::invalid_argument,
813 "Transpose indices array not large enough." );
814#else
815 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
816#endif
817 // Count the number of entries in each column
818 Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
819 ind_end = indices.end();
820 for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
821 ++(count[(*ind_it) + 1]);
822 }
823 // Accumulate
824 typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
825 cnt_end = count.end();
826 for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
827 *cnt_it = *cnt_it + *(cnt_it - 1);
828 }
829 // This becomes the array of column pointers
830 trans_ptr.assign(count);
831
832 /* Move the nonzero values into their final place in nzval, based on the
833 * counts found previously.
834 *
835 * This sequence deviates from Knuth's algorithm a bit, following more
836 * closely the description presented in Gustavson, Fred G. "Two Fast
837 * Algorithms for Sparse Matrices: Multiplication and Permuted
838 * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
839 * 250--269, http://doi.acm.org/10.1145/355791.355796.
840 *
841 * The output indices end up in sorted order
842 */
843
844 GlobalSizeT size = ptr.size();
845 for( GlobalSizeT i = 0; i < size - 1; ++i ){
846 GlobalOrdinal u = ptr[i];
847 GlobalOrdinal v = ptr[i + 1];
848 for( GlobalOrdinal j = u; j < v; ++j ){
849 GlobalOrdinal k = count[indices[j]];
850 trans_vals[k] = vals[j];
851 trans_indices[k] = i;
852 ++(count[indices[j]]);
853 }
854 }
855 }
856
857
858 template <typename Scalar1, typename Scalar2>
859 void
860 scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
861 size_t ld, Teuchos::ArrayView<Scalar2> s)
862 {
863 size_t vals_size = vals.size();
864#ifdef HAVE_AMESOS2_DEBUG
865 size_t s_size = s.size();
866 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
867 std::invalid_argument,
868 "Scale vector must have length at least that of the vector" );
869#endif
870 size_t i, s_i;
871 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
872 if( s_i == l ){
873 // bring i to the next multiple of ld
874 i += ld - s_i;
875 s_i = 0;
876 }
877 vals[i] *= s[s_i];
878 }
879 }
880
881 template <typename Scalar1, typename Scalar2, class BinaryOp>
882 void
883 scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
884 size_t ld, Teuchos::ArrayView<Scalar2> s,
885 BinaryOp binary_op)
886 {
887 size_t vals_size = vals.size();
888#ifdef HAVE_AMESOS2_DEBUG
889 size_t s_size = s.size();
890 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
891 std::invalid_argument,
892 "Scale vector must have length at least that of the vector" );
893#endif
894 size_t i, s_i;
895 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
896 if( s_i == l ){
897 // bring i to the next multiple of ld
898 i += ld - s_i;
899 s_i = 0;
900 }
901 vals[i] = binary_op(vals[i], s[s_i]);
902 }
903 }
904
905 template<class row_ptr_view_t, class cols_view_t, class per_view_t>
906 void
907 reorder(row_ptr_view_t & row_ptr, cols_view_t & cols,
908 per_view_t & perm, per_view_t & peri, size_t & nnz,
909 bool permute_matrix)
910 {
911 #ifndef HAVE_AMESOS2_METIS
912 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
913 "Cannot reorder for cuSolver because no METIS is available.");
914 #else
915 typedef typename cols_view_t::value_type ordinal_type;
916 typedef typename row_ptr_view_t::value_type size_type;
917
918 // begin on host where we'll run metis reorder
919 auto host_row_ptr = Kokkos::create_mirror_view(row_ptr);
920 auto host_cols = Kokkos::create_mirror_view(cols);
921 Kokkos::deep_copy(host_row_ptr, row_ptr);
922 Kokkos::deep_copy(host_cols, cols);
923
924 // strip out the diagonals - metis will just crash with them included.
925 // make space for the stripped version
926 typedef Kokkos::View<idx_t*, Kokkos::HostSpace> host_metis_array;
927 const ordinal_type size = row_ptr.size() - 1;
928 size_type max_nnz = host_row_ptr(size);
929 host_metis_array host_strip_diag_row_ptr(
930 Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_row_ptr"),
931 size+1);
932 host_metis_array host_strip_diag_cols(
933 Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_cols"),
934 max_nnz);
935
936 size_type new_nnz = 0;
937 for(ordinal_type i = 0; i < size; ++i) {
938 host_strip_diag_row_ptr(i) = new_nnz;
939 for(size_type j = host_row_ptr(i); j < host_row_ptr(i+1); ++j) {
940 if (i != host_cols(j)) {
941 host_strip_diag_cols(new_nnz++) = host_cols(j);
942 }
943 }
944 }
945 host_strip_diag_row_ptr(size) = new_nnz;
946
947 // we'll get original permutations on host
948 host_metis_array host_perm(
949 Kokkos::ViewAllocateWithoutInitializing("host_perm"), size);
950 host_metis_array host_peri(
951 Kokkos::ViewAllocateWithoutInitializing("host_peri"), size);
952
953 // If we want to remove metis.h included in this header we can move this
954 // to the cpp, but we need to decide how to handle the idx_t declaration.
955 idx_t metis_size = size;
956 int err = METIS_NodeND(&metis_size, host_strip_diag_row_ptr.data(), host_strip_diag_cols.data(),
957 NULL, NULL, host_perm.data(), host_peri.data());
958
959 TEUCHOS_TEST_FOR_EXCEPTION(err != METIS_OK, std::runtime_error,
960 "METIS_NodeND failed to sort matrix.");
961
962 // put the permutations on our saved device ptrs
963 // these will be used to permute x and b when we solve
964 typedef typename cols_view_t::execution_space exec_space_t;
965 auto device_perm = Kokkos::create_mirror_view(exec_space_t(), host_perm);
966 auto device_peri = Kokkos::create_mirror_view(exec_space_t(), host_peri);
967 deep_copy(device_perm, host_perm);
968 deep_copy(device_peri, host_peri);
969
970 // also set the permutation which may need to convert the type from
971 // metis to the native ordinal_type
972 deep_copy_or_assign_view(perm, device_perm);
973 deep_copy_or_assign_view(peri, device_peri);
974
975 if (permute_matrix) {
976 // we'll permute matrix on device to a new set of arrays
977 row_ptr_view_t new_row_ptr(
978 Kokkos::ViewAllocateWithoutInitializing("new_row_ptr"), row_ptr.size());
979 cols_view_t new_cols(
980 Kokkos::ViewAllocateWithoutInitializing("new_cols"), cols.size() - new_nnz/2);
981
982 // permute row indices
983 Kokkos::RangePolicy<exec_space_t> policy_row(0, row_ptr.size());
984 Kokkos::parallel_scan(policy_row, KOKKOS_LAMBDA(
985 ordinal_type i, size_type & update, const bool &final) {
986 if(final) {
987 new_row_ptr(i) = update;
988 }
989 if(i < size) {
990 ordinal_type count = 0;
991 const ordinal_type row = device_perm(i);
992 for(ordinal_type k = row_ptr(row); k < row_ptr(row + 1); ++k) {
993 const ordinal_type j = device_peri(cols(k));
994 count += (i >= j);
995 }
996 update += count;
997 }
998 });
999
1000 // permute col indices
1001 Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1002 Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1003 const ordinal_type kbeg = new_row_ptr(i);
1004 const ordinal_type row = device_perm(i);
1005 const ordinal_type col_beg = row_ptr(row);
1006 const ordinal_type col_end = row_ptr(row + 1);
1007 const ordinal_type nk = col_end - col_beg;
1008 for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1009 const ordinal_type tk = kbeg + t;
1010 const ordinal_type sk = col_beg + k;
1011 const ordinal_type j = device_peri(cols(sk));
1012 if(i >= j) {
1013 new_cols(tk) = j;
1014 ++t;
1015 }
1016 }
1017 });
1018
1019 // finally set the inputs to the new sorted arrays
1020 row_ptr = new_row_ptr;
1021 cols = new_cols;
1022 }
1023
1024 nnz = new_nnz;
1025 #endif // HAVE_AMESOS2_METIS
1026 }
1027
1028 template<class values_view_t, class row_ptr_view_t,
1029 class cols_view_t, class per_view_t>
1030 void
1031 reorder_values(values_view_t & values, const row_ptr_view_t & orig_row_ptr,
1032 const row_ptr_view_t & new_row_ptr,
1033 const cols_view_t & orig_cols, const per_view_t & perm, const per_view_t & peri,
1034 size_t nnz)
1035 {
1036 typedef typename cols_view_t::value_type ordinal_type;
1037 typedef typename cols_view_t::execution_space exec_space_t;
1038
1039 auto device_perm = Kokkos::create_mirror_view(exec_space_t(), perm);
1040 auto device_peri = Kokkos::create_mirror_view(exec_space_t(), peri);
1041 deep_copy(device_perm, perm);
1042 deep_copy(device_peri, peri);
1043
1044 const ordinal_type size = orig_row_ptr.size() - 1;
1045
1046 auto host_orig_row_ptr = Kokkos::create_mirror_view(orig_row_ptr);
1047 auto new_nnz = host_orig_row_ptr(size); // TODO: Maybe optimize this by caching
1048
1049 values_view_t new_values(
1050 Kokkos::ViewAllocateWithoutInitializing("new_values"), values.size() - new_nnz/2);
1051
1052 // permute col indices
1053 Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1054 Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1055 const ordinal_type kbeg = new_row_ptr(i);
1056 const ordinal_type row = device_perm(i);
1057 const ordinal_type col_beg = orig_row_ptr(row);
1058 const ordinal_type col_end = orig_row_ptr(row + 1);
1059 const ordinal_type nk = col_end - col_beg;
1060 for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1061 const ordinal_type tk = kbeg + t;
1062 const ordinal_type sk = col_beg + k;
1063 const ordinal_type j = device_peri(orig_cols(sk));
1064 if(i >= j) {
1065 new_values(tk) = values(sk);
1066 ++t;
1067 }
1068 }
1069 });
1070
1071 values = new_values;
1072 }
1073
1074 template<class array_view_t, class per_view_t>
1075 void
1076 apply_reorder_permutation(const array_view_t & array,
1077 array_view_t & permuted_array, const per_view_t & permutation) {
1078 if(permuted_array.extent(0) != array.extent(0) || permuted_array.extent(1) != array.extent(1)) {
1079 permuted_array = array_view_t(
1080 Kokkos::ViewAllocateWithoutInitializing("permuted_array"),
1081 array.extent(0), array.extent(1));
1082 }
1083 typedef typename array_view_t::execution_space exec_space_t;
1084 Kokkos::RangePolicy<exec_space_t> policy(0, array.extent(0));
1085 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(size_t i) {
1086 for(size_t j = 0; j < array.extent(1); ++j) {
1087 permuted_array(i, j) = array(permutation(i), j);
1088 }
1089 });
1090 }
1091
1094 } // end namespace Util
1095
1096} // end namespace Amesos2
1097
1098#endif // #ifndef AMESOS2_UTIL_HPP
Copy or assign views based on memory spaces.
Provides some simple meta-programming utilities for Amesos2.
Enum and other types declarations for Amesos2.
@ DISTRIBUTED
Definition Amesos2_TypeDecl.hpp:124
@ GLOBALLY_REPLICATED
Definition Amesos2_TypeDecl.hpp:126
@ DISTRIBUTED_NO_OVERLAP
Definition Amesos2_TypeDecl.hpp:125
@ ROOTED
Definition Amesos2_TypeDecl.hpp:127
@ CONTIGUOUS_AND_ROOTED
Definition Amesos2_TypeDecl.hpp:128
EDistribution
Definition Amesos2_TypeDecl.hpp:123
EStorage_Ordering
Definition Amesos2_TypeDecl.hpp:141
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition Amesos2_Util.cpp:119
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:762
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
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 RCP< const Teuchos::Comm< int > > to_teuchos_comm(RCP< const Epetra_Comm > c)
Transform an Epetra_Comm object into a Teuchos::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:674
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:729
A generic helper class for getting a CCS representation of a Matrix.
Definition Amesos2_Util.hpp:652
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition Amesos2_Util.hpp:663
A generic base class for the CRS and CCS helpers.
Definition Amesos2_Util.hpp:271