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 #include <Tpetra_DistObject_decl.hpp>
64 #include <Tpetra_ComputeGatherMap.hpp> // added for gather map... where is the best place??
65 
66 #include "Amesos2_TypeDecl.hpp"
67 #include "Amesos2_Meta.hpp"
68 
69 #ifdef HAVE_TPETRA_INST_INT_INT
70 #ifdef HAVE_AMESOS2_EPETRA
71 #include <Epetra_Map.h>
72 #endif
73 #endif
74 
75 
76 namespace Amesos2 {
77 
78  namespace Util {
79 
86  using Teuchos::RCP;
87  using Teuchos::ArrayView;
88 
89  using Meta::is_same;
90  using Meta::if_then_else;
91 
108  template <typename LO, typename GO, typename GS, typename Node>
109  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
110  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map );
111 
112 
113  template <typename LO, typename GO, typename GS, typename Node>
114  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
115  getDistributionMap(EDistribution distribution,
116  GS num_global_elements,
117  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
118  GO indexBase = 0,
119  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map = Teuchos::null);
120 
121 
122 #ifdef HAVE_TPETRA_INST_INT_INT
123 #ifdef HAVE_AMESOS2_EPETRA
124 
130  template <typename LO, typename GO, typename GS, typename Node>
131  RCP<Tpetra::Map<LO,GO,Node> >
132  epetra_map_to_tpetra_map(const Epetra_BlockMap& map);
133 
139  template <typename LO, typename GO, typename GS, typename Node>
140  RCP<Epetra_Map>
141  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map);
142 
148  const RCP<const Teuchos::Comm<int> > to_teuchos_comm(RCP<const Epetra_Comm> c);
149 
155  const RCP<const Epetra_Comm> to_epetra_comm(RCP<const Teuchos::Comm<int> > c);
156 #endif // HAVE_AMESOS2_EPETRA
157 #endif // HAVE_TPETRA_INST_INT_INT
158 
164  template <typename Scalar,
165  typename GlobalOrdinal,
166  typename GlobalSizeT>
167  void transpose(ArrayView<Scalar> vals,
168  ArrayView<GlobalOrdinal> indices,
169  ArrayView<GlobalSizeT> ptr,
170  ArrayView<Scalar> trans_vals,
171  ArrayView<GlobalOrdinal> trans_indices,
172  ArrayView<GlobalSizeT> trans_ptr);
173 
187  template <typename Scalar1, typename Scalar2>
188  void scale(ArrayView<Scalar1> vals, size_t l,
189  size_t ld, ArrayView<Scalar2> s);
190 
209  template <typename Scalar1, typename Scalar2, class BinaryOp>
210  void scale(ArrayView<Scalar1> vals, size_t l,
211  size_t ld, ArrayView<Scalar2> s, BinaryOp binary_op);
212 
213 
215  void printLine( Teuchos::FancyOStream &out );
216 
217 
219  // Matrix/MultiVector Utilities //
221 
222 #ifndef DOXYGEN_SHOULD_SKIP_THIS
223  /*
224  * The following represents a general way of getting a CRS or CCS
225  * representation of a matrix with implicit type conversions. The
226  * \c get_crs_helper and \c get_ccs_helper classes are templated
227  * on 4 types:
228  *
229  * - A Matrix type (conforming to the Amesos2 MatrixAdapter interface)
230  * - A scalar type
231  * - A global ordinal type
232  * - A global size type
233  *
234  * The last three template types correspond to the input argument
235  * types. For example, if the scalar type is \c double , then we
236  * require that the \c nzvals argument is a \c
237  * Teuchos::ArrayView<double> type.
238  *
239  * These helpers perform any type conversions that must be
240  * performed to go between the Matrix's types and the input types.
241  * If no conversions are necessary the static functions can be
242  * effectively inlined.
243  */
244 
245  template <class M, typename S, typename GO, typename GS, class Op>
246  struct same_gs_helper
247  {
248  static void do_get(const Teuchos::Ptr<const M> mat,
249  const ArrayView<typename M::scalar_t> nzvals,
250  const ArrayView<typename M::global_ordinal_t> indices,
251  const ArrayView<GS> pointers,
252  GS& nnz,
253  const Teuchos::Ptr<
254  const Tpetra::Map<typename M::local_ordinal_t,
255  typename M::global_ordinal_t,
256  typename M::node_t> > map,
257  EDistribution distribution,
258  EStorage_Ordering ordering)
259  {
260  Op::apply(mat, nzvals, indices, pointers, nnz, map, distribution, ordering);
261  }
262  };
263 
264  template <class M, typename S, typename GO, typename GS, class Op>
265  struct diff_gs_helper
266  {
267  static void do_get(const Teuchos::Ptr<const M> mat,
268  const ArrayView<typename M::scalar_t> nzvals,
269  const ArrayView<typename M::global_ordinal_t> indices,
270  const ArrayView<GS> pointers,
271  GS& nnz,
272  const Teuchos::Ptr<
273  const Tpetra::Map<typename M::local_ordinal_t,
274  typename M::global_ordinal_t,
275  typename M::node_t> > map,
276  EDistribution distribution,
277  EStorage_Ordering ordering)
278  {
279  typedef typename M::global_size_t mat_gs_t;
280  typename ArrayView<GS>::size_type i, size = pointers.size();
281  Teuchos::Array<mat_gs_t> pointers_tmp(size);
282  mat_gs_t nnz_tmp = 0;
283 
284  Op::apply(mat, nzvals, indices, pointers_tmp, nnz_tmp, map, distribution, ordering);
285 
286  for (i = 0; i < size; ++i){
287  pointers[i] = Teuchos::as<GS>(pointers_tmp[i]);
288  }
289  nnz = Teuchos::as<GS>(nnz_tmp);
290  }
291  };
292 
293  template <class M, typename S, typename GO, typename GS, class Op>
294  struct same_go_helper
295  {
296  static void do_get(const Teuchos::Ptr<const M> mat,
297  const ArrayView<typename M::scalar_t> nzvals,
298  const ArrayView<GO> indices,
299  const ArrayView<GS> pointers,
300  GS& nnz,
301  const Teuchos::Ptr<
302  const Tpetra::Map<typename M::local_ordinal_t,
303  typename M::global_ordinal_t,
304  typename M::node_t> > map,
305  EDistribution distribution,
306  EStorage_Ordering ordering)
307  {
308  typedef typename M::global_size_t mat_gs_t;
309  if_then_else<is_same<GS,mat_gs_t>::value,
310  same_gs_helper<M,S,GO,GS,Op>,
311  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
312  pointers, nnz, map,
313  distribution, ordering);
314  }
315  };
316 
317  template <class M, typename S, typename GO, typename GS, class Op>
318  struct diff_go_helper
319  {
320  static void do_get(const Teuchos::Ptr<const M> mat,
321  const ArrayView<typename M::scalar_t> nzvals,
322  const ArrayView<GO> indices,
323  const ArrayView<GS> pointers,
324  GS& nnz,
325  const Teuchos::Ptr<
326  const Tpetra::Map<typename M::local_ordinal_t,
327  typename M::global_ordinal_t,
328  typename M::node_t> > map,
329  EDistribution distribution,
330  EStorage_Ordering ordering)
331  {
332  typedef typename M::global_ordinal_t mat_go_t;
333  typedef typename M::global_size_t mat_gs_t;
334  typename ArrayView<GO>::size_type i, size = indices.size();
335  Teuchos::Array<mat_go_t> indices_tmp(size);
336 
337  if_then_else<is_same<GS,mat_gs_t>::value,
338  same_gs_helper<M,S,GO,GS,Op>,
339  diff_gs_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices_tmp,
340  pointers, nnz, map,
341  distribution, ordering);
342 
343  for (i = 0; i < size; ++i){
344  indices[i] = Teuchos::as<GO>(indices_tmp[i]);
345  }
346  }
347  };
348 
349  template <class M, typename S, typename GO, typename GS, class Op>
350  struct same_scalar_helper
351  {
352  static void do_get(const Teuchos::Ptr<const M> mat,
353  const ArrayView<S> nzvals,
354  const ArrayView<GO> indices,
355  const ArrayView<GS> pointers,
356  GS& nnz,
357  const Teuchos::Ptr<
358  const Tpetra::Map<typename M::local_ordinal_t,
359  typename M::global_ordinal_t,
360  typename M::node_t> > map,
361  EDistribution distribution,
362  EStorage_Ordering ordering)
363  {
364  typedef typename M::global_ordinal_t mat_go_t;
365  if_then_else<is_same<GO,mat_go_t>::value,
366  same_go_helper<M,S,GO,GS,Op>,
367  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals, indices,
368  pointers, nnz, map,
369  distribution, ordering);
370  }
371  };
372 
373  template <class M, typename S, typename GO, typename GS, class Op>
374  struct diff_scalar_helper
375  {
376  static void do_get(const Teuchos::Ptr<const M> mat,
377  const ArrayView<S> nzvals,
378  const ArrayView<GO> indices,
379  const ArrayView<GS> pointers,
380  GS& nnz,
381  const Teuchos::Ptr<
382  const Tpetra::Map<typename M::local_ordinal_t,
383  typename M::global_ordinal_t,
384  typename M::node_t> > map,
385  EDistribution distribution,
386  EStorage_Ordering ordering)
387  {
388  typedef typename M::scalar_t mat_scalar_t;
389  typedef typename M::global_ordinal_t mat_go_t;
390  typename ArrayView<S>::size_type i, size = nzvals.size();
391  Teuchos::Array<mat_scalar_t> nzvals_tmp(size);
392 
393  if_then_else<is_same<GO,mat_go_t>::value,
394  same_go_helper<M,S,GO,GS,Op>,
395  diff_go_helper<M,S,GO,GS,Op> >::type::do_get(mat, nzvals_tmp, indices,
396  pointers, nnz, map,
397  distribution, ordering);
398 
399  for (i = 0; i < size; ++i){
400  nzvals[i] = Teuchos::as<S>(nzvals_tmp[i]);
401  }
402  }
403  };
404 #endif // DOXYGEN_SHOULD_SKIP_THIS
405 
418  template<class Matrix, typename S, typename GO, typename GS, class Op>
419  struct get_cxs_helper
420  {
421  static void do_get(const Teuchos::Ptr<const Matrix> mat,
422  const ArrayView<S> nzvals,
423  const ArrayView<GO> indices,
424  const ArrayView<GS> pointers,
425  GS& nnz,
426  EDistribution distribution,
427  EStorage_Ordering ordering=ARBITRARY,
428  GO indexBase = 0)
429  {
430  typedef typename Matrix::local_ordinal_t lo_t;
431  typedef typename Matrix::global_ordinal_t go_t;
432  typedef typename Matrix::global_size_t gs_t;
433  typedef typename Matrix::node_t node_t;
434 
435  const Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
436  = getDistributionMap<lo_t,go_t,gs_t,node_t>(distribution,
437  Op::get_dimension(mat),
438  mat->getComm(),
439  indexBase,
440  Op::getMapFromMatrix(mat) //getMap must be the map returned, NOT rowmap or colmap
441  );
442 
443  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
444  }
445 
450  static void do_get(const Teuchos::Ptr<const Matrix> mat,
451  const ArrayView<S> nzvals,
452  const ArrayView<GO> indices,
453  const ArrayView<GS> pointers,
454  GS& nnz,
455  EDistribution distribution, // Does this one need a distribution argument??
456  EStorage_Ordering ordering=ARBITRARY)
457  {
458  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
459  typename Matrix::global_ordinal_t,
460  typename Matrix::node_t> > map
461  = Op::getMap(mat);
462  do_get(mat, nzvals, indices, pointers, nnz, Teuchos::ptrInArg(*map), distribution, ordering);
463  }
464 
469  static void do_get(const Teuchos::Ptr<const Matrix> mat,
470  const ArrayView<S> nzvals,
471  const ArrayView<GO> indices,
472  const ArrayView<GS> pointers,
473  GS& nnz,
474  const Teuchos::Ptr<
475  const Tpetra::Map<typename Matrix::local_ordinal_t,
476  typename Matrix::global_ordinal_t,
477  typename Matrix::node_t> > map,
478  EDistribution distribution,
479  EStorage_Ordering ordering=ARBITRARY)
480  {
481  typedef typename Matrix::scalar_t mat_scalar;
482  if_then_else<is_same<mat_scalar,S>::value,
483  same_scalar_helper<Matrix,S,GO,GS,Op>,
484  diff_scalar_helper<Matrix,S,GO,GS,Op> >::type::do_get(mat,
485  nzvals, indices,
486  pointers, nnz,
487  map,
488  distribution, ordering);
489  }
490  };
491 
492 #ifndef DOXYGEN_SHOULD_SKIP_THIS
493  /*
494  * These two function-like classes are meant to be used as the \c
495  * Op template parameter for the \c get_cxs_helper template class.
496  */
497  template<class Matrix>
498  struct get_ccs_func
499  {
500  static void apply(const Teuchos::Ptr<const Matrix> mat,
501  const ArrayView<typename Matrix::scalar_t> nzvals,
502  const ArrayView<typename Matrix::global_ordinal_t> rowind,
503  const ArrayView<typename Matrix::global_size_t> colptr,
504  typename Matrix::global_size_t& nnz,
505  const Teuchos::Ptr<
506  const Tpetra::Map<typename Matrix::local_ordinal_t,
507  typename Matrix::global_ordinal_t,
508  typename Matrix::node_t> > map,
509  EDistribution distribution,
510  EStorage_Ordering ordering)
511  {
512  mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering, distribution);
513  //mat->getCcs(nzvals, rowind, colptr, nnz, map, ordering);
514  }
515 
516  static
517  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
518  typename Matrix::global_ordinal_t,
519  typename Matrix::node_t> >
520  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
521  {
522  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
523  }
524 
525  static
526  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
527  typename Matrix::global_ordinal_t,
528  typename Matrix::node_t> >
529  getMap(const Teuchos::Ptr<const Matrix> mat)
530  {
531  return mat->getColMap();
532  }
533 
534  static
535  typename Matrix::global_size_t
536  get_dimension(const Teuchos::Ptr<const Matrix> mat)
537  {
538  return mat->getGlobalNumCols();
539  }
540  };
541 
542  template<class Matrix>
543  struct get_crs_func
544  {
545  static void apply(const Teuchos::Ptr<const Matrix> mat,
546  const ArrayView<typename Matrix::scalar_t> nzvals,
547  const ArrayView<typename Matrix::global_ordinal_t> colind,
548  const ArrayView<typename Matrix::global_size_t> rowptr,
549  typename Matrix::global_size_t& nnz,
550  const Teuchos::Ptr<
551  const Tpetra::Map<typename Matrix::local_ordinal_t,
552  typename Matrix::global_ordinal_t,
553  typename Matrix::node_t> > map,
554  EDistribution distribution,
555  EStorage_Ordering ordering)
556  {
557  mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering, distribution);
558  //mat->getCrs(nzvals, colind, rowptr, nnz, map, ordering);
559  }
560 
561  static
562  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
563  typename Matrix::global_ordinal_t,
564  typename Matrix::node_t> >
565  getMapFromMatrix(const Teuchos::Ptr<const Matrix> mat)
566  {
567  return mat->getMap(); // returns Teuchos::null if mat is Epetra_CrsMatrix
568  }
569 
570  static
571  const Teuchos::RCP<const Tpetra::Map<typename Matrix::local_ordinal_t,
572  typename Matrix::global_ordinal_t,
573  typename Matrix::node_t> >
574  getMap(const Teuchos::Ptr<const Matrix> mat)
575  {
576  return mat->getRowMap();
577  }
578 
579  static
580  typename Matrix::global_size_t
581  get_dimension(const Teuchos::Ptr<const Matrix> mat)
582  {
583  return mat->getGlobalNumRows();
584  }
585  };
586 #endif // DOXYGEN_SHOULD_SKIP_THIS
587 
625  template<class Matrix, typename S, typename GO, typename GS>
626  struct get_ccs_helper : get_cxs_helper<Matrix,S,GO,GS,get_ccs_func<Matrix> >
627  {};
628 
636  template<class Matrix, typename S, typename GO, typename GS>
637  struct get_crs_helper : get_cxs_helper<Matrix,S,GO,GS,get_crs_func<Matrix> >
638  {};
639 
640  /* End Matrix/MultiVector Utilities */
641 
642 
644  // Definitions //
646 
647 
648  template <typename LO, typename GO, typename GS, typename Node>
649  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
650  getGatherMap( const Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > &map )
651  {
652  //RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream( Teuchos::null ); // may need to pass an osstream to computeGatherMap for debugging cases...
653  Teuchos::RCP< const Tpetra::Map<LO,GO,Node> > gather_map = Tpetra::Details::computeGatherMap(map, Teuchos::null);
654  return gather_map;
655  }
656 
657 
658  template <typename LO, typename GO, typename GS, typename Node>
659  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
660  getDistributionMap(EDistribution distribution,
661  GS num_global_elements,
662  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
663  GO indexBase,
664  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >& map)
665  {
666  // TODO: Need to add indexBase to cases other than ROOTED
667  // We do not support these maps in any solver now.
668  switch( distribution ){
669  case DISTRIBUTED:
670  case DISTRIBUTED_NO_OVERLAP:
671  return Tpetra::createUniformContigMapWithNode<LO,GO, Node>(num_global_elements, comm);
672  case GLOBALLY_REPLICATED:
673  return Tpetra::createLocalMapWithNode<LO,GO, Node>(num_global_elements, comm);
674  case ROOTED:
675  {
676  int rank = Teuchos::rank(*comm);
677  size_t my_num_elems = Teuchos::OrdinalTraits<size_t>::zero();
678  if( rank == 0 ) { my_num_elems = num_global_elements; }
679 
680  return rcp(new Tpetra::Map<LO,GO, Node>(num_global_elements,
681  my_num_elems, indexBase, comm));
682  }
683  case CONTIGUOUS_AND_ROOTED:
684  {
685  const Teuchos::RCP<const Tpetra::Map<LO,GO,Node> > gathermap
686  = getGatherMap<LO,GO,GS,Node>( map ); //getMap must be the map returned, NOT rowmap or colmap
687  return gathermap;
688  }
689  default:
690  TEUCHOS_TEST_FOR_EXCEPTION( true,
691  std::logic_error,
692  "Control should never reach this point. "
693  "Please contact the Amesos2 developers." );
694  }
695  }
696 
697 
698 #ifdef HAVE_TPETRA_INST_INT_INT
699 #ifdef HAVE_AMESOS2_EPETRA
700 
701  //#pragma message "include 3"
702  //#include <Epetra_Map.h>
703 
704  template <typename LO, typename GO, typename GS, typename Node>
705  Teuchos::RCP<Tpetra::Map<LO,GO,Node> >
706  epetra_map_to_tpetra_map(const Epetra_BlockMap& map)
707  {
708  using Teuchos::as;
709  using Teuchos::rcp;
710 
711  int num_my_elements = map.NumMyElements();
712  Teuchos::Array<int> my_global_elements(num_my_elements);
713  map.MyGlobalElements(my_global_elements.getRawPtr());
714 
715  typedef Tpetra::Map<LO,GO,Node> map_t;
716  RCP<map_t> tmap = rcp(new map_t(Teuchos::OrdinalTraits<GS>::invalid(),
717  my_global_elements(),
718  as<GO>(map.IndexBase()),
719  to_teuchos_comm(Teuchos::rcpFromRef(map.Comm()))));
720  return tmap;
721  }
722 
723  template <typename LO, typename GO, typename GS, typename Node>
724  Teuchos::RCP<Epetra_Map>
725  tpetra_map_to_epetra_map(const Tpetra::Map<LO,GO,Node>& map)
726  {
727  using Teuchos::as;
728 
729  Teuchos::Array<GO> elements_tmp;
730  elements_tmp = map.getNodeElementList();
731  int num_my_elements = elements_tmp.size();
732  Teuchos::Array<int> my_global_elements(num_my_elements);
733  for (int i = 0; i < num_my_elements; ++i){
734  my_global_elements[i] = as<int>(elements_tmp[i]);
735  }
736 
737  using Teuchos::rcp;
738  RCP<Epetra_Map> emap = rcp(new Epetra_Map(-1,
739  num_my_elements,
740  my_global_elements.getRawPtr(),
741  as<GO>(map.getIndexBase()),
742  *to_epetra_comm(map.getComm())));
743  return emap;
744  }
745 #endif // HAVE_AMESOS2_EPETRA
746 #endif // HAVE_TPETRA_INST_INT_INT
747 
748  template <typename Scalar,
749  typename GlobalOrdinal,
750  typename GlobalSizeT>
751  void transpose(Teuchos::ArrayView<Scalar> vals,
752  Teuchos::ArrayView<GlobalOrdinal> indices,
753  Teuchos::ArrayView<GlobalSizeT> ptr,
754  Teuchos::ArrayView<Scalar> trans_vals,
755  Teuchos::ArrayView<GlobalOrdinal> trans_indices,
756  Teuchos::ArrayView<GlobalSizeT> trans_ptr)
757  {
758  /* We have a compressed-row storage format of this matrix. We
759  * transform this into a compressed-column format using a
760  * distribution-counting sort algorithm, which is described by
761  * D.E. Knuth in TAOCP Vol 3, 2nd ed pg 78.
762  */
763 
764 #ifdef HAVE_AMESOS2_DEBUG
765  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_begin, ind_end;
766  ind_begin = indices.begin();
767  ind_end = indices.end();
768  size_t min_trans_ptr_size = *std::max_element(ind_begin, ind_end) + 1;
769  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::as<size_t>(trans_ptr.size()) < min_trans_ptr_size,
770  std::invalid_argument,
771  "Transpose pointer size not large enough." );
772  TEUCHOS_TEST_FOR_EXCEPTION( trans_vals.size() < vals.size(),
773  std::invalid_argument,
774  "Transpose values array not large enough." );
775  TEUCHOS_TEST_FOR_EXCEPTION( trans_indices.size() < indices.size(),
776  std::invalid_argument,
777  "Transpose indices array not large enough." );
778 #else
779  typename Teuchos::ArrayView<GlobalOrdinal>::iterator ind_it, ind_end;
780 #endif
781  // Count the number of entries in each column
782  Teuchos::Array<GlobalSizeT> count(trans_ptr.size(), 0);
783  ind_end = indices.end();
784  for( ind_it = indices.begin(); ind_it != ind_end; ++ind_it ){
785  ++(count[(*ind_it) + 1]);
786  }
787  // Accumulate
788  typename Teuchos::Array<GlobalSizeT>::iterator cnt_it, cnt_end;
789  cnt_end = count.end();
790  for( cnt_it = count.begin() + 1; cnt_it != cnt_end; ++cnt_it ){
791  *cnt_it = *cnt_it + *(cnt_it - 1);
792  }
793  // This becomes the array of column pointers
794  trans_ptr.assign(count);
795 
796  /* Move the nonzero values into their final place in nzval, based on the
797  * counts found previously.
798  *
799  * This sequence deviates from Knuth's algorithm a bit, following more
800  * closely the description presented in Gustavson, Fred G. "Two Fast
801  * Algorithms for Sparse Matrices: Multiplication and Permuted
802  * Transposition" ACM Trans. Math. Softw. volume 4, number 3, 1978, pages
803  * 250--269, http://doi.acm.org/10.1145/355791.355796.
804  *
805  * The output indices end up in sorted order
806  */
807 
808  GlobalSizeT size = ptr.size();
809  for( GlobalSizeT i = 0; i < size - 1; ++i ){
810  GlobalOrdinal u = ptr[i];
811  GlobalOrdinal v = ptr[i + 1];
812  for( GlobalOrdinal j = u; j < v; ++j ){
813  GlobalOrdinal k = count[indices[j]];
814  trans_vals[k] = vals[j];
815  trans_indices[k] = i;
816  ++(count[indices[j]]);
817  }
818  }
819  }
820 
821 
822  template <typename Scalar1, typename Scalar2>
823  void
824  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
825  size_t ld, Teuchos::ArrayView<Scalar2> s)
826  {
827  size_t vals_size = vals.size();
828 #ifdef HAVE_AMESOS2_DEBUG
829  size_t s_size = s.size();
830  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
831  std::invalid_argument,
832  "Scale vector must have length at least that of the vector" );
833 #endif
834  size_t i, s_i;
835  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
836  if( s_i == l ){
837  // bring i to the next multiple of ld
838  i += ld - s_i;
839  s_i = 0;
840  }
841  vals[i] *= s[s_i];
842  }
843  }
844 
845  template <typename Scalar1, typename Scalar2, class BinaryOp>
846  void
847  scale(Teuchos::ArrayView<Scalar1> vals, size_t l,
848  size_t ld, Teuchos::ArrayView<Scalar2> s,
849  BinaryOp binary_op)
850  {
851  size_t vals_size = vals.size();
852 #ifdef HAVE_AMESOS2_DEBUG
853  size_t s_size = s.size();
854  TEUCHOS_TEST_FOR_EXCEPTION( s_size < l,
855  std::invalid_argument,
856  "Scale vector must have length at least that of the vector" );
857 #endif
858  size_t i, s_i;
859  for( i = 0, s_i = 0; i < vals_size; ++i, ++s_i ){
860  if( s_i == l ){
861  // bring i to the next multiple of ld
862  i += ld - s_i;
863  s_i = 0;
864  }
865  vals[i] = binary_op(vals[i], s[s_i]);
866  }
867  }
868 
871  } // end namespace Util
872 
873 } // end namespace Amesos2
874 
875 #endif // #ifndef AMESOS2_UTIL_HPP
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::get_crs_helper
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:637
Amesos2_Meta.hpp
Provides some simple meta-programming utilities for Amesos2.
Amesos2::Util::get_cxs_helper
A generic base class for the CRS and CCS helpers.
Definition: Amesos2_Util.hpp:419
Amesos2::ARBITRARY
Definition: Amesos2_TypeDecl.hpp:143
Amesos2::EStorage_Ordering
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
Amesos2::DISTRIBUTED_NO_OVERLAP
Definition: Amesos2_TypeDecl.hpp:125
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::get_cxs_helper::do_get
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, EDistribution distribution, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:469
Amesos2
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
Amesos2::Util::printLine
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 "-"s on std::cout.
Definition: Amesos2_Util.cpp:123
Amesos2::Util::get_ccs_helper
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:626
Amesos2::GLOBALLY_REPLICATED
Definition: Amesos2_TypeDecl.hpp:126
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:650
Amesos2::Util::get_cxs_helper::do_get
static void do_get(const Teuchos::Ptr< const Matrix > mat, const ArrayView< S > nzvals, const ArrayView< GO > indices, const ArrayView< GS > pointers, GS &nnz, EDistribution distribution, EStorage_Ordering ordering=ARBITRARY)
Definition: Amesos2_Util.hpp:450
Amesos2::DISTRIBUTED
Definition: Amesos2_TypeDecl.hpp:124
Amesos2_TypeDecl.hpp
Enum and other types declarations for Amesos2.
Amesos2::ROOTED
Definition: Amesos2_TypeDecl.hpp:127
Amesos2::EDistribution
EDistribution
Definition: Amesos2_TypeDecl.hpp:123
Amesos2::CONTIGUOUS_AND_ROOTED
Definition: Amesos2_TypeDecl.hpp:128
Generated by   doxygen 1.8.13