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"
69#include "Amesos2_Kokkos_View_Copy_Assign.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#ifndef DOXYGEN_SHOULD_SKIP_THIS
255 /*
256 * The following represents a general way of getting a CRS or CCS
257 * representation of a matrix with implicit type conversions. The
258 * \c get_crs_helper and \c get_ccs_helper classes are templated
259 * on 4 types:
260 *
261 * - A Matrix type (conforming to the Amesos2 MatrixAdapter interface)
262 * - A scalar type
263 * - A global ordinal type
264 * - A global size type
265 *
266 * The last three template types correspond to the input argument
267 * types. For example, if the scalar type is \c double , then we
268 * require that the \c nzvals argument is a \c
269 * Teuchos::ArrayView<double> type.
270 *
271 * These helpers perform any type conversions that must be
272 * performed to go between the Matrix's types and the input types.
273 * If no conversions are necessary the static functions can be
274 * effectively inlined.
275 */
276
277 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
278 template <class M, typename S, typename GO, typename GS, class Op>
279 struct same_gs_helper
280 {
281 static void do_get(const Teuchos::Ptr<const M> mat,
282 const ArrayView<typename M::scalar_t> nzvals,
283 const ArrayView<typename M::global_ordinal_t> indices,
284 const ArrayView<GS> pointers,
285 GS& nnz,
286 const Teuchos::Ptr<
287 const Tpetra::Map<typename M::local_ordinal_t,
288 typename M::global_ordinal_t,
289 typename M::node_t> > map,
290 EDistribution distribution,
291 EStorage_Ordering ordering)
292 {
293 Op::apply(mat, nzvals, indices, pointers, nnz, map, distribution, ordering);
294 }
295 };
296
297 template <class M, typename S, typename GO, typename GS, class Op>
298 struct diff_gs_helper
299 {
300 static void do_get(const Teuchos::Ptr<const M> mat,
301 const ArrayView<typename M::scalar_t> nzvals,
302 const ArrayView<typename M::global_ordinal_t> indices,
303 const ArrayView<GS> pointers,
304 GS& nnz,
305 const Teuchos::Ptr<
306 const Tpetra::Map<typename M::local_ordinal_t,
307 typename M::global_ordinal_t,
308 typename M::node_t> > map,
309 EDistribution distribution,
310 EStorage_Ordering ordering)
311 {
312 typedef typename M::global_size_t mat_gs_t;
313 typename ArrayView<GS>::size_type i, size = pointers.size();
314 Teuchos::Array<mat_gs_t> pointers_tmp(size);
315 mat_gs_t nnz_tmp = 0;
316
317 Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, distribution, ordering);
318
319 for (i = 0; i < size; ++i){
320 pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
321 }
322 nnz = Teuchos::as<GS>(nnz_tmp);
323 }
324 };
325
326 template <class M, typename S, typename GO, typename GS, class Op>
327 struct same_go_helper
328 {
329 static void do_get(const Teuchos::Ptr<const M> mat,
330 const ArrayView<typename M::scalar_t> nzvals,
331 const ArrayView<GO> indices,
332 const ArrayView<GS> pointers,
333 GS& nnz,
334 const Teuchos::Ptr<
335 const Tpetra::Map<typename M::local_ordinal_t,
336 typename M::global_ordinal_t,
337 typename M::node_t> > map,
338 EDistribution distribution,
339 EStorage_Ordering ordering)
340 {
341 typedef typename M::global_size_t mat_gs_t;
342 if_then_else<is_same<GS,mat_gs_t>::value,
343 same_gs_helper<M,S,GO,GS,Op>,
344 diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
345 pointers, nnz, map,
346 distribution, ordering);
347 }
348 };
349
350 template <class M, typename S, typename GO, typename GS, class Op>
351 struct diff_go_helper
352 {
353 static void do_get(const Teuchos::Ptr<const M> mat,
354 const ArrayView<typename M::scalar_t> nzvals,
355 const ArrayView<GO> indices,
356 const ArrayView<GS> pointers,
357 GS& nnz,
358 const Teuchos::Ptr<
359 const Tpetra::Map<typename M::local_ordinal_t,
360 typename M::global_ordinal_t,
361 typename M::node_t> > map,
362 EDistribution distribution,
363 EStorage_Ordering ordering)
364 {
365 typedef typename M::global_ordinal_t mat_go_t;
366 typedef typename M::global_size_t mat_gs_t;
367 typename ArrayView<GO>::size_type i, size = indices.size();
368 Teuchos::Array<mat_go_t> indices_tmp(size);
369
370 if_then_else<is_same<GS,mat_gs_t>::value,
371 same_gs_helper<M,S,GO,GS,Op>,
372 diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
373 pointers, nnz, map,
374 distribution, ordering);
375
376 for (i = 0; i < size; ++i){
377 indices[i] = Teuchos::as<GO>(indices_tmp[i]);
378 }
379 }
380 };
381
382 template <class M, typename S, typename GO, typename GS, class Op>
383 struct same_scalar_helper
384 {
385 static void do_get(const Teuchos::Ptr<const M> mat,
386 const ArrayView<S> nzvals,
387 const ArrayView<GO> indices,
388 const ArrayView<GS> pointers,
389 GS& nnz,
390 const Teuchos::Ptr<
391 const Tpetra::Map<typename M::local_ordinal_t,
392 typename M::global_ordinal_t,
393 typename M::node_t> > map,
394 EDistribution distribution,
395 EStorage_Ordering ordering)
396 {
397 typedef typename M::global_ordinal_t mat_go_t;
398 if_then_else<is_same<GO,mat_go_t>::value,
399 same_go_helper<M,S,GO,GS,Op>,
400 diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
401 pointers, nnz, map,
402 distribution, ordering);
403 }
404 };
405
406 template <class M, typename S, typename GO, typename GS, class Op>
407 struct diff_scalar_helper
408 {
409 static void do_get(const Teuchos::Ptr<const M> mat,
410 const ArrayView<S> nzvals,
411 const ArrayView<GO> indices,
412 const ArrayView<GS> pointers,
413 GS& 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::scalar_t mat_scalar_t;
422 typedef typename M::global_ordinal_t mat_go_t;
423 typename ArrayView<S>::size_type i, size = nzvals.size();
424 Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
425
426 if_then_else<is_same<GO,mat_go_t>::value,
427 same_go_helper<M,S,GO,GS,Op>,
428 diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
429 pointers, nnz, map,
430 distribution, ordering);
431
432 for (i = 0; i < size; ++i){
433 nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
434 }
435 }
436 };
437 #endif
438
439#endif // DOXYGEN_SHOULD_SKIP_THIS
440
453 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
454 template<class Matrix, typename S, typename GO, typename GS, class Op>
455 struct get_cxs_helper
456 {
457 static void do_get(const Teuchos::Ptr<const Matrix> mat,
458 const ArrayView<S> nzvals,
459 const ArrayView<GO> indices,
460 const ArrayView<GS> pointers,
461 GS& nnz,
462 EDistribution distribution,
463 EStorage_Ordering ordering=ARBITRARY,
464 GO indexBase = 0)
465 {
466 typedef typename Matrix::local_ordinal_t lo_t;
467 typedef typename Matrix::global_ordinal_t go_t;
468 typedef typename Matrix::global_size_t gs_t;
469 typedef typename Matrix::node_t node_t;
470
471 const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
472 = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
473 Op::get_dimension(mat),
474 mat->getComm(),
475 indexBase,
476 Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
477 );
478
479 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
480 }
481
486 static void do_get(const Teuchos::Ptr<const Matrix> mat,
487 const ArrayView<S> nzvals,
488 const ArrayView<GO> indices,
489 const ArrayView<GS> pointers,
490 GS& nnz,
491 EDistribution distribution, // Does this one need a distribution argument??
492 EStorage_Ordering ordering=ARBITRARY)
493 {
494 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
495 typename Matrix::global_ordinal_t,
496 typename Matrix::node_t> > map
497 = Op::getMap(mat);
498 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
499 }
500
505 static void do_get(const Teuchos::Ptr<const Matrix> mat,
506 const ArrayView<S> nzvals,
507 const ArrayView<GO> indices,
508 const ArrayView<GS> pointers,
509 GS& nnz,
510 const Teuchos::Ptr<
511 const Tpetra::Map<typename Matrix::local_ordinal_t,
512 typename Matrix::global_ordinal_t,
513 typename Matrix::node_t> > map,
514 EDistribution distribution,
515 EStorage_Ordering ordering=ARBITRARY)
516 {
517 typedef typename Matrix::scalar_t mat_scalar;
518 if_then_else<is_same<mat_scalar,S>::value,
519 same_scalar_helper<Matrix,S,GO,GS,Op>,
520 diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
521 nzvals, indices,
522 pointers, nnz,
523 map,
524 distribution, ordering);
525 }
526 };
527 #endif
528
529
530 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
531 struct same_gs_helper_kokkos_view
532 {
533 static void do_get(const Teuchos::Ptr<const M> mat,
534 KV_S& nzvals,
535 KV_GO& indices,
536 KV_GS& pointers,
537 typename KV_GS::value_type& nnz,
538 const Teuchos::Ptr<
539 const Tpetra::Map<typename M::local_ordinal_t,
540 typename M::global_ordinal_t,
541 typename M::node_t> > map,
542 EDistribution distribution,
543 EStorage_Ordering ordering)
544 {
545 Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
546 indices, pointers, nnz, map, distribution, ordering);
547 }
548 };
549
550 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
551 struct diff_gs_helper_kokkos_view
552 {
553 static void do_get(const Teuchos::Ptr<const M> mat,
554 KV_S& nzvals,
555 KV_GO& indices,
556 KV_GS& pointers,
557 typename KV_GS::value_type& nnz,
558 const Teuchos::Ptr<
559 const Tpetra::Map<typename M::local_ordinal_t,
560 typename M::global_ordinal_t,
561 typename M::node_t> > map,
562 EDistribution distribution,
563 EStorage_Ordering ordering)
564 {
565 typedef typename M::global_size_t mat_gs_t;
566 typedef typename Kokkos::View<mat_gs_t*, Kokkos::HostSpace> KV_TMP;
567 size_t i, size = pointers.extent(0);
568 KV_TMP pointers_tmp(Kokkos::ViewAllocateWithoutInitializing("pointers_tmp"), size);
569
570 mat_gs_t nnz_tmp = 0;
571 Op::template apply_kokkos_view<KV_S, KV_GO, KV_TMP>(mat, nzvals,
572 indices, pointers_tmp, nnz_tmp, Teuchos::ptrInArg(*map), distribution, ordering);
573 nnz = Teuchos::as<typename KV_GS::value_type>(nnz_tmp);
574
575 typedef typename KV_GS::value_type view_gs_t;
576 for (i = 0; i < size; ++i){
577 pointers(i) = Teuchos::as<view_gs_t>(pointers_tmp(i));
578 }
579 nnz = Teuchos::as<view_gs_t>(nnz_tmp);
580 }
581 };
582
583 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
584 struct same_go_helper_kokkos_view
585 {
586 static void do_get(const Teuchos::Ptr<const M> mat,
587 KV_S& nzvals,
588 KV_GO& indices,
589 KV_GS& pointers,
590 typename KV_GS::value_type& nnz,
591 const Teuchos::Ptr<
592 const Tpetra::Map<typename M::local_ordinal_t,
593 typename M::global_ordinal_t,
594 typename M::node_t> > map,
595 EDistribution distribution,
596 EStorage_Ordering ordering)
597 {
598 typedef typename M::global_size_t mat_gs_t;
599 typedef typename KV_GS::value_type view_gs_t;
600 if_then_else<is_same<view_gs_t,mat_gs_t>::value,
601 same_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
602 diff_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
603 pointers, nnz, map,
604 distribution, ordering);
605 }
606 };
607
608 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
609 struct diff_go_helper_kokkos_view
610 {
611 static void do_get(const Teuchos::Ptr<const M> mat,
612 KV_S& nzvals,
613 KV_GO& indices,
614 KV_GS& pointers,
615 typename KV_GS::value_type& nnz,
616 const Teuchos::Ptr<
617 const Tpetra::Map<typename M::local_ordinal_t,
618 typename M::global_ordinal_t,
619 typename M::node_t> > map,
620 EDistribution distribution,
621 EStorage_Ordering ordering)
622 {
623 typedef typename M::global_ordinal_t mat_go_t;
624 typedef typename M::global_size_t mat_gs_t;
625 typedef typename Kokkos::View<mat_go_t*, Kokkos::HostSpace> KV_TMP;
626 size_t i, size = indices.extent(0);
627 KV_TMP indices_tmp(Kokkos::ViewAllocateWithoutInitializing("indices_tmp"), size);
628
629 typedef typename KV_GO::value_type view_go_t;
630 typedef typename KV_GS::value_type view_gs_t;
631 if_then_else<is_same<view_gs_t,mat_gs_t>::value,
632 same_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op>,
633 diff_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op> >::type::do_get(mat, nzvals, indices_tmp,
634 pointers, nnz, map,
635 distribution, ordering);
636 for (i = 0; i < size; ++i){
637 indices(i) = Teuchos::as<view_go_t>(indices_tmp(i));
638 }
639 }
640 };
641
642 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
643 struct same_scalar_helper_kokkos_view
644 {
645 static void do_get(const Teuchos::Ptr<const M> mat,
646 KV_S& nzvals,
647 KV_GO& indices,
648 KV_GS& pointers,
649 typename KV_GS::value_type& nnz,
650 const Teuchos::Ptr<
651 const Tpetra::Map<typename M::local_ordinal_t,
652 typename M::global_ordinal_t,
653 typename M::node_t> > map,
654 EDistribution distribution,
655 EStorage_Ordering ordering)
656 {
657 typedef typename M::global_ordinal_t mat_go_t;
658 typedef typename KV_GO::value_type view_go_t;
659 if_then_else<is_same<view_go_t,mat_go_t>::value,
660 same_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
661 diff_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals, indices,
662 pointers, nnz, map,
663 distribution, ordering);
664 }
665 };
666
667 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
668 struct diff_scalar_helper_kokkos_view
669 {
670 static void do_get(const Teuchos::Ptr<const M> mat,
671 KV_S& nzvals,
672 KV_GO& indices,
673 KV_GS& pointers,
674 typename KV_GS::value_type& nnz,
675 const Teuchos::Ptr<
676 const Tpetra::Map<typename M::local_ordinal_t,
677 typename M::global_ordinal_t,
678 typename M::node_t> > map,
679 EDistribution distribution,
680 EStorage_Ordering ordering)
681 {
682 typedef typename M::global_ordinal_t mat_go_t;
683 typedef typename Kokkos::ArithTraits<typename M::scalar_t>::val_type mat_scalar_t;
684 typedef typename Kokkos::View<mat_scalar_t*, Kokkos::HostSpace> KV_TMP;
685 size_t i, size = nzvals.extent(0);
686 KV_TMP nzvals_tmp(Kokkos::ViewAllocateWithoutInitializing("nzvals_tmp"), size);
687
688 typedef typename KV_S::value_type view_scalar_t;
689 typedef typename KV_GO::value_type view_go_t;
690 if_then_else<is_same<view_go_t,mat_go_t>::value,
691 same_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op>,
692 diff_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op> >::type::do_get(mat, nzvals_tmp, indices,
693 pointers, nnz, map,
694 distribution, ordering);
695
696 for (i = 0; i < size; ++i){
697 nzvals(i) = Teuchos::as<view_scalar_t>(nzvals_tmp(i));
698 }
699 }
700 };
701
702
703 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS, class Op>
704 struct get_cxs_helper_kokkos_view
705 {
706 static void do_get(const Teuchos::Ptr<const Matrix> mat,
707 KV_S& nzvals,
708 KV_GO& indices,
709 KV_GS& pointers,
710 typename KV_GS::value_type& nnz,
711 EDistribution distribution,
712 EStorage_Ordering ordering=ARBITRARY,
713 typename KV_GO::value_type indexBase = 0)
714 {
715 typedef typename Matrix::local_ordinal_t lo_t;
716 typedef typename Matrix::global_ordinal_t go_t;
717 typedef typename Matrix::global_size_t gs_t;
718 typedef typename Matrix::node_t node_t;
719
720 const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
721 = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
722 Op::get_dimension(mat),
723 mat->getComm(),
724 indexBase,
725 Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
726 );
727 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
728 }
729
734 static void do_get(const Teuchos::Ptr<const Matrix> mat,
735 KV_S& nzvals,
736 KV_GO& indices,
737 KV_GS& pointers,
738 typename KV_GS::value_type& nnz,
739 EDistribution distribution, // Does this one need a distribution argument??
740 EStorage_Ordering ordering=ARBITRARY)
741 {
742 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
743 typename Matrix::global_ordinal_t,
744 typename Matrix::node_t> > map
745 = Op::getMap(mat);
746 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
747 }
748
753 static void do_get(const Teuchos::Ptr<const Matrix> mat,
754 KV_S& nzvals,
755 KV_GO& indices,
756 KV_GS& pointers,
757 typename KV_GS::value_type& nnz,
758 const Teuchos::Ptr<
759 const Tpetra::Map<typename Matrix::local_ordinal_t,
760 typename Matrix::global_ordinal_t,
761 typename Matrix::node_t> > map,
762 EDistribution distribution,
763 EStorage_Ordering ordering=ARBITRARY)
764 {
765 typedef typename Matrix::scalar_t mat_scalar;
766 typedef typename KV_S::value_type view_scalar_t;
767
768 if_then_else<is_same<mat_scalar,view_scalar_t>::value,
769 same_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op>,
770 diff_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op> >::type::do_get(mat,
771 nzvals, indices,
772 pointers, nnz,
773 map,
774 distribution, ordering);
775 }
776 };
777
778#ifndef DOXYGEN_SHOULD_SKIP_THIS
779 /*
780 * These two function-like classes are meant to be used as the \c
781 * Op template parameter for the \c get_cxs_helper template class.
782 */
783 template<class Matrix>
784 struct get_ccs_func
785 {
786 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
787 static void apply(const Teuchos::Ptr<const Matrix> mat,
788 const ArrayView<typename Matrix::scalar_t> nzvals,
789 const ArrayView<typename Matrix::global_ordinal_t> rowind,
790 const ArrayView<typename Matrix::global_size_t> colptr,
791 typename Matrix::global_size_t& nnz,
792 const Teuchos::Ptr<
793 const Tpetra::Map<typename Matrix::local_ordinal_t,
794 typename Matrix::global_ordinal_t,
795 typename Matrix::node_t> > map,
796 EDistribution distribution,
797 EStorage_Ordering ordering)
798 {
799 mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering, distribution);
800 //mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
801 }
802 #endif
803
804 template<typename KV_S, typename KV_GO, typename KV_GS>
805 static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
806 KV_S& nzvals,
807 KV_GO& rowind,
808 KV_GS& colptr,
809 typename Matrix::global_size_t& nnz,
810 const Teuchos::Ptr<
811 const Tpetra::Map<typename Matrix::local_ordinal_t,
812 typename Matrix::global_ordinal_t,
813 typename Matrix::node_t> > map,
814 EDistribution distribution,
815 EStorage_Ordering ordering)
816 {
817 mat->getCcs_kokkos_view(nzvals, rowind, colptr, nnz, map, ordering, distribution);
818 }
819
820 static
821 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
822 typename Matrix::global_ordinal_t,
823 typename Matrix::node_t> >
824 getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
825 {
826 return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
827 }
828
829 static
830 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
831 typename Matrix::global_ordinal_t,
832 typename Matrix::node_t> >
833 getMap(const Teuchos::Ptr<const Matrix> mat)
834 {
835 return mat->getColMap();
836 }
837
838 static
839 typename Matrix::global_size_t
840 get_dimension(const Teuchos::Ptr<const Matrix> mat)
841 {
842 return mat->getGlobalNumCols();
843 }
844 };
845
846 template<class Matrix>
847 struct get_crs_func
848 {
849 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
850 static void apply(const Teuchos::Ptr<const Matrix> mat,
851 const ArrayView<typename Matrix::scalar_t> nzvals,
852 const ArrayView<typename Matrix::global_ordinal_t> colind,
853 const ArrayView<typename Matrix::global_size_t> rowptr,
854 typename Matrix::global_size_t& nnz,
855 const Teuchos::Ptr<
856 const Tpetra::Map<typename Matrix::local_ordinal_t,
857 typename Matrix::global_ordinal_t,
858 typename Matrix::node_t> > map,
859 EDistribution distribution,
860 EStorage_Ordering ordering)
861 {
862 mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering, distribution);
863 }
864 #endif
865
866 template<typename KV_S, typename KV_GO, typename KV_GS>
867 static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
868 KV_S& nzvals,
869 KV_GO& colind,
870 KV_GS& rowptr,
871 typename Matrix::global_size_t& nnz,
872 const Teuchos::Ptr<
873 const Tpetra::Map<typename Matrix::local_ordinal_t,
874 typename Matrix::global_ordinal_t,
875 typename Matrix::node_t> > map,
876 EDistribution distribution,
877 EStorage_Ordering ordering)
878 {
879 mat->getCrs_kokkos_view(nzvals, colind, rowptr, nnz, map, ordering, distribution);
880 }
881
882 static
883 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
884 typename Matrix::global_ordinal_t,
885 typename Matrix::node_t> >
886 getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
887 {
888 return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
889 }
890
891 static
892 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
893 typename Matrix::global_ordinal_t,
894 typename Matrix::node_t> >
895 getMap(const Teuchos::Ptr<const Matrix> mat)
896 {
897 return mat->getRowMap();
898 }
899
900 static
901 typename Matrix::global_size_t
902 get_dimension(const Teuchos::Ptr<const Matrix> mat)
903 {
904 return mat->getGlobalNumRows();
905 }
906 };
907#endif // DOXYGEN_SHOULD_SKIP_THIS
908
946 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
947 template<class Matrix, typename S, typename GO, typename GS>
948 struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
949 {};
950 #endif
951
952 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
953 struct get_ccs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_ccs_func<Matrix> >
954 {};
955
963 #ifdef TPETRA_ENABLE_DEPRECATED_CODE
964 template<class Matrix, typename S, typename GO, typename GS>
965 struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
966 {};
967 #endif
968
969 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
970 struct get_crs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_crs_func<Matrix> >
971 {};
972 /* End Matrix/MultiVector Utilities */
973
974
976 // Definitions //
978
979
980 template <typename LO, typename GO, typename GS, typename Node>
981 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
982 getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
983 {
984 //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
985 Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
986 return gather_map;
987 }
988
989
990 template <typename LO, typename GO, typename GS, typename Node>
991 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
992 getDistributionMap(EDistribution distribution,
993 GS num_global_elements,
994 const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
995 GO indexBase,
996 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
997 {
998 // TODO: Need to add indexBase to cases other than ROOTED
999 // We do not support these maps in any solver now.
1000 switch( distribution ){
1001 case DISTRIBUTED:
1002 case DISTRIBUTED_NO_OVERLAP:
1003 return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
1004 case GLOBALLY_REPLICATED:
1005 return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
1006 case ROOTED:
1007 {
1008 int rank = Teuchos::rank(*comm);
1009 size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
1010 if( rank == 0 ) { my_num_elems = num_global_elements; }
1011
1012 return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
1013 my_num_elems, indexBase, comm));
1014 }
1015 case CONTIGUOUS_AND_ROOTED:
1016 {
1017 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
1018 = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
1019 return gathermap;
1020 }
1021 default:
1022 TEUCHOS_TEST_FOR_EXCEPTION( true,
1023 std::logic_error,
1024 "Control should never reach this point. "
1025 "Please contact the Amesos2 developers." );
1026 }
1027 }
1028
1029
1030#ifdef HAVE_AMESOS2_EPETRA
1031
1032 //#pragma message "include 3"
1033 //#include <Epetra_Map.h>
1034
1035 template <typename LO, typename GO, typename GS, typename Node>
1036 Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
1037 epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
1038 {
1039 using Teuchos::as;
1040 using Teuchos::rcp;
1041
1042 int num_my_elements = map.NumMyElements();
1043 Teuchos::Array<int> my_global_elements(num_my_elements);
1044 map.MyGlobalElements(my_global_elements.getRawPtr());
1045
1046 Teuchos::Array<GO> my_gbl_inds_buf;
1047 Teuchos::ArrayView<GO> my_gbl_inds;
1048 if (! std::is_same<int, GO>::value) {
1049 my_gbl_inds_buf.resize (num_my_elements);
1050 my_gbl_inds = my_gbl_inds_buf ();
1051 for (int k = 0; k < num_my_elements; ++k) {
1052 my_gbl_inds[k] = static_cast<GO> (my_global_elements[k]);
1053 }
1054 }
1055 else {
1056 using Teuchos::av_reinterpret_cast;
1057 my_gbl_inds = av_reinterpret_cast<GO> (my_global_elements ());
1058 }
1059
1060 typedef Tpetra::Map<LO,GO,Node> map_t;
1061 RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
1062 my_gbl_inds(),
1063 as<GO>(map.IndexBase()),
1064 to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
1065 return tmap;
1066 }
1067
1068 template <typename LO, typename GO, typename GS, typename Node>
1069 Teuchos::RCP<Epetra_Map>
1070 tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
1071 {
1072 using Teuchos::as;
1073
1074 Teuchos::Array<GO> elements_tmp;
1075 elements_tmp = map.getNodeElementList();
1076 int num_my_elements = elements_tmp.size();
1077 Teuchos::Array<int> my_global_elements(num_my_elements);
1078 for (int i = 0; i < num_my_elements; ++i){
1079 my_global_elements[i] = as<int>(elements_tmp[i]);
1080 }
1081
1082 using Teuchos::rcp;
1083 RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
1084 num_my_elements,
1085 my_global_elements.getRawPtr(),
1086 as<GO>(map.getIndexBase()),
1087 *to_epetra_comm(map.getComm())));
1088 return emap;
1089 }
1090#endif // HAVE_AMESOS2_EPETRA
1091
1092 template <typename Scalar,
1093 typename GlobalOrdinal,
1094 typename GlobalSizeT>
1095 void transpose(Teuchos::ArrayView<Scalar> vals,
1096 Teuchos::ArrayView<GlobalOrdinal> indices,
1097 Teuchos::ArrayView<GlobalSizeT> ptr,
1098 Teuchos::ArrayView<Scalar> trans_vals,
1099 Teuchos::ArrayView<GlobalOrdinal> trans_indices,
1100 Teuchos::ArrayView<GlobalSizeT> trans_ptr)
1101 {
1102 /* We have a compressed-row storage format of this matrix. We
1103 * transform this into a compressed-column format using a
1104 * distribution-counting sort algorithm, which is described by
1105 * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
1106 */
1107
1108#ifdef HAVE_AMESOS2_DEBUG
1109 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
1110 ind_begin = indices.begin();
1111 ind_end = indices.end();
1112 size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
1113 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
1114 std::invalid_argument,
1115 "Transpose pointer size not large enough." );
1116 TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
1117 std::invalid_argument,
1118 "Transpose values array not large enough." );
1119 TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
1120 std::invalid_argument,
1121 "Transpose indices array not large enough." );
1122#else
1123 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
1124#endif
1125 // Count the number of entries in each column
1126 Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
1127 ind_end = indices.end();
1128 for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
1129 ++(count[(*ind_it) + 1]);
1130 }
1131 // Accumulate
1132 typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
1133 cnt_end = count.end();
1134 for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
1135 *cnt_it = *cnt_it + *(cnt_it - 1);
1136 }
1137 // This becomes the array of column pointers
1138 trans_ptr.assign(count);
1139
1140 /* Move the nonzero values into their final place in nzval, based on the
1141 * counts found previously.
1142 *
1143 * This sequence deviates from Knuth's algorithm a bit, following more
1144 * closely the description presented in Gustavson, Fred G. "Two Fast
1145 * Algorithms for Sparse Matrices: Multiplication and Permuted
1146 * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
1147 * 250--269, http://doi.acm.org/10.1145/355791.355796.
1148 *
1149 * The output indices end up in sorted order
1150 */
1151
1152 GlobalSizeT size = ptr.size();
1153 for( GlobalSizeT i = 0; i < size - 1; ++i ){
1154 GlobalOrdinal u = ptr[i];
1155 GlobalOrdinal v = ptr[i + 1];
1156 for( GlobalOrdinal j = u; j < v; ++j ){
1157 GlobalOrdinal k = count[indices[j]];
1158 trans_vals[k] = vals[j];
1159 trans_indices[k] = i;
1160 ++(count[indices[j]]);
1161 }
1162 }
1163 }
1164
1165
1166 template <typename Scalar1, typename Scalar2>
1167 void
1168 scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
1169 size_t ld, Teuchos::ArrayView<Scalar2> s)
1170 {
1171 size_t vals_size = vals.size();
1172#ifdef HAVE_AMESOS2_DEBUG
1173 size_t s_size = s.size();
1174 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
1175 std::invalid_argument,
1176 "Scale vector must have length at least that of the vector" );
1177#endif
1178 size_t i, s_i;
1179 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
1180 if( s_i == l ){
1181 // bring i to the next multiple of ld
1182 i += ld - s_i;
1183 s_i = 0;
1184 }
1185 vals[i] *= s[s_i];
1186 }
1187 }
1188
1189 template <typename Scalar1, typename Scalar2, class BinaryOp>
1190 void
1191 scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
1192 size_t ld, Teuchos::ArrayView<Scalar2> s,
1193 BinaryOp binary_op)
1194 {
1195 size_t vals_size = vals.size();
1196#ifdef HAVE_AMESOS2_DEBUG
1197 size_t s_size = s.size();
1198 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
1199 std::invalid_argument,
1200 "Scale vector must have length at least that of the vector" );
1201#endif
1202 size_t i, s_i;
1203 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
1204 if( s_i == l ){
1205 // bring i to the next multiple of ld
1206 i += ld - s_i;
1207 s_i = 0;
1208 }
1209 vals[i] = binary_op(vals[i], s[s_i]);
1210 }
1211 }
1212
1213 template<class row_ptr_view_t, class cols_view_t, class per_view_t>
1214 void
1215 reorder(row_ptr_view_t & row_ptr, cols_view_t & cols,
1216 per_view_t & perm, per_view_t & peri, size_t & nnz,
1217 bool permute_matrix)
1218 {
1219 #ifndef HAVE_AMESOS2_METIS
1220 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
1221 "Cannot reorder for cuSolver because no METIS is available.");
1222 #else
1223 typedef typename cols_view_t::value_type ordinal_type;
1224 typedef typename row_ptr_view_t::value_type size_type;
1225
1226 // begin on host where we'll run metis reorder
1227 auto host_row_ptr = Kokkos::create_mirror_view(row_ptr);
1228 auto host_cols = Kokkos::create_mirror_view(cols);
1229 Kokkos::deep_copy(host_row_ptr, row_ptr);
1230 Kokkos::deep_copy(host_cols, cols);
1231
1232 // strip out the diagonals - metis will just crash with them included.
1233 // make space for the stripped version
1234 typedef Kokkos::View<idx_t*, Kokkos::HostSpace> host_metis_array;
1235 const ordinal_type size = row_ptr.size() - 1;
1236 size_type max_nnz = host_row_ptr(size);
1237 host_metis_array host_strip_diag_row_ptr(
1238 Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_row_ptr"),
1239 size+1);
1240 host_metis_array host_strip_diag_cols(
1241 Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_cols"),
1242 max_nnz);
1243
1244 size_type new_nnz = 0;
1245 for(ordinal_type i = 0; i < size; ++i) {
1246 host_strip_diag_row_ptr(i) = new_nnz;
1247 for(size_type j = host_row_ptr(i); j < host_row_ptr(i+1); ++j) {
1248 if (i != host_cols(j)) {
1249 host_strip_diag_cols(new_nnz++) = host_cols(j);
1250 }
1251 }
1252 }
1253 host_strip_diag_row_ptr(size) = new_nnz;
1254
1255 // we'll get original permutations on host
1256 host_metis_array host_perm(
1257 Kokkos::ViewAllocateWithoutInitializing("host_perm"), size);
1258 host_metis_array host_peri(
1259 Kokkos::ViewAllocateWithoutInitializing("host_peri"), size);
1260
1261 // If we want to remove metis.h included in this header we can move this
1262 // to the cpp, but we need to decide how to handle the idx_t declaration.
1263 idx_t metis_size = size;
1264 int err = METIS_NodeND(&metis_size, host_strip_diag_row_ptr.data(), host_strip_diag_cols.data(),
1265 NULL, NULL, host_perm.data(), host_peri.data());
1266
1267 TEUCHOS_TEST_FOR_EXCEPTION(err != METIS_OK, std::runtime_error,
1268 "METIS_NodeND failed to sort matrix.");
1269
1270 // put the permutations on our saved device ptrs
1271 // these will be used to permute x and b when we solve
1272 typedef typename cols_view_t::execution_space exec_space_t;
1273 auto device_perm = Kokkos::create_mirror_view(exec_space_t(), host_perm);
1274 auto device_peri = Kokkos::create_mirror_view(exec_space_t(), host_peri);
1275 deep_copy(device_perm, host_perm);
1276 deep_copy(device_peri, host_peri);
1277
1278 // also set the permutation which may need to convert the type from
1279 // metis to the native ordinal_type
1280 deep_copy_or_assign_view(perm, device_perm);
1281 deep_copy_or_assign_view(peri, device_peri);
1282
1283 if (permute_matrix) {
1284 // we'll permute matrix on device to a new set of arrays
1285 row_ptr_view_t new_row_ptr(
1286 Kokkos::ViewAllocateWithoutInitializing("new_row_ptr"), row_ptr.size());
1287 cols_view_t new_cols(
1288 Kokkos::ViewAllocateWithoutInitializing("new_cols"), cols.size() - new_nnz/2);
1289
1290 // permute row indices
1291 Kokkos::RangePolicy<exec_space_t> policy_row(0, row_ptr.size());
1292 Kokkos::parallel_scan(policy_row, KOKKOS_LAMBDA(
1293 ordinal_type i, size_type & update, const bool &final) {
1294 if(final) {
1295 new_row_ptr(i) = update;
1296 }
1297 if(i < size) {
1298 ordinal_type count = 0;
1299 const ordinal_type row = device_perm(i);
1300 for(ordinal_type k = row_ptr(row); k < row_ptr(row + 1); ++k) {
1301 const ordinal_type j = device_peri(cols(k));
1302 count += (i >= j);
1303 }
1304 update += count;
1305 }
1306 });
1307
1308 // permute col indices
1309 Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1310 Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1311 const ordinal_type kbeg = new_row_ptr(i);
1312 const ordinal_type row = device_perm(i);
1313 const ordinal_type col_beg = row_ptr(row);
1314 const ordinal_type col_end = row_ptr(row + 1);
1315 const ordinal_type nk = col_end - col_beg;
1316 for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1317 const ordinal_type tk = kbeg + t;
1318 const ordinal_type sk = col_beg + k;
1319 const ordinal_type j = device_peri(cols(sk));
1320 if(i >= j) {
1321 new_cols(tk) = j;
1322 ++t;
1323 }
1324 }
1325 });
1326
1327 // finally set the inputs to the new sorted arrays
1328 row_ptr = new_row_ptr;
1329 cols = new_cols;
1330 }
1331
1332 nnz = new_nnz;
1333 #endif // HAVE_AMESOS2_METIS
1334 }
1335
1336 template<class values_view_t, class row_ptr_view_t,
1337 class cols_view_t, class per_view_t>
1338 void
1339 reorder_values(values_view_t & values, const row_ptr_view_t & orig_row_ptr,
1340 const row_ptr_view_t & new_row_ptr,
1341 const cols_view_t & orig_cols, const per_view_t & perm, const per_view_t & peri,
1342 size_t nnz)
1343 {
1344 typedef typename cols_view_t::value_type ordinal_type;
1345 typedef typename cols_view_t::execution_space exec_space_t;
1346
1347 auto device_perm = Kokkos::create_mirror_view(exec_space_t(), perm);
1348 auto device_peri = Kokkos::create_mirror_view(exec_space_t(), peri);
1349 deep_copy(device_perm, perm);
1350 deep_copy(device_peri, peri);
1351
1352 const ordinal_type size = orig_row_ptr.size() - 1;
1353
1354 auto host_orig_row_ptr = Kokkos::create_mirror_view(orig_row_ptr);
1355 auto new_nnz = host_orig_row_ptr(size); // TODO: Maybe optimize this by caching
1356
1357 values_view_t new_values(
1358 Kokkos::ViewAllocateWithoutInitializing("new_values"), values.size() - new_nnz/2);
1359
1360 // permute col indices
1361 Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
1362 Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
1363 const ordinal_type kbeg = new_row_ptr(i);
1364 const ordinal_type row = device_perm(i);
1365 const ordinal_type col_beg = orig_row_ptr(row);
1366 const ordinal_type col_end = orig_row_ptr(row + 1);
1367 const ordinal_type nk = col_end - col_beg;
1368 for(ordinal_type k = 0, t = 0; k < nk; ++k) {
1369 const ordinal_type tk = kbeg + t;
1370 const ordinal_type sk = col_beg + k;
1371 const ordinal_type j = device_peri(orig_cols(sk));
1372 if(i >= j) {
1373 new_values(tk) = values(sk);
1374 ++t;
1375 }
1376 }
1377 });
1378
1379 values = new_values;
1380 }
1381
1382 template<class array_view_t, class per_view_t>
1383 void
1384 apply_reorder_permutation(const array_view_t & array,
1385 array_view_t & permuted_array, const per_view_t & permutation) {
1386 if(permuted_array.extent(0) != array.extent(0) || permuted_array.extent(1) != array.extent(1)) {
1387 permuted_array = array_view_t(
1388 Kokkos::ViewAllocateWithoutInitializing("permuted_array"),
1389 array.extent(0), array.extent(1));
1390 }
1391 typedef typename array_view_t::execution_space exec_space_t;
1392 Kokkos::RangePolicy<exec_space_t> policy(0, array.extent(0));
1393 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(size_t i) {
1394 for(size_t j = 0; j < array.extent(1); ++j) {
1395 permuted_array(i, j) = array(permutation(i), j);
1396 }
1397 });
1398 }
1399
1402 } // end namespace Util
1403
1404} // end namespace Amesos2
1405
1406#endif // #ifndef AMESOS2_UTIL_HPP
Amesos2_Kokkos_View_Copy_Assign.hpp
Copy or assign views based on memory spaces.
Amesos2_Meta.hpp
Provides some simple meta-programming utilities for Amesos2.
Amesos2_TypeDecl.hpp
Enum and other types declarations for Amesos2.
Amesos2::DISTRIBUTED
@ DISTRIBUTED
Definition: Amesos2_TypeDecl.hpp:124
Amesos2::GLOBALLY_REPLICATED
@ GLOBALLY_REPLICATED
Definition: Amesos2_TypeDecl.hpp:126
Amesos2::DISTRIBUTED_NO_OVERLAP
@ DISTRIBUTED_NO_OVERLAP
Definition: Amesos2_TypeDecl.hpp:125
Amesos2::ROOTED
@ ROOTED
Definition: Amesos2_TypeDecl.hpp:127
Amesos2::CONTIGUOUS_AND_ROOTED
@ CONTIGUOUS_AND_ROOTED
Definition: Amesos2_TypeDecl.hpp:128
Amesos2::EDistribution
EDistribution
Definition: Amesos2_TypeDecl.hpp:123
Amesos2::EStorage_Ordering
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
Amesos2::Util::printLine
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition: Amesos2_Util.cpp:119
Amesos2::Util::tpetra_map_to_epetra_map
Teuchos::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:1070
Amesos2::Util::scale
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
Amesos2::Util::transpose
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
Amesos2::Util::to_epetra_comm
const RCP< const Epetra_Comm > to_epetra_comm(RCP< const Teuchos::Comm< int > > c)
Transfrom a Teuchos::Comm object into an Epetra_Comm object.
Amesos2::Util::to_teuchos_comm
const RCP< const Teuchos::Comm< int > > to_teuchos_comm(RCP< const Epetra_Comm > c)
Transform an Epetra_Comm object into a Teuchos::Comm object.
Amesos2::Util::getGatherMap
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:982
Amesos2::Util::epetra_map_to_tpetra_map
Teuchos::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:1037
Amesos2::Util::get_ccs_helper_kokkos_view
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:954
Amesos2::Util::get_crs_helper_kokkos_view
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:971
Amesos2::Util::same_gs_helper_kokkos_view
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:532
Generated by doxygen 1.9.4