Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Util.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Amesos2: Templated Direct Sparse Solver Package
4//
5// Copyright 2011 NTESS and the Amesos2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
17
18#ifndef AMESOS2_UTIL_HPP
19#define AMESOS2_UTIL_HPP
20
21#include <cstdio>
22#include <fstream>
23#include <iostream>
24
25#include "Amesos2_config.h"
26
27#include "Teuchos_RCP.hpp"
28#include "Teuchos_BLAS_types.hpp"
29#include "Teuchos_Array.hpp"
30#include "Teuchos_ArrayView.hpp"
31#include "Teuchos_FancyOStream.hpp"
32
33#include <Tpetra_Map.hpp>
34#include <Tpetra_DistObject_decl.hpp>
35#include <Tpetra_ComputeGatherMap.hpp> // added for gather map... where is the best place??
36
37#include "Amesos2_TypeDecl.hpp"
38#include "Amesos2_Meta.hpp"
39#include "Amesos2_Kokkos_View_Copy_Assign.hpp"
40
41#ifdef HAVE_AMESOS2_METIS
42#include "metis.h" // to discuss, remove from header?
43#endif
44
45namespace Amesos2 {
46
47 namespace Util {
48
54
55 using Teuchos::RCP;
56 using Teuchos::ArrayView;
57
73
74 template <typename LO, typename GO, typename GS, typename Node>
75 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
76 getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map );
77
78
79 template <typename LO, typename GO, typename GS, typename Node>
80 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
81 getDistributionMap(EDistribution distribution,
82 GS num_global_elements,
83 const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
84 GO indexBase = 0,
85 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map = Teuchos::null);
86
92 template <typename Scalar,
93 typename GlobalOrdinal,
94 typename GlobalSizeT>
95 void transpose(ArrayView<Scalar> vals,
96 ArrayView<GlobalOrdinal> indices,
97 ArrayView<GlobalSizeT> ptr,
98 ArrayView<Scalar> trans_vals,
99 ArrayView<GlobalOrdinal> trans_indices,
100 ArrayView<GlobalSizeT> trans_ptr);
101
115 template <typename Scalar1, typename Scalar2>
116 void scale(ArrayView<Scalar1> vals, size_t l,
117 size_t ld, ArrayView<Scalar2> s);
118
137 template <typename Scalar1, typename Scalar2, class BinaryOp>
138 void scale(ArrayView<Scalar1> vals, size_t l,
139 size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
140
141
143 void printLine( Teuchos::FancyOStream &out );
144
145 // Helper function used to convert Kokkos::complex pointer
146 // to std::complex pointer; needed for optimized code path
147 // when retrieving the CRS raw pointers
148 template < class T0, class T1 >
149 struct getStdCplxType
150 {
151 using common_type = typename std::common_type<T0,T1>::type;
152 using type = common_type;
153 };
154
155 template < class T0, class T1 >
156 struct getStdCplxType< T0, T1* >
157 {
158 using common_type = typename std::common_type<T0,T1>::type;
159 using type = common_type;
160 };
161
162#if defined(HAVE_TEUCHOS_COMPLEX) && defined(HAVE_AMESOS2_KOKKOS)
163 template < class T0 >
164 struct getStdCplxType< T0, Kokkos::complex<T0>* >
165 {
166 using type = std::complex<T0>;
167 };
168
169 template < class T0 , class T1 >
170 struct getStdCplxType< T0, Kokkos::complex<T1>* >
171 {
172 using common_type = typename std::common_type<T0,T1>::type;
173 using type = std::complex<common_type>;
174 };
175#endif
176
178 // Matrix/MultiVector Utilities //
180
181
182
195
196 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
197 struct same_gs_helper_kokkos_view
198 {
199 static void do_get(const Teuchos::Ptr<const M> mat,
200 KV_S& nzvals,
201 KV_GO& indices,
202 KV_GS& pointers,
203 typename KV_GS::value_type& nnz,
204 const Teuchos::Ptr<
205 const Tpetra::Map<typename M::local_ordinal_t,
206 typename M::global_ordinal_t,
207 typename M::node_t> > map,
208 EDistribution distribution,
209 EStorage_Ordering ordering)
210 {
211 Op::template apply_kokkos_view<KV_S, KV_GO, KV_GS>(mat, nzvals,
212 indices, pointers, nnz, map, distribution, ordering);
213 }
214 };
215
216 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
217 struct diff_gs_helper_kokkos_view
218 {
219 static void do_get(const Teuchos::Ptr<const M> mat,
220 KV_S& nzvals,
221 KV_GO& indices,
222 KV_GS& pointers,
223 typename KV_GS::value_type& nnz,
224 const Teuchos::Ptr<
225 const Tpetra::Map<typename M::local_ordinal_t,
226 typename M::global_ordinal_t,
227 typename M::node_t> > map,
228 EDistribution distribution,
229 EStorage_Ordering ordering)
230 {
231 typedef typename M::global_size_t mat_gs_t;
232 typedef typename Kokkos::View<mat_gs_t*, Kokkos::HostSpace> KV_TMP;
233 size_t i, size = (pointers.extent(0) > 0 ? pointers.extent(0) : 1); // making sure it is at least 1, even for empty local matrix
234 KV_TMP pointers_tmp(Kokkos::ViewAllocateWithoutInitializing("pointers_tmp"), size);
235
236 mat_gs_t nnz_tmp = 0;
237 Op::template apply_kokkos_view<KV_S, KV_GO, KV_TMP>(mat, nzvals,
238 indices, pointers_tmp, nnz_tmp, Teuchos::ptrInArg(*map), distribution, ordering);
239 nnz = Teuchos::as<typename KV_GS::value_type>(nnz_tmp);
240
241 typedef typename KV_GS::value_type view_gs_t;
242 if (pointers.extent(0) == 1) {
243 Kokkos::deep_copy(pointers, 0);
244 } else {
245 auto host_pointers = Kokkos::create_mirror_view(pointers);
246 for (i = 0; i < pointers.extent(0); ++i){
247 host_pointers(i) = Teuchos::as<view_gs_t>(pointers_tmp(i));
248 }
249 Kokkos::deep_copy(pointers, host_pointers);
250 }
251 nnz = Teuchos::as<view_gs_t>(nnz_tmp);
252 }
253 };
254
255 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
256 struct same_go_helper_kokkos_view
257 {
258 static void do_get(const Teuchos::Ptr<const M> mat,
259 KV_S& nzvals,
260 KV_GO& indices,
261 KV_GS& pointers,
262 typename KV_GS::value_type& nnz,
263 const Teuchos::Ptr<
264 const Tpetra::Map<typename M::local_ordinal_t,
265 typename M::global_ordinal_t,
266 typename M::node_t> > map,
267 EDistribution distribution,
268 EStorage_Ordering ordering)
269 {
270 typedef typename M::global_size_t mat_gs_t;
271 typedef typename KV_GS::value_type view_gs_t;
272 std::conditional_t<std::is_same_v<view_gs_t,mat_gs_t>,
273 same_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
274 diff_gs_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::do_get(mat, nzvals, indices,
275 pointers, nnz, map,
276 distribution, ordering);
277 }
278 };
279
280 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
281 struct diff_go_helper_kokkos_view
282 {
283 static void do_get(const Teuchos::Ptr<const M> mat,
284 KV_S& nzvals,
285 KV_GO& indices,
286 KV_GS& pointers,
287 typename KV_GS::value_type& nnz,
288 const Teuchos::Ptr<
289 const Tpetra::Map<typename M::local_ordinal_t,
290 typename M::global_ordinal_t,
291 typename M::node_t> > map,
292 EDistribution distribution,
293 EStorage_Ordering ordering)
294 {
295 typedef typename M::global_ordinal_t mat_go_t;
296 typedef typename M::global_size_t mat_gs_t;
297 typedef typename Kokkos::View<mat_go_t*, Kokkos::HostSpace> KV_TMP;
298 size_t i, size = indices.extent(0);
299 KV_TMP indices_tmp(Kokkos::ViewAllocateWithoutInitializing("indices_tmp"), size);
300
301 typedef typename KV_GO::value_type view_go_t;
302 typedef typename KV_GS::value_type view_gs_t;
303 std::conditional_t<std::is_same_v<view_gs_t,mat_gs_t>,
304 same_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op>,
305 diff_gs_helper_kokkos_view<M, KV_S, KV_TMP, KV_GS, Op> >::do_get(mat, nzvals, indices_tmp,
306 pointers, nnz, map,
307 distribution, ordering);
308 auto host_indices = Kokkos::create_mirror_view(indices);
309 for (i = 0; i < size; ++i){
310 host_indices(i) = Teuchos::as<view_go_t>(indices_tmp(i));
311 }
312 Kokkos::deep_copy(indices, host_indices);
313 }
314 };
315
316 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
317 struct same_scalar_helper_kokkos_view
318 {
319 static void do_get(const Teuchos::Ptr<const M> mat,
320 KV_S& nzvals,
321 KV_GO& indices,
322 KV_GS& pointers,
323 typename KV_GS::value_type& nnz,
324 const Teuchos::Ptr<
325 const Tpetra::Map<typename M::local_ordinal_t,
326 typename M::global_ordinal_t,
327 typename M::node_t> > map,
328 EDistribution distribution,
329 EStorage_Ordering ordering)
330 {
331 typedef typename M::global_ordinal_t mat_go_t;
332 typedef typename KV_GO::value_type view_go_t;
333 std::conditional_t<std::is_same_v<view_go_t, mat_go_t>,
334 same_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op>,
335 diff_go_helper_kokkos_view<M, KV_S, KV_GO, KV_GS, Op> >::do_get(mat, nzvals, indices,
336 pointers, nnz, map,
337 distribution, ordering);
338 }
339 };
340
341 template<class M, typename KV_S, typename KV_GO, typename KV_GS, class Op>
342 struct diff_scalar_helper_kokkos_view
343 {
344 static void do_get(const Teuchos::Ptr<const M> mat,
345 KV_S& nzvals,
346 KV_GO& indices,
347 KV_GS& pointers,
348 typename KV_GS::value_type& nnz,
349 const Teuchos::Ptr<
350 const Tpetra::Map<typename M::local_ordinal_t,
351 typename M::global_ordinal_t,
352 typename M::node_t> > map,
353 EDistribution distribution,
354 EStorage_Ordering ordering)
355 {
356 typedef typename M::global_ordinal_t mat_go_t;
357#if KOKKOS_VERSION >= 40799
358 typedef typename KokkosKernels::ArithTraits<typename M::scalar_t>::val_type mat_scalar_t;
359#else
360 typedef typename Kokkos::ArithTraits<typename M::scalar_t>::val_type mat_scalar_t;
361#endif
362 typedef typename Kokkos::View<mat_scalar_t*, Kokkos::HostSpace> KV_TMP;
363 size_t i, size = nzvals.extent(0);
364 KV_TMP nzvals_tmp(Kokkos::ViewAllocateWithoutInitializing("nzvals_tmp"), size);
365
366 typedef typename KV_S::value_type view_scalar_t;
367 typedef typename KV_GO::value_type view_go_t;
368 std::conditional_t<std::is_same_v<view_go_t, mat_go_t>,
369 same_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op>,
370 diff_go_helper_kokkos_view<M, KV_TMP, KV_GO, KV_GS, Op> >::do_get(mat, nzvals_tmp, indices,
371 pointers, nnz, map,
372 distribution, ordering);
373
374 auto host_nzvals = Kokkos::create_mirror_view(nzvals);
375 for (i = 0; i < size; ++i){
376 host_nzvals(i) = Teuchos::as<view_scalar_t>(nzvals_tmp(i));
377 }
378 Kokkos::deep_copy(nzvals, host_nzvals);
379 }
380 };
381
382
383 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS, class Op>
384 struct get_cxs_helper_kokkos_view
385 {
386 static void do_get(const Teuchos::Ptr<const Matrix> mat,
387 KV_S& nzvals,
388 KV_GO& indices,
389 KV_GS& pointers,
390 typename KV_GS::value_type& nnz,
391 EDistribution distribution,
392 EStorage_Ordering ordering=ARBITRARY,
393 typename KV_GO::value_type indexBase = 0)
394 {
395 typedef typename Matrix::local_ordinal_t lo_t;
396 typedef typename Matrix::global_ordinal_t go_t;
397 typedef typename Matrix::global_size_t gs_t;
398 typedef typename Matrix::node_t node_t;
399
400 const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
401 = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
402 Op::get_dimension(mat),
403 mat->getComm(),
404 indexBase,
405 Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
406 );
407 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
408 }
409
414 static void do_get(const Teuchos::Ptr<const Matrix> mat,
415 KV_S& nzvals,
416 KV_GO& indices,
417 KV_GS& pointers,
418 typename KV_GS::value_type& nnz,
419 EDistribution distribution, // Does this one need a distribution argument??
420 EStorage_Ordering ordering=ARBITRARY)
421 {
422 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
423 typename Matrix::global_ordinal_t,
424 typename Matrix::node_t> > map
425 = Op::getMap(mat);
426 do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
427 }
428
433 static void do_get(const Teuchos::Ptr<const Matrix> mat,
434 KV_S& nzvals,
435 KV_GO& indices,
436 KV_GS& pointers,
437 typename KV_GS::value_type& nnz,
438 const Teuchos::Ptr<
439 const Tpetra::Map<typename Matrix::local_ordinal_t,
440 typename Matrix::global_ordinal_t,
441 typename Matrix::node_t> > map,
442 EDistribution distribution,
443 EStorage_Ordering ordering=ARBITRARY)
444 {
445 typedef typename Matrix::scalar_t mat_scalar;
446 typedef typename KV_S::value_type view_scalar_t;
447
448 std::conditional_t<std::is_same_v<mat_scalar,view_scalar_t>,
449 same_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op>,
450 diff_scalar_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,Op> >::do_get(mat,
451 nzvals, indices,
452 pointers, nnz,
453 map,
454 distribution, ordering);
455 }
456 };
457
458#ifndef DOXYGEN_SHOULD_SKIP_THIS
459 /*
460 * These two function-like classes are meant to be used as the \c
461 * Op template parameter for the \c get_cxs_helper template class.
462 */
463 template<class Matrix>
464 struct get_ccs_func
465 {
466 template<typename KV_S, typename KV_GO, typename KV_GS>
467 static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
468 KV_S& nzvals,
469 KV_GO& rowind,
470 KV_GS& colptr,
471 typename Matrix::global_size_t& nnz,
472 const Teuchos::Ptr<
473 const Tpetra::Map<typename Matrix::local_ordinal_t,
474 typename Matrix::global_ordinal_t,
475 typename Matrix::node_t> > map,
476 EDistribution distribution,
477 EStorage_Ordering ordering)
478 {
479 mat->getCcs_kokkos_view(nzvals, rowind, colptr, nnz, map, ordering, distribution);
480 }
481
482 static
483 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
484 typename Matrix::global_ordinal_t,
485 typename Matrix::node_t> >
486 getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
487 {
488 return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
489 }
490
491 static
492 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
493 typename Matrix::global_ordinal_t,
494 typename Matrix::node_t> >
495 getMap(const Teuchos::Ptr<const Matrix> mat)
496 {
497 return mat->getColMap();
498 }
499
500 static
501 typename Matrix::global_size_t
502 get_dimension(const Teuchos::Ptr<const Matrix> mat)
503 {
504 return mat->getGlobalNumCols();
505 }
506 };
507
508 template<class Matrix>
509 struct get_crs_func
510 {
511 template<typename KV_S, typename KV_GO, typename KV_GS>
512 static void apply_kokkos_view(const Teuchos::Ptr<const Matrix> mat,
513 KV_S& nzvals,
514 KV_GO& colind,
515 KV_GS& rowptr,
516 typename Matrix::global_size_t& nnz,
517 const Teuchos::Ptr<
518 const Tpetra::Map<typename Matrix::local_ordinal_t,
519 typename Matrix::global_ordinal_t,
520 typename Matrix::node_t> > map,
521 EDistribution distribution,
522 EStorage_Ordering ordering)
523 {
524 mat->getCrs_kokkos_view(nzvals, colind, rowptr, nnz, map, ordering, distribution);
525 }
526
527 static
528 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
529 typename Matrix::global_ordinal_t,
530 typename Matrix::node_t> >
531 getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
532 {
533 return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
534 }
535
536 static
537 const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
538 typename Matrix::global_ordinal_t,
539 typename Matrix::node_t> >
540 getMap(const Teuchos::Ptr<const Matrix> mat)
541 {
542 return mat->getRowMap();
543 }
544
545 static
546 typename Matrix::global_size_t
547 get_dimension(const Teuchos::Ptr<const Matrix> mat)
548 {
549 return mat->getGlobalNumRows();
550 }
551 };
552#endif // DOXYGEN_SHOULD_SKIP_THIS
553
591 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
592 struct get_ccs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_ccs_func<Matrix> >
593 {};
594
602 template<class Matrix, typename KV_S, typename KV_GO, typename KV_GS>
603 struct get_crs_helper_kokkos_view : get_cxs_helper_kokkos_view<Matrix,KV_S,KV_GO,KV_GS,get_crs_func<Matrix> >
604 {};
605 /* End Matrix/MultiVector Utilities */
606
607
609 // Definitions //
611
612
613 template <typename LO, typename GO, typename GS, typename Node>
614 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
615 getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
616 {
617 //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
618 Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
619 return gather_map;
620 }
621
622
623 template <typename LO, typename GO, typename GS, typename Node>
624 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
625 getDistributionMap(EDistribution distribution,
626 GS num_global_elements,
627 const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
628 GO indexBase,
629 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
630 {
631 // TODO: Need to add indexBase to cases other than ROOTED
632 // We do not support these maps in any solver now.
633 switch( distribution ){
634 case DISTRIBUTED:
635 case DISTRIBUTED_NO_OVERLAP:
636 return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
637 case GLOBALLY_REPLICATED:
638 return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
639 case ROOTED:
640 {
641 int rank = Teuchos::rank(*comm);
642 size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
643 if( rank == 0 ) { my_num_elems = num_global_elements; }
644
645 return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
646 my_num_elems, indexBase, comm));
647 }
648 case CONTIGUOUS_AND_ROOTED:
649 {
650 const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
651 = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
652 return gathermap;
653 }
654 default:
655 TEUCHOS_TEST_FOR_EXCEPTION( true,
656 std::logic_error,
657 "Control should never reach this point. "
658 "Please contact the Amesos2 developers." );
659 }
660 }
661
662 template <typename Scalar,
663 typename GlobalOrdinal,
664 typename GlobalSizeT>
665 void transpose(Teuchos::ArrayView<Scalar> vals,
666 Teuchos::ArrayView<GlobalOrdinal> indices,
667 Teuchos::ArrayView<GlobalSizeT> ptr,
668 Teuchos::ArrayView<Scalar> trans_vals,
669 Teuchos::ArrayView<GlobalOrdinal> trans_indices,
670 Teuchos::ArrayView<GlobalSizeT> trans_ptr)
671 {
672 /* We have a compressed-row storage format of this matrix. We
673 * transform this into a compressed-column format using a
674 * distribution-counting sort algorithm, which is described by
675 * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
676 */
677
678#ifdef HAVE_AMESOS2_DEBUG
679 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
680 ind_begin = indices.begin();
681 ind_end = indices.end();
682 size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
683 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
684 std::invalid_argument,
685 "Transpose pointer size not large enough." );
686 TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
687 std::invalid_argument,
688 "Transpose values array not large enough." );
689 TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
690 std::invalid_argument,
691 "Transpose indices array not large enough." );
692#else
693 typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
694#endif
695 // Count the number of entries in each column
696 Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
697 ind_end = indices.end();
698 for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
699 ++(count[(*ind_it) + 1]);
700 }
701 // Accumulate
702 typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
703 cnt_end = count.end();
704 for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
705 *cnt_it = *cnt_it + *(cnt_it - 1);
706 }
707 // This becomes the array of column pointers
708 trans_ptr.assign(count);
709
710 /* Move the nonzero values into their final place in nzval, based on the
711 * counts found previously.
712 *
713 * This sequence deviates from Knuth's algorithm a bit, following more
714 * closely the description presented in Gustavson, Fred G. "Two Fast
715 * Algorithms for Sparse Matrices: Multiplication and Permuted
716 * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
717 * 250--269, http://doi.acm.org/10.1145/355791.355796.
718 *
719 * The output indices end up in sorted order
720 */
721
722 GlobalSizeT size = ptr.size();
723 for( GlobalSizeT i = 0; i < size - 1; ++i ){
724 GlobalOrdinal u = ptr[i];
725 GlobalOrdinal v = ptr[i + 1];
726 for( GlobalOrdinal j = u; j < v; ++j ){
727 GlobalOrdinal k = count[indices[j]];
728 trans_vals[k] = vals[j];
729 trans_indices[k] = i;
730 ++(count[indices[j]]);
731 }
732 }
733 }
734
735
736 template <typename Scalar1, typename Scalar2>
737 void
738 scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
739 size_t ld, Teuchos::ArrayView<Scalar2> s)
740 {
741 size_t vals_size = vals.size();
742#ifdef HAVE_AMESOS2_DEBUG
743 size_t s_size = s.size();
744 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
745 std::invalid_argument,
746 "Scale vector must have length at least that of the vector" );
747#endif
748 size_t i, s_i;
749 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
750 if( s_i == l ){
751 // bring i to the next multiple of ld
752 i += ld - s_i;
753 s_i = 0;
754 }
755 vals[i] *= s[s_i];
756 }
757 }
758
759 template <typename Scalar1, typename Scalar2, class BinaryOp>
760 void
761 scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
762 size_t ld, Teuchos::ArrayView<Scalar2> s,
763 BinaryOp binary_op)
764 {
765 size_t vals_size = vals.size();
766#ifdef HAVE_AMESOS2_DEBUG
767 size_t s_size = s.size();
768 TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
769 std::invalid_argument,
770 "Scale vector must have length at least that of the vector" );
771#endif
772 size_t i, s_i;
773 for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
774 if( s_i == l ){
775 // bring i to the next multiple of ld
776 i += ld - s_i;
777 s_i = 0;
778 }
779 vals[i] = binary_op(vals[i], s[s_i]);
780 }
781 }
782
783 template<class row_ptr_view_t, class cols_view_t, class per_view_t>
784 void
785 reorder(row_ptr_view_t & row_ptr, cols_view_t & cols,
786 per_view_t & perm, per_view_t & peri, size_t & nnz,
787 bool permute_matrix)
788 {
789 #ifndef HAVE_AMESOS2_METIS
790 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
791 "Cannot reorder for cuSolver because no METIS is available.");
792 #else
793 typedef typename cols_view_t::value_type ordinal_type;
794 typedef typename row_ptr_view_t::value_type size_type;
795
796 // begin on host where we'll run metis reorder
797 auto host_row_ptr = Kokkos::create_mirror_view(row_ptr);
798 auto host_cols = Kokkos::create_mirror_view(cols);
799 Kokkos::deep_copy(host_row_ptr, row_ptr);
800 Kokkos::deep_copy(host_cols, cols);
801
802 // strip out the diagonals - metis will just crash with them included.
803 // make space for the stripped version
804 typedef Kokkos::View<idx_t*, Kokkos::HostSpace> host_metis_array;
805 const ordinal_type size = row_ptr.size() - 1;
806 size_type max_nnz = host_row_ptr(size);
807 host_metis_array host_strip_diag_row_ptr(
808 Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_row_ptr"),
809 size+1);
810 host_metis_array host_strip_diag_cols(
811 Kokkos::ViewAllocateWithoutInitializing("host_strip_diag_cols"),
812 max_nnz);
813
814 size_type new_nnz = 0;
815 for(ordinal_type i = 0; i < size; ++i) {
816 host_strip_diag_row_ptr(i) = new_nnz;
817 for(size_type j = host_row_ptr(i); j < host_row_ptr(i+1); ++j) {
818 if (i != host_cols(j)) {
819 host_strip_diag_cols(new_nnz++) = host_cols(j);
820 }
821 }
822 }
823 host_strip_diag_row_ptr(size) = new_nnz;
824
825 // we'll get original permutations on host
826 host_metis_array host_perm(
827 Kokkos::ViewAllocateWithoutInitializing("host_perm"), size);
828 host_metis_array host_peri(
829 Kokkos::ViewAllocateWithoutInitializing("host_peri"), size);
830
831 // If we want to remove metis.h included in this header we can move this
832 // to the cpp, but we need to decide how to handle the idx_t declaration.
833 idx_t metis_size = size;
834 int err = METIS_NodeND(&metis_size, host_strip_diag_row_ptr.data(), host_strip_diag_cols.data(),
835 NULL, NULL, host_perm.data(), host_peri.data());
836
837 TEUCHOS_TEST_FOR_EXCEPTION(err != METIS_OK, std::runtime_error,
838 "METIS_NodeND failed to sort matrix.");
839
840 // put the permutations on our saved device ptrs
841 // these will be used to permute x and b when we solve
842 typedef typename cols_view_t::execution_space exec_space_t;
843 auto device_perm = Kokkos::create_mirror_view(exec_space_t(), host_perm);
844 auto device_peri = Kokkos::create_mirror_view(exec_space_t(), host_peri);
845 deep_copy(device_perm, host_perm);
846 deep_copy(device_peri, host_peri);
847
848 // also set the permutation which may need to convert the type from
849 // metis to the native ordinal_type
850 deep_copy_or_assign_view(perm, device_perm);
851 deep_copy_or_assign_view(peri, device_peri);
852
853 if (permute_matrix) {
854 // we'll permute matrix on device to a new set of arrays
855 row_ptr_view_t new_row_ptr(
856 Kokkos::ViewAllocateWithoutInitializing("new_row_ptr"), row_ptr.size());
857 cols_view_t new_cols(
858 Kokkos::ViewAllocateWithoutInitializing("new_cols"), cols.size() - new_nnz/2);
859
860 // permute row indices
861 Kokkos::RangePolicy<exec_space_t> policy_row(0, row_ptr.size());
862 Kokkos::parallel_scan(policy_row, KOKKOS_LAMBDA(
863 ordinal_type i, size_type & update, const bool &final) {
864 if(final) {
865 new_row_ptr(i) = update;
866 }
867 if(i < size) {
868 ordinal_type count = 0;
869 const ordinal_type row = device_perm(i);
870 for(ordinal_type k = row_ptr(row); k < row_ptr(row + 1); ++k) {
871 const ordinal_type j = device_peri(cols(k));
872 count += (i >= j);
873 }
874 update += count;
875 }
876 });
877
878 // permute col indices
879 Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
880 Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
881 const ordinal_type kbeg = new_row_ptr(i);
882 const ordinal_type row = device_perm(i);
883 const ordinal_type col_beg = row_ptr(row);
884 const ordinal_type col_end = row_ptr(row + 1);
885 const ordinal_type nk = col_end - col_beg;
886 for(ordinal_type k = 0, t = 0; k < nk; ++k) {
887 const ordinal_type tk = kbeg + t;
888 const ordinal_type sk = col_beg + k;
889 const ordinal_type j = device_peri(cols(sk));
890 if(i >= j) {
891 new_cols(tk) = j;
892 ++t;
893 }
894 }
895 });
896
897 // finally set the inputs to the new sorted arrays
898 row_ptr = new_row_ptr;
899 cols = new_cols;
900 }
901
902 nnz = new_nnz;
903 #endif // HAVE_AMESOS2_METIS
904 }
905
906 template<class values_view_t, class row_ptr_view_t,
907 class cols_view_t, class per_view_t>
908 void
909 reorder_values(values_view_t & values, const row_ptr_view_t & orig_row_ptr,
910 const row_ptr_view_t & new_row_ptr,
911 const cols_view_t & orig_cols, const per_view_t & perm, const per_view_t & peri,
912 size_t nnz)
913 {
914 typedef typename cols_view_t::value_type ordinal_type;
915 typedef typename cols_view_t::execution_space exec_space_t;
916
917 auto device_perm = Kokkos::create_mirror_view(exec_space_t(), perm);
918 auto device_peri = Kokkos::create_mirror_view(exec_space_t(), peri);
919 deep_copy(device_perm, perm);
920 deep_copy(device_peri, peri);
921
922 const ordinal_type size = orig_row_ptr.size() - 1;
923
924 auto host_orig_row_ptr = Kokkos::create_mirror_view(orig_row_ptr);
925 auto new_nnz = host_orig_row_ptr(size); // TODO: Maybe optimize this by caching
926
927 values_view_t new_values(
928 Kokkos::ViewAllocateWithoutInitializing("new_values"), values.size() - new_nnz/2);
929
930 // permute col indices
931 Kokkos::RangePolicy<exec_space_t> policy_col(0, size);
932 Kokkos::parallel_for(policy_col, KOKKOS_LAMBDA(ordinal_type i) {
933 const ordinal_type kbeg = new_row_ptr(i);
934 const ordinal_type row = device_perm(i);
935 const ordinal_type col_beg = orig_row_ptr(row);
936 const ordinal_type col_end = orig_row_ptr(row + 1);
937 const ordinal_type nk = col_end - col_beg;
938 for(ordinal_type k = 0, t = 0; k < nk; ++k) {
939 const ordinal_type tk = kbeg + t;
940 const ordinal_type sk = col_beg + k;
941 const ordinal_type j = device_peri(orig_cols(sk));
942 if(i >= j) {
943 new_values(tk) = values(sk);
944 ++t;
945 }
946 }
947 });
948
949 values = new_values;
950 }
951
952 template<class array_view_t, class per_view_t>
953 void
954 apply_reorder_permutation(const array_view_t & array,
955 array_view_t & permuted_array, const per_view_t & permutation) {
956 if(permuted_array.extent(0) != array.extent(0) || permuted_array.extent(1) != array.extent(1)) {
957 permuted_array = array_view_t(
958 Kokkos::ViewAllocateWithoutInitializing("permuted_array"),
959 array.extent(0), array.extent(1));
960 }
961 typedef typename array_view_t::execution_space exec_space_t;
962 Kokkos::RangePolicy<exec_space_t> policy(0, array.extent(0));
963 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(size_t i) {
964 for(size_t j = 0; j < array.extent(1); ++j) {
965 permuted_array(i, j) = array(permutation(i), j);
966 }
967 });
968 }
969
970
971 // used to read matrix with gapped GIDs for test/example
972 template <typename GO, typename Scalar>
973 bool
974 readEntryFromFile (GO& gblRowInd, GO& gblColInd, Scalar& val, const std::string& s)
975 {
976 if (s.size () == 0 || s.find ("%") != std::string::npos) {
977 return false; // empty line or comment line
978 }
979 std::istringstream in (s);
980
981 if (! in) {
982 return false;
983 }
984 in >> gblRowInd;
985 if (! in) {
986 return false;
987 }
988 in >> gblColInd;
989 if (! in) {
990 return false;
991 }
992 in >> val;
993 return true;
994 }
995
996 template<class map_type, class MAT>
997 Teuchos::RCP<MAT>
998 readCrsMatrixFromFile (const std::string& matrixFilename,
999 Teuchos::RCP<Teuchos::FancyOStream> & fos,
1000 const Teuchos::RCP<const map_type>& rowMap,
1001 const Teuchos::RCP<const map_type>& domainMap,
1002 const Teuchos::RCP<const map_type>& rangeMap,
1003 const bool convert_to_zero_base,
1004 const int header_size)
1005 {
1006 using Scalar = typename MAT::scalar_type;
1007 using GO = typename MAT::global_ordinal_type;
1008
1009 using counter_type = std::map<GO, size_t>;
1010 using pair_type = std::pair<const GO, size_t>;
1011 using Teuchos::RCP;
1012
1013 auto comm = rowMap->getComm ();
1014 const int myRank = comm->getRank ();
1015
1016 std::ifstream inFile;
1017 int opened = 0;
1018 if (myRank == 0)
1019 {
1020 try {
1021 inFile.open (matrixFilename);
1022 if (inFile) {
1023 opened = 1;
1024 }
1025 }
1026 catch (...) {
1027 opened = 0;
1028 }
1029 }
1030 Teuchos::broadcast<int, int> (*comm, 0, Teuchos::outArg (opened));
1031 TEUCHOS_TEST_FOR_EXCEPTION
1032 (opened == 0, std::runtime_error, "readCrsMatrixFromFile: "
1033 "Failed to open file \"" << matrixFilename << "\" on Process 0.");
1034
1035 RCP<MAT> A;
1036 if (myRank == 0)
1037 {
1038 std::string line;
1039
1040 // Skip the first N lines. This is a hack, specific to the input file in question.
1041 //*fos << " Reading matrix market file. Skip " << header_size << " header lines" << std::endl;
1042 for ( int i = 0; i < header_size; ++i ) {
1043 std::getline (inFile, line);
1044 }
1045
1046 counter_type counts;
1047 Teuchos::Array<Scalar> vals;
1048 Teuchos::Array<GO> gblRowInds;
1049 Teuchos::Array<GO> gblColInds;
1050 while (inFile) {
1051 std::getline (inFile, line);
1052 GO gblRowInd {};
1053 GO gblColInd {};
1054 Scalar val {};
1055 const bool gotLine = readEntryFromFile (gblRowInd, gblColInd, val, line);
1056 if (gotLine) {
1057 //*fos << " read mtx rank: " << myRank << " | gblRowInd = " << gblRowInd << " gblColInd = " << gblColInd << std::endl;
1058 if ( convert_to_zero_base ) {
1059 gblRowInd -= 1 ;
1060 gblColInd -= 1 ;
1061 }
1062 counts[gblRowInd]++;
1063 vals.push_back(val);
1064 gblRowInds.push_back(gblRowInd);
1065 gblColInds.push_back(gblColInd);
1066 }
1067 }
1068
1069 // Max number of entries in any row
1070 auto pr = std::max_element(
1071 std::begin(counts),
1072 std::end(counts),
1073 [] (pair_type const& p1, pair_type const& p2){ return p1.second < p2.second; }
1074 );
1075 size_t maxCount = (counts.empty()) ? size_t(0) : pr->second;
1076 A = Teuchos::rcp(new MAT(rowMap, maxCount));
1077 for (typename Teuchos::Array<GO>::size_type i=0; i<gblRowInds.size(); i++) {
1078 A->insertGlobalValues (gblRowInds[i], gblColInds(i,1), vals(i,1));
1079 }
1080 } else {
1081 A = Teuchos::rcp(new MAT(rowMap, 0));
1082 }
1083
1084 A->fillComplete (domainMap, rangeMap);
1085 return A;
1086 }
1088
1089 } // end namespace Util
1090
1091} // end namespace Amesos2
1092
1093#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:90
Amesos2::GLOBALLY_REPLICATED
@ GLOBALLY_REPLICATED
Definition Amesos2_TypeDecl.hpp:92
Amesos2::DISTRIBUTED_NO_OVERLAP
@ DISTRIBUTED_NO_OVERLAP
Definition Amesos2_TypeDecl.hpp:91
Amesos2::ROOTED
@ ROOTED
Definition Amesos2_TypeDecl.hpp:93
Amesos2::CONTIGUOUS_AND_ROOTED
@ CONTIGUOUS_AND_ROOTED
Definition Amesos2_TypeDecl.hpp:94
Amesos2::EDistribution
EDistribution
Definition Amesos2_TypeDecl.hpp:89
Amesos2::EStorage_Ordering
EStorage_Ordering
Definition Amesos2_TypeDecl.hpp:107
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::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:615
size
const int size
Definition klu2_simple.cpp:50
Amesos2::Util::get_ccs_helper_kokkos_view
A generic helper class for getting a CCS representation of a Matrix.
Definition Amesos2_Util.hpp:593
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:604
Amesos2::Util::same_gs_helper_kokkos_view
Scales a 1-D representation of a multivector.
Definition Amesos2_Util.hpp:198
Generated by doxygen 1.15.0