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_ArrayView.hpp>
60 #include <Teuchos_FancyOStream.hpp>
61 
62 #include <Tpetra_Map.hpp>
63 
64 #include "Amesos2_TypeDecl.hpp"
65 #include "Amesos2_Meta.hpp"
66 
67 #ifdef HAVE_TPETRA_INST_INT_INT
68 #ifdef HAVE_AMESOS2_EPETRA
69 #include <Epetra_Map.h>
70 #endif
71 #endif
72 
73 
74 namespace Amesos2 {
75 
76  namespace Util {
77 
84  using Teuchos::RCP;
85  using Teuchos::ArrayView;
86 
87  using Meta::is_same;
88  using Meta::if_then_else;
89 
105  template <typename LO, typename GO, typename GS, typename Node>
106  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
107  getDistributionMap(EDistribution distribution,
108  GS num_global_elements,
109  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
110  GO indexBase = 0);
111 
112 
113 #ifdef HAVE_TPETRA_INST_INT_INT
114 #ifdef HAVE_AMESOS2_EPETRA
115 
121  template <typename LO, typename GO, typename GS, typename Node>
122  RCP<Tpetra::Map<LO,GO,Node> >
123  epetra_map_to_tpetra_map(const Epetra_BlockMap& map);
124 
130  template <typename LO, typename GO, typename GS, typename Node>
131  RCP<Epetra_Map>
132  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map);
133 
139  const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
140 
146  const RCP<const Epetra_Comm> to_epetra_comm(RCP<const Teuchos::Comm<int> > c);
147 #endif // HAVE_AMESOS2_EPETRA
148 #endif // HAVE_TPETRA_INST_INT_INT
149 
155  template <typename Scalar,
156  typename GlobalOrdinal,
157  typename GlobalSizeT>
158  void transpose(ArrayView<Scalar> vals,
159  ArrayView<GlobalOrdinal> indices,
160  ArrayView<GlobalSizeT> ptr,
161  ArrayView<Scalar> trans_vals,
162  ArrayView<GlobalOrdinal> trans_indices,
163  ArrayView<GlobalSizeT> trans_ptr);
164 
178  template <typename Scalar1, typename Scalar2>
179  void scale(ArrayView<Scalar1> vals, size_t l,
180  size_t ld, ArrayView<Scalar2> s);
181 
200  template <typename Scalar1, typename Scalar2, class BinaryOp>
201  void scale(ArrayView<Scalar1> vals, size_t l,
202  size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
203 
204 
206  void printLine( Teuchos::FancyOStream &out );
207 
208 
210  // Matrix/MultiVector Utilities //
212 
213 #ifndef DOXYGEN_SHOULD_SKIP_THIS
214  /*
215  * The following represents a general way of getting a CRS or CCS
216  * representation of a matrix with implicit type conversions. The
217  * \c get_crs_helper and \c get_ccs_helper classes are templated
218  * on 4 types:
219  *
220  * - A Matrix type (conforming to the Amesos2 MatrixAdapter interface)
221  * - A scalar type
222  * - A global ordinal type
223  * - A global size type
224  *
225  * The last three template types correspond to the input argument
226  * types. For example, if the scalar type is \c double , then we
227  * require that the \c nzvals argument is a \c
228  * Teuchos::ArrayView<double> type.
229  *
230  * These helpers perform any type conversions that must be
231  * performed to go between the Matrix's types and the input types.
232  * If no conversions are necessary the static functions can be
233  * effectively inlined.
234  */
235 
236  template <class M, typename S, typename GO, typename GS, class Op>
237  struct same_gs_helper
238  {
239  static void do_get(const Teuchos::Ptr<const M> mat,
240  const ArrayView<typename M::scalar_t> nzvals,
241  const ArrayView<typename M::global_ordinal_t> indices,
242  const ArrayView<GS> pointers,
243  GS& nnz,
244  const Teuchos::Ptr<
245  const Tpetra::Map<typename M::local_ordinal_t,
246  typename M::global_ordinal_t,
247  typename M::node_t> > map,
248  EStorage_Ordering ordering)
249  {
250  Op::apply(mat, nzvals, indices, pointers, nnz, map, ordering);
251  }
252  };
253 
254  template <class M, typename S, typename GO, typename GS, class Op>
255  struct diff_gs_helper
256  {
257  static void do_get(const Teuchos::Ptr<const M> mat,
258  const ArrayView<typename M::scalar_t> nzvals,
259  const ArrayView<typename M::global_ordinal_t> indices,
260  const ArrayView<GS> pointers,
261  GS& nnz,
262  const Teuchos::Ptr<
263  const Tpetra::Map<typename M::local_ordinal_t,
264  typename M::global_ordinal_t,
265  typename M::node_t> > map,
266  EStorage_Ordering ordering)
267  {
268  typedef typename M::global_size_t mat_gs_t;
269  typename ArrayView<GS>::size_type i, size = pointers.size();
270  Teuchos::Array<mat_gs_t> pointers_tmp(size);
271  mat_gs_t nnz_tmp = 0;
272 
273  Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, ordering);
274 
275  for (i = 0; i < size; ++i){
276  pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
277  }
278  nnz = Teuchos::as<GS>(nnz_tmp);
279  }
280  };
281 
282  template <class M, typename S, typename GO, typename GS, class Op>
283  struct same_go_helper
284  {
285  static void do_get(const Teuchos::Ptr<const M> mat,
286  const ArrayView<typename M::scalar_t> nzvals,
287  const ArrayView<GO> indices,
288  const ArrayView<GS> pointers,
289  GS& nnz,
290  const Teuchos::Ptr<
291  const Tpetra::Map<typename M::local_ordinal_t,
292  typename M::global_ordinal_t,
293  typename M::node_t> > map,
294  EStorage_Ordering ordering)
295  {
296  typedef typename M::global_size_t mat_gs_t;
297  if_then_else<is_same<GS,mat_gs_t>::value,
298  same_gs_helper<M,S,GO,GS,Op>,
299  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
300  pointers, nnz, map,
301  ordering);
302  }
303  };
304 
305  template <class M, typename S, typename GO, typename GS, class Op>
306  struct diff_go_helper
307  {
308  static void do_get(const Teuchos::Ptr<const M> mat,
309  const ArrayView<typename M::scalar_t> nzvals,
310  const ArrayView<GO> indices,
311  const ArrayView<GS> pointers,
312  GS& nnz,
313  const Teuchos::Ptr<
314  const Tpetra::Map<typename M::local_ordinal_t,
315  typename M::global_ordinal_t,
316  typename M::node_t> > map,
317  EStorage_Ordering ordering)
318  {
319  typedef typename M::global_ordinal_t mat_go_t;
320  typedef typename M::global_size_t mat_gs_t;
321  typename ArrayView<GO>::size_type i, size = indices.size();
322  Teuchos::Array<mat_go_t> indices_tmp(size);
323 
324  if_then_else<is_same<GS,mat_gs_t>::value,
325  same_gs_helper<M,S,GO,GS,Op>,
326  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
327  pointers, nnz, map,
328  ordering);
329 
330  for (i = 0; i < size; ++i){
331  indices[i] = Teuchos::as<GO>(indices_tmp[i]);
332  }
333  }
334  };
335 
336  template <class M, typename S, typename GO, typename GS, class Op>
337  struct same_scalar_helper
338  {
339  static void do_get(const Teuchos::Ptr<const M> mat,
340  const ArrayView<S> nzvals,
341  const ArrayView<GO> indices,
342  const ArrayView<GS> pointers,
343  GS& nnz,
344  const Teuchos::Ptr<
345  const Tpetra::Map<typename M::local_ordinal_t,
346  typename M::global_ordinal_t,
347  typename M::node_t> > map,
348  EStorage_Ordering ordering)
349  {
350  typedef typename M::global_ordinal_t mat_go_t;
351  if_then_else<is_same<GO,mat_go_t>::value,
352  same_go_helper<M,S,GO,GS,Op>,
353  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
354  pointers, nnz, map,
355  ordering);
356  }
357  };
358 
359  template <class M, typename S, typename GO, typename GS, class Op>
360  struct diff_scalar_helper
361  {
362  static void do_get(const Teuchos::Ptr<const M> mat,
363  const ArrayView<S> nzvals,
364  const ArrayView<GO> indices,
365  const ArrayView<GS> pointers,
366  GS& nnz,
367  const Teuchos::Ptr<
368  const Tpetra::Map<typename M::local_ordinal_t,
369  typename M::global_ordinal_t,
370  typename M::node_t> > map,
371  EStorage_Ordering ordering)
372  {
373  typedef typename M::scalar_t mat_scalar_t;
374  typedef typename M::global_ordinal_t mat_go_t;
375  typename ArrayView<S>::size_type i, size = nzvals.size();
376  Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
377 
378  if_then_else<is_same<GO,mat_go_t>::value,
379  same_go_helper<M,S,GO,GS,Op>,
380  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
381  pointers, nnz, map,
382  ordering);
383 
384  for (i = 0; i < size; ++i){
385  nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
386  }
387  }
388  };
389 #endif // DOXYGEN_SHOULD_SKIP_THIS
390 
403  template<class Matrix, typename S, typename GO, typename GS, class Op>
405  {
406  static void do_get(const Teuchos::Ptr<const Matrix> mat,
407  const ArrayView<S> nzvals,
408  const ArrayView<GO> indices,
409  const ArrayView<GS> pointers,
410  GS& nnz, EDistribution distribution,
411  EStorage_Ordering ordering=ARBITRARY,
412  GO indexBase = 0)
413  {
414  typedef typename Matrix::local_ordinal_t lo_t;
415  typedef typename Matrix::global_ordinal_t go_t;
416  typedef typename Matrix::global_size_t gs_t;
417  typedef typename Matrix::node_t node_t;
418 
419  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
420  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
421  Op::get_dimension(mat),
422  mat->getComm(), indexBase);
423  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
424  }
425 
430  static void do_get(const Teuchos::Ptr<const Matrix> mat,
431  const ArrayView<S> nzvals,
432  const ArrayView<GO> indices,
433  const ArrayView<GS> pointers,
434  GS& nnz, EStorage_Ordering ordering=ARBITRARY)
435  {
436  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
437  typename Matrix::global_ordinal_t,
438  typename Matrix::node_t> > map
439  = Op::getMap(mat);
440  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), ordering);
441  }
442 
447  static void do_get(const Teuchos::Ptr<const Matrix> mat,
448  const ArrayView<S> nzvals,
449  const ArrayView<GO> indices,
450  const ArrayView<GS> pointers,
451  GS& nnz,
452  const Teuchos::Ptr<
453  const Tpetra::Map<typename Matrix::local_ordinal_t,
454  typename Matrix::global_ordinal_t,
455  typename Matrix::node_t> > map,
456  EStorage_Ordering ordering=ARBITRARY)
457  {
458  typedef typename Matrix::scalar_t mat_scalar;
459  if_then_else<is_same<mat_scalar,S>::value,
460  same_scalar_helper<Matrix,S,GO,GS,Op>,
461  diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
462  nzvals, indices,
463  pointers, nnz,
464  map, ordering);
465  }
466  };
467 
468 #ifndef DOXYGEN_SHOULD_SKIP_THIS
469  /*
470  * These two function-like classes are meant to be used as the \c
471  * Op template parameter for the \c get_cxs_helper template class.
472  */
473  template<class Matrix>
474  struct get_ccs_func
475  {
476  static void apply(const Teuchos::Ptr<const Matrix> mat,
477  const ArrayView<typename Matrix::scalar_t> nzvals,
478  const ArrayView<typename Matrix::global_ordinal_t> rowind,
479  const ArrayView<typename Matrix::global_size_t> colptr,
480  typename Matrix::global_size_t& nnz,
481  const Teuchos::Ptr<
482  const Tpetra::Map<typename Matrix::local_ordinal_t,
483  typename Matrix::global_ordinal_t,
484  typename Matrix::node_t> > map,
485  EStorage_Ordering ordering)
486  {
487  mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
488  }
489 
490  static
491  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
492  typename Matrix::global_ordinal_t,
493  typename Matrix::node_t> >
494  getMap(const Teuchos::Ptr<const Matrix> mat)
495  {
496  return mat->getColMap();
497  }
498 
499  static
500  typename Matrix::global_size_t
501  get_dimension(const Teuchos::Ptr<const Matrix> mat)
502  {
503  return mat->getGlobalNumCols();
504  }
505  };
506 
507  template<class Matrix>
508  struct get_crs_func
509  {
510  static void apply(const Teuchos::Ptr<const Matrix> mat,
511  const ArrayView<typename Matrix::scalar_t> nzvals,
512  const ArrayView<typename Matrix::global_ordinal_t> colind,
513  const ArrayView<typename Matrix::global_size_t> rowptr,
514  typename Matrix::global_size_t& nnz,
515  const Teuchos::Ptr<
516  const Tpetra::Map<typename Matrix::local_ordinal_t,
517  typename Matrix::global_ordinal_t,
518  typename Matrix::node_t> > map,
519  EStorage_Ordering ordering)
520  {
521  mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering);
522  }
523 
524  static
525  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
526  typename Matrix::global_ordinal_t,
527  typename Matrix::node_t> >
528  getMap(const Teuchos::Ptr<const Matrix> mat)
529  {
530  return mat->getRowMap();
531  }
532 
533  static
534  typename Matrix::global_size_t
535  get_dimension(const Teuchos::Ptr<const Matrix> mat)
536  {
537  return mat->getGlobalNumRows();
538  }
539  };
540 #endif // DOXYGEN_SHOULD_SKIP_THIS
541 
579  template<class Matrix, typename S, typename GO, typename GS>
580  struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
581  {};
582 
590  template<class Matrix, typename S, typename GO, typename GS>
591  struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
592  {};
593 
594  /* End Matrix/MultiVector Utilities */
595 
596 
598  // Definitions //
600 
601  template <typename LO, typename GO, typename GS, typename Node>
602  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
604  GS num_global_elements,
605  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
606  GO indexBase)
607  {
608  // TODO: Need to add indexBase to cases other than ROOTED
609  // We do not support these maps in any solver now.
610  switch( distribution ){
611  case DISTRIBUTED:
613  return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
614  break;
615  case GLOBALLY_REPLICATED:
616  return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
617  break;
618  case ROOTED:
619  {
620  int rank = Teuchos::rank(*comm);
621  size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
622  if( rank == 0 ) my_num_elems = num_global_elements;
623  return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
624  my_num_elems, indexBase, comm));
625  break;
626  }
627  default:
628  TEUCHOS_TEST_FOR_EXCEPTION( true,
629  std::logic_error,
630  "Control should never reach this point. "
631  "Please contact the Amesos2 developers." );
632  break;
633  }
634  }
635 
636 
637 #ifdef HAVE_TPETRA_INST_INT_INT
638 #ifdef HAVE_AMESOS2_EPETRA
639 
640  //#pragma message "include 3"
641  //#include <Epetra_Map.h>
642 
643  template <typename LO, typename GO, typename GS, typename Node>
644  Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
645  epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
646  {
647  using Teuchos::as;
648  using Teuchos::rcp;
649 
650  int num_my_elements = map.NumMyElements();
651  Teuchos::Array<int> my_global_elements(num_my_elements);
652  map.MyGlobalElements(my_global_elements.getRawPtr());
653 
654  typedef Tpetra::Map<LO,GO,Node> map_t;
655  RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
656  my_global_elements(),
657  as<GO>(map.IndexBase()),
658  to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
659  return tmap;
660  }
661 
662  template <typename LO, typename GO, typename GS, typename Node>
663  Teuchos::RCP<Epetra_Map>
664  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
665  {
666  using Teuchos::as;
667 
668  Teuchos::Array<GO> elements_tmp;
669  elements_tmp = map.getNodeElementList();
670  int num_my_elements = elements_tmp.size();
671  Teuchos::Array<int> my_global_elements(num_my_elements);
672  for (int i = 0; i < num_my_elements; ++i){
673  my_global_elements[i] = as<int>(elements_tmp[i]);
674  }
675 
676  using Teuchos::rcp;
677  RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
678  num_my_elements,
679  my_global_elements.getRawPtr(),
680  as<GO>(map.getIndexBase()),
681  *to_epetra_comm(map.getComm())));
682  return emap;
683  }
684 #endif // HAVE_AMESOS2_EPETRA
685 #endif // HAVE_TPETRA_INST_INT_INT
686 
687  template <typename Scalar,
688  typename GlobalOrdinal,
689  typename GlobalSizeT>
690  void transpose(Teuchos::ArrayView<Scalar> vals,
691  Teuchos::ArrayView<GlobalOrdinal> indices,
692  Teuchos::ArrayView<GlobalSizeT> ptr,
693  Teuchos::ArrayView<Scalar> trans_vals,
694  Teuchos::ArrayView<GlobalOrdinal> trans_indices,
695  Teuchos::ArrayView<GlobalSizeT> trans_ptr)
696  {
697  /* We have a compressed-row storage format of this matrix. We
698  * transform this into a compressed-column format using a
699  * distribution-counting sort algorithm, which is described by
700  * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
701  */
702 
703 #ifdef HAVE_AMESOS2_DEBUG
704  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
705  ind_begin = indices.begin();
706  ind_end = indices.end();
707  size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
708  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
709  std::invalid_argument,
710  "Transpose pointer size not large enough." );
711  TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
712  std::invalid_argument,
713  "Transpose values array not large enough." );
714  TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
715  std::invalid_argument,
716  "Transpose indices array not large enough." );
717 #else
718  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
719 #endif
720 
721  // Count the number of entries in each column
722  Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
723  ind_end = indices.end();
724  for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
725  ++(count[(*ind_it) + 1]);
726  }
727 
728  // Accumulate
729  typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
730  cnt_end = count.end();
731  for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
732  *cnt_it = *cnt_it + *(cnt_it - 1);
733  }
734  // This becomes the array of column pointers
735  trans_ptr.assign(count);
736 
737  /* Move the nonzero values into their final place in nzval, based on the
738  * counts found previously.
739  *
740  * This sequence deviates from Knuth's algorithm a bit, following more
741  * closely the description presented in Gustavson, Fred G. "Two Fast
742  * Algorithms for Sparse Matrices: Multiplication and Permuted
743  * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
744  * 250--269, http://doi.acm.org/10.1145/355791.355796.
745  *
746  * The output indices end up in sorted order
747  */
748 
749  GlobalSizeT size = ptr.size();
750  for( GlobalSizeT i = 0; i < size - 1; ++i ){
751  GlobalOrdinal u = ptr[i];
752  GlobalOrdinal v = ptr[i + 1];
753  for( GlobalOrdinal j = u; j < v; ++j ){
754  GlobalOrdinal k = count[indices[j]];
755  trans_vals[k] = vals[j];
756  trans_indices[k] = i;
757  ++(count[indices[j]]);
758  }
759  }
760  }
761 
762 
763  template <typename Scalar1, typename Scalar2>
764  void
765  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
766  size_t ld, Teuchos::ArrayView<Scalar2> s)
767  {
768  size_t vals_size = vals.size();
769 #ifdef HAVE_AMESOS2_DEBUG
770  size_t s_size = s.size();
771  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
772  std::invalid_argument,
773  "Scale vector must have length at least that of the vector" );
774 #endif
775  size_t i, s_i;
776  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
777  if( s_i == l ){
778  // bring i to the next multiple of ld
779  i += ld - s_i;
780  s_i = 0;
781  }
782  vals[i] *= s[s_i];
783  }
784  }
785 
786  template <typename Scalar1, typename Scalar2, class BinaryOp>
787  void
788  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
789  size_t ld, Teuchos::ArrayView<Scalar2> s,
790  BinaryOp binary_op)
791  {
792  size_t vals_size = vals.size();
793 #ifdef HAVE_AMESOS2_DEBUG
794  size_t s_size = s.size();
795  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
796  std::invalid_argument,
797  "Scale vector must have length at least that of the vector" );
798 #endif
799  size_t i, s_i;
800  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
801  if( s_i == l ){
802  // bring i to the next multiple of ld
803  i += ld - s_i;
804  s_i = 0;
805  }
806  vals[i] = binary_op(vals[i], s[s_i]);
807  }
808  }
809 
812  } // end namespace Util
813 
814 } // end namespace Amesos2
815 
816 #endif // #ifndef AMESOS2_UTIL_HPP
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
static void do_get(const Teuchos::Ptr< const Matrix > mat, const ArrayView< S > nzvals, const ArrayView< GO > indices, const ArrayView< GS > pointers, GS &nnz, const Teuchos::Ptr< const Tpetra::Map< typename Matrix::local_ordinal_t, typename Matrix::global_ordinal_t, typename Matrix::node_t > > map, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:447
const Teuchos::RCP< const Tpetra::Map< LO, GO, Node > > getDistributionMap(EDistribution distribution, GS num_global_elements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, GO indexBase=0)
Gets a Tpetra::Map described by the EDistribution.
Definition: Amesos2_Util.hpp:603
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:591
Provides some simple meta-programming utilities for Amesos2.
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:404
Definition: Amesos2_TypeDecl.hpp:142
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:140
Definition: Amesos2_TypeDecl.hpp:125
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition: Amesos2_Util.cpp:123
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
Definition: Amesos2_TypeDecl.hpp:126
static void do_get(const Teuchos::Ptr< const Matrix > mat, const ArrayView< S > nzvals, const ArrayView< GO > indices, const ArrayView< GS > pointers, GS &nnz, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:430
Definition: Amesos2_TypeDecl.hpp:124
Enum and other types declarations for Amesos2.
Definition: Amesos2_TypeDecl.hpp:127
EDistribution
Definition: Amesos2_TypeDecl.hpp:123