Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_MatrixAdapter_def.hpp
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 
44 
45 #ifndef AMESOS2_MATRIXADAPTER_DEF_HPP
46 #define AMESOS2_MATRIXADAPTER_DEF_HPP
47 #include <Tpetra_Util.hpp>
48 #include "Amesos2_MatrixAdapter_decl.hpp"
49 #include "Amesos2_ConcreteMatrixAdapter_def.hpp"
50 //#include "Amesos2_ConcreteMatrixAdapter.hpp"
51 
52 #define TESTING_AMESOS2_WITH_TPETRA_REMOVE_UVM
53 #if defined(TESTING_AMESOS2_WITH_TPETRA_REMOVE_UVM)
54 #include "KokkosKernels_SparseUtils.hpp"
55 #endif
56 
57 namespace Amesos2 {
58 
59 
60  template < class Matrix >
61  MatrixAdapter<Matrix>::MatrixAdapter(const Teuchos::RCP<Matrix> m)
62  : mat_(m)
63  {
64  comm_ = static_cast<const adapter_t*>(this)->getComm_impl();
65  col_map_ = static_cast<const adapter_t*>(this)->getColMap_impl();
66  row_map_ = static_cast<const adapter_t*>(this)->getRowMap_impl();
67  }
68 
69  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
70  template < class Matrix >
71  void
72  MatrixAdapter<Matrix>::getCrs(const Teuchos::ArrayView<scalar_t> nzval,
73  const Teuchos::ArrayView<global_ordinal_t> colind,
74  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
75  typename MatrixAdapter<Matrix>::global_size_t& nnz,
76  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t, global_ordinal_t, node_t> > rowmap,
77  EStorage_Ordering ordering,
78  EDistribution distribution) const
79  {
80  help_getCrs(nzval, colind, rowptr,
81  nnz, rowmap, distribution, ordering,
82  typename adapter_t::get_crs_spec());
83  }
84  #endif
85 
86  template < class Matrix >
87  template<typename KV_S, typename KV_GO, typename KV_GS>
88  void
89  MatrixAdapter<Matrix>::getCrs_kokkos_view(KV_S & nzval,
90  KV_GO & colind,
91  KV_GS & rowptr,
92  typename MatrixAdapter<Matrix>::global_size_t& nnz,
93  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t, global_ordinal_t, node_t> > rowmap,
94  EStorage_Ordering ordering,
95  EDistribution distribution) const
96  {
97  help_getCrs_kokkos_view(nzval, colind, rowptr,
98  nnz, rowmap, distribution, ordering,
99  typename adapter_t::get_crs_spec());
100  }
101 
102  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
103  template < class Matrix >
104  void
105  MatrixAdapter<Matrix>::getCrs(const Teuchos::ArrayView<scalar_t> nzval,
106  const Teuchos::ArrayView<global_ordinal_t> colind,
107  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
108  typename MatrixAdapter<Matrix>::global_size_t& nnz,
109  EDistribution distribution,
110  EStorage_Ordering ordering) const
111  {
112  const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap
113  = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
114  this->getGlobalNumRows(),
115  this->getComm());
116  this->getCrs(nzval, colind, rowptr, nnz, Teuchos::ptrInArg(*rowmap), ordering, distribution);
117  }
118  #endif
119 
120  template < class Matrix >
121  template<typename KV_S, typename KV_GO, typename KV_GS>
122  void
123  MatrixAdapter<Matrix>::getCrs_kokkos_view(KV_S & nzval,
124  KV_GO & colind,
125  KV_GS & rowptr,
126  typename MatrixAdapter<Matrix>::global_size_t& nnz,
127  EDistribution distribution,
128  EStorage_Ordering ordering) const
129  {
130  const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap
131  = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
132  this->getGlobalNumRows(),
133  this->getComm());
134  this->getCrs_kokkos_view(nzval, colind, rowptr, nnz, Teuchos::ptrInArg(*rowmap), ordering, distribution);
135  }
136 
137  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
138  template < class Matrix >
139  void
140  MatrixAdapter<Matrix>::getCcs(const Teuchos::ArrayView<scalar_t> nzval,
141  const Teuchos::ArrayView<global_ordinal_t> rowind,
142  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
143  typename MatrixAdapter<Matrix>::global_size_t& nnz,
144  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t, global_ordinal_t, node_t> > colmap,
145  EStorage_Ordering ordering,
146  EDistribution distribution) const
147  {
148  help_getCcs(nzval, rowind, colptr,
149  nnz, colmap, distribution, ordering,
150  typename adapter_t::get_ccs_spec());
151  }
152  #endif
153 
154  template < class Matrix >
155  template<typename KV_S, typename KV_GO, typename KV_GS>
156  void
157  MatrixAdapter<Matrix>::getCcs_kokkos_view(KV_S & nzval,
158  KV_GO & rowind,
159  KV_GS & colptr,
160  typename MatrixAdapter<Matrix>::global_size_t& nnz,
161  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t, global_ordinal_t, node_t> > colmap,
162  EStorage_Ordering ordering,
163  EDistribution distribution) const
164  {
165  help_getCcs_kokkos_view(nzval, rowind, colptr,
166  nnz, colmap, distribution, ordering,
167  typename adapter_t::get_ccs_spec());
168  }
169 
170  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
171  template < class Matrix >
172  void
173  MatrixAdapter<Matrix>::getCcs(const Teuchos::ArrayView<scalar_t> nzval,
174  const Teuchos::ArrayView<global_ordinal_t> rowind,
175  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
176  typename MatrixAdapter<Matrix>::global_size_t& nnz,
177  EDistribution distribution,
178  EStorage_Ordering ordering) const
179  {
180  const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap
181  = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
182  this->getGlobalNumCols(),
183  this->getComm());
184  this->getCcs(nzval, rowind, colptr, nnz, Teuchos::ptrInArg(*colmap), ordering, distribution);
185  }
186  #endif
187 
188  template < class Matrix >
189  template<typename KV_S, typename KV_GO, typename KV_GS>
190  void
191  MatrixAdapter<Matrix>::getCcs_kokkos_view(KV_S & nzval,
192  KV_GO & rowind,
193  KV_GS & colptr,
194  typename MatrixAdapter<Matrix>::global_size_t& nnz,
195  EDistribution distribution,
196  EStorage_Ordering ordering) const
197  {
198  const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap
199  = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
200  this->getGlobalNumCols(),
201  this->getComm());
202  this->getCcs_kokkos_view(nzval, rowind, colptr, nnz, Teuchos::ptrInArg(*colmap), ordering, distribution);
203  }
204 
205 
206  template < class Matrix >
207  typename MatrixAdapter<Matrix>::global_size_t
209  {
210  return static_cast<const adapter_t*>(this)->getGlobalNumRows_impl();
211  }
212 
213  template < class Matrix >
214  typename MatrixAdapter<Matrix>::global_size_t
216  {
217  return static_cast<const adapter_t*>(this)->getGlobalNumCols_impl();
218  }
219 
220  template < class Matrix >
221  typename MatrixAdapter<Matrix>::global_size_t
223  {
224  // Kokkos adapter is for serial only testing right now and will not
225  // create row_map_
226  return (row_map_ != Teuchos::null) ? row_map_->getIndexBase() : 0;
227  }
228 
229  template < class Matrix >
230  typename MatrixAdapter<Matrix>::global_size_t
232  {
233  // Kokkos adapter is for serial only testing right now and will not
234  // create col_map_
235  return (col_map_ != Teuchos::null) ? col_map_->getIndexBase() : 0;
236  }
237 
238  template < class Matrix >
239  typename MatrixAdapter<Matrix>::global_size_t
241  {
242  return static_cast<const adapter_t*>(this)->getGlobalNNZ_impl();
243  }
244 
245  template < class Matrix >
246  size_t
248  {
249  return row_map_->getNodeNumElements(); // TODO: verify
250  }
251 
252  template < class Matrix >
253  size_t
255  {
256  return col_map_->getNodeNumElements();
257  }
258 
259  template < class Matrix >
260  size_t
262  {
263  return static_cast<const adapter_t*>(this)->getLocalNNZ_impl();
264  }
265 
266  // NDE: This is broken for Epetra_CrsMatrix
267  template < class Matrix >
268  std::string
270  {
271  std::ostringstream oss;
272  oss << "Amesos2::MatrixAdapter wrapping: ";
273  oss << mat_->description(); // NDE: This is not defined in Epetra_CrsMatrix, only in Tpetra::CrsMatrix
274  return oss.str();
275  }
276 
277  template < class Matrix >
278  void
279  MatrixAdapter<Matrix>::describe(Teuchos::FancyOStream &out,
280  const Teuchos::EVerbosityLevel verbLevel) const
281  {}
282 
283  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
284  template < class Matrix >
287  {
288  return static_cast<const adapter_t*>(this)->getSparseRowPtr();
289  }
290 
291  template < class Matrix >
292  typename MatrixAdapter<Matrix>::spmtx_idx_t
293  MatrixAdapter<Matrix>::returnColInd() const
294  {
295  return static_cast<const adapter_t*>(this)->getSparseColInd();
296  }
297 
298  template < class Matrix >
299  typename MatrixAdapter<Matrix>::spmtx_vals_t
300  MatrixAdapter<Matrix>::returnValues() const
301  {
302  return static_cast<const adapter_t*>(this)->getSparseValues();
303  }
304  #endif
305 
306  template < class Matrix >
307  template < class KV >
309  {
310  return static_cast<const adapter_t*>(this)->getSparseRowPtr_kokkos_view(view);
311  }
312 
313  template < class Matrix >
314  template < class KV >
316  {
317  return static_cast<const adapter_t*>(this)->getSparseColInd_kokkos_view(view);
318  }
319 
320  template < class Matrix >
321  template < class KV >
323  {
324  return static_cast<const adapter_t*>(this)->getSparseValues_kokkos_view(view);
325  }
326 
327  /******************************
328  * Private method definitions *
329  ******************************/
330 
331  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
332  template < class Matrix >
333  void
334  MatrixAdapter<Matrix>::help_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
335  const Teuchos::ArrayView<global_ordinal_t> colind,
336  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
337  typename MatrixAdapter<Matrix>::global_size_t& nnz,
338  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
339  EDistribution distribution,
340  EStorage_Ordering ordering,
341  has_special_impl) const
342  {
343  static_cast<const adapter_t*>(this)->getCrs_spec(nzval, colind, rowptr,
344  nnz, rowmap, ordering);
345  }
346 
347  template < class Matrix >
348  void
349  MatrixAdapter<Matrix>::help_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
350  const Teuchos::ArrayView<global_ordinal_t> colind,
351  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
352  typename MatrixAdapter<Matrix>::global_size_t& nnz,
353  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
354  EDistribution distribution,
355  EStorage_Ordering ordering,
356  no_special_impl nsi) const
357  {
358 
359  //Added void to remove parameter not used warning
360  ((void)nsi);
361  do_getCrs(nzval, colind, rowptr,
362  nnz, rowmap, distribution, ordering,
363  typename adapter_t::major_access());
364  }
365  #endif
366 
367  template < class Matrix >
368  template<typename KV_S, typename KV_GO, typename KV_GS>
369  void
370  MatrixAdapter<Matrix>::help_getCrs_kokkos_view(KV_S & nzval,
371  KV_GO & colind,
372  KV_GS & rowptr,
373  typename MatrixAdapter<Matrix>::global_size_t& nnz,
374  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
375  EDistribution distribution,
376  EStorage_Ordering ordering,
377  no_special_impl nsi) const
378  {
379 
380  //Added void to remove parameter not used warning
381  ((void)nsi);
382  do_getCrs_kokkos_view(nzval, colind, rowptr,
383  nnz, rowmap, distribution, ordering,
384  typename adapter_t::major_access());
385  }
386 
387  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
388  template < class Matrix >
389  void
390  MatrixAdapter<Matrix>::do_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
391  const Teuchos::ArrayView<global_ordinal_t> colind,
392  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
393  typename MatrixAdapter<Matrix>::global_size_t& nnz,
394  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
395  EDistribution distribution,
396  EStorage_Ordering ordering,
397  row_access ra) const
398  {
399  using Teuchos::rcp;
400  using Teuchos::RCP;
401  using Teuchos::ArrayView;
402  using Teuchos::OrdinalTraits;
403 
404 
405  ((void) ra);
406 
407  RCP<const type> get_mat;
408  if( *rowmap == *this->row_map_ && distribution != CONTIGUOUS_AND_ROOTED ){
409  // No need to redistribute
410  get_mat = rcp(this,false); // non-owning
411  } else {
412  get_mat = get(rowmap, distribution);
413  }
414  // RCP<const type> get_mat = get(rowmap);
415 
416  // rmap may not necessarily check same as rowmap because rmap may
417  // have been constructued with Tpetra's "expert" constructor,
418  // which assumes that the map points are non-contiguous.
419  //
420  // TODO: There may be some more checking between the row map
421  // compatibility, but things are working fine now.
422 
423  RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rmap = get_mat->getRowMap();
424  ArrayView<const global_ordinal_t> node_elements = rmap->getNodeElementList();
425  if( node_elements.size() == 0 ) return; // no more contribution
426 
427  typename ArrayView<const global_ordinal_t>::iterator row_it, row_end;
428  row_end = node_elements.end();
429 
430  size_t rowptr_ind = OrdinalTraits<size_t>::zero();
431  global_ordinal_t rowInd = OrdinalTraits<global_ordinal_t>::zero();
432  for( row_it = node_elements.begin(); row_it != row_end; ++row_it ){
433  rowptr[rowptr_ind++] = rowInd;
434  size_t rowNNZ = get_mat->getGlobalRowNNZ(*row_it);
435  size_t nnzRet = OrdinalTraits<size_t>::zero();
436  ArrayView<global_ordinal_t> colind_view = colind.view(rowInd,rowNNZ);
437  ArrayView<scalar_t> nzval_view = nzval.view(rowInd,rowNNZ);
438 
439  get_mat->getGlobalRowCopy(*row_it, colind_view, nzval_view, nnzRet);
440  for (size_t rr = 0; rr < nnzRet ; rr++)
441  {
442  colind_view[rr] = colind_view[rr] - rmap->getIndexBase();
443  }
444 
445  // It was suggested that instead of sorting each row's indices
446  // individually, that we instead do a double-transpose at the
447  // end, which would also lead to the indices being sorted.
448  if( ordering == SORTED_INDICES ){
449  Tpetra::sort2(colind_view.begin(), colind_view.end(), nzval_view.begin());
450  }
451 
452  TEUCHOS_TEST_FOR_EXCEPTION( rowNNZ != nnzRet,
453  std::runtime_error,
454  "Number of values returned different from "
455  "number of values reported");
456  rowInd += rowNNZ;
457  }
458  rowptr[rowptr_ind] = nnz = rowInd;
459  }
460  #endif
461 
462  template < class Matrix >
463  template<typename KV_S, typename KV_GO, typename KV_GS>
464  void
465  MatrixAdapter<Matrix>::do_getCrs_kokkos_view(KV_S & nzval,
466  KV_GO & colind,
467  KV_GS & rowptr,
468  typename MatrixAdapter<Matrix>::global_size_t& nnz,
469  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
470  EDistribution distribution,
471  EStorage_Ordering ordering,
472  row_access ra) const
473  {
474  // Kokkos adapter will be serial and won't have the rowmap.
475  // Tacho for example wouldn't ever call this in serial but Cholmod will
476  // call ccs and want to convert using this.
477  // If the kokkos adapter is extended to multiple ranks then this will
478  // need to change.
479  if(this->row_map_ == Teuchos::null) {
480  this->returnValues_kokkos_view(nzval);
481  this->returnRowPtr_kokkos_view(rowptr);
482  this->returnColInd_kokkos_view(colind);
483  nnz = nzval.size();
484  return;
485  }
486 
487  using Teuchos::rcp;
488  using Teuchos::RCP;
489  using Teuchos::ArrayView;
490  using Teuchos::OrdinalTraits;
491 
492  ((void) ra);
493 
494  RCP<const type> get_mat;
495  if( *rowmap == *this->row_map_ && distribution != CONTIGUOUS_AND_ROOTED ){
496  // No need to redistribute
497  get_mat = rcp(this,false); // non-owning
498  } else {
499  get_mat = get(rowmap, distribution);
500  }
501  // RCP<const type> get_mat = get(rowmap);
502 
503  // rmap may not necessarily check same as rowmap because rmap may
504  // have been constructued with Tpetra's "expert" constructor,
505  // which assumes that the map points are non-contiguous.
506  //
507  // TODO: There may be some more checking between the row map
508  // compatibility, but things are working fine now.
509 
510  RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rmap = get_mat->getRowMap();
511  ArrayView<const global_ordinal_t> node_elements = rmap->getNodeElementList();
512  //if( node_elements.size() == 0 ) return; // no more contribution
513  typename ArrayView<const global_ordinal_t>::iterator row_it, row_end;
514  row_end = node_elements.end();
515 
516  size_t rowptr_ind = OrdinalTraits<size_t>::zero();
517  global_ordinal_t rowInd = OrdinalTraits<global_ordinal_t>::zero();
518 
519  // For rowptr we can just make a mirror and deep_copy at the end
520  typename KV_GS::HostMirror host_rowptr = Kokkos::create_mirror_view(rowptr);
521 
522  #if !defined(TESTING_AMESOS2_WITH_TPETRA_REMOVE_UVM)
523  // Note nzval, colind, and rowptr will not all be in the same memory space.
524  // Currently only Cholmod exercises this code which has all the arrays on host,
525  // so this will need extension and testing when we have a solver using device here.
526  Kokkos::View<scalar_t*, Kokkos::HostSpace>
527  mat_nzval(Kokkos::ViewAllocateWithoutInitializing("mat_nzval"), nzval.size());
528 
529  Kokkos::View<global_ordinal_t*, Kokkos::HostSpace>
530  mat_colind(Kokkos::ViewAllocateWithoutInitializing("mat_colind"), colind.size());
531 
532  ArrayView<scalar_t> nzval_arrayview(mat_nzval.data(), nzval.size());
533  ArrayView<global_ordinal_t> colind_arrayview(mat_colind.data(), colind.size());
534 
535  for( row_it = node_elements.begin(); row_it != row_end; ++row_it ){
536  host_rowptr(rowptr_ind++) = rowInd;
537  size_t rowNNZ = get_mat->getGlobalRowNNZ(*row_it);
538  size_t nnzRet = OrdinalTraits<size_t>::zero();
539  ArrayView<global_ordinal_t> colind_view = colind_arrayview.view(rowInd,rowNNZ);
540  ArrayView<scalar_t> nzval_view = nzval_arrayview.view(rowInd,rowNNZ);
541 
542  get_mat->getGlobalRowCopy(*row_it, colind_view, nzval_view, nnzRet);
543 
544  for (size_t rr = 0; rr < nnzRet ; rr++) {
545  colind_view[rr] -= rmap->getIndexBase();
546  }
547 
548  // It was suggested that instead of sorting each row's indices
549  // individually, that we instead do a double-transpose at the
550  // end, which would also lead to the indices being sorted.
551  if( ordering == SORTED_INDICES ) {
552  Tpetra::sort2(colind_view.begin(), colind_view.end(), nzval_view.begin());
553  }
554 
555  TEUCHOS_TEST_FOR_EXCEPTION( rowNNZ != nnzRet,
556  std::runtime_error,
557  "Number of values returned different from "
558  "number of values reported");
559  rowInd += rowNNZ;
560  }
561  host_rowptr(rowptr_ind) = nnz = rowInd;
562 
563  deep_copy_or_assign_view(nzval, mat_nzval);
564  deep_copy_or_assign_view(colind, mat_colind);
565  deep_copy_or_assign_view(rowptr, host_rowptr);
566  #else
567  // create temporary views to hold colind and nvals (TODO: allocate as much as needed, also for rowptr)
568  global_host_idx_t mat_colind(Kokkos::ViewAllocateWithoutInitializing("mat_colind"), nzval.size());
569  global_host_val_t mat_nzvals(Kokkos::ViewAllocateWithoutInitializing("mat_nzvals"), colind.size());
570 
571  auto host_colind = Kokkos::create_mirror_view(colind);
572  auto host_nzval = Kokkos::create_mirror_view(nzval);
573 
574  // load crs (on host)
575  for( row_it = node_elements.begin(); row_it != row_end; ++row_it ){
576  size_t rowNNZ = get_mat->getGlobalRowNNZ(*row_it);
577  size_t nnzRet = OrdinalTraits<size_t>::zero();
578  //using range_type = Kokkos::pair<int, int>;
579  //auto colind_view = Kokkos::subview(mat_colind, range_type(rowInd, rowInd+rowNNZ));
580  //auto nzval_view = Kokkos::subview(mat_nzvals, range_type(rowInd, rowInd+rowNNZ));
581  global_host_idx_t colind_view (&(mat_colind(rowInd)), rowNNZ);
582  global_host_val_t nzvals_view (&(mat_nzvals(rowInd)), rowNNZ);
583 
584  global_ordinal_t row_id = *row_it;
585  get_mat->getGlobalRowCopy_kokkos_view(row_id, colind_view, nzvals_view, nnzRet);
586 
587  TEUCHOS_TEST_FOR_EXCEPTION( rowNNZ != nnzRet,
588  std::runtime_error,
589  "Number of values returned different from "
590  "number of values reported");
591  host_rowptr(rowptr_ind++) = rowInd;
592  rowInd += rowNNZ;
593  }
594  host_rowptr(rowptr_ind) = nnz = rowInd;
595 
596  // fix index-base
597  if (rmap->getIndexBase() != 0) {
598  for (size_t k = 0; k < mat_colind.extent(0); k++) {
599  mat_colind(k) -= rmap->getIndexBase();
600  }
601  }
602 
603  // copy to device (note: everything in the vectors are copied, though they may not be used)
604  deep_copy_or_assign_view(nzval, mat_nzvals);
605  deep_copy_or_assign_view(colind, mat_colind);
606  deep_copy_or_assign_view(rowptr, host_rowptr);
607 
608  // sort
609  if( ordering == SORTED_INDICES ) {
610  using execution_space = typename KV_GS::execution_space;
611  KokkosKernels::Impl::sort_crs_matrix <execution_space, KV_GS, KV_GO, KV_S>
612  (rowptr, colind, nzval);
613  }
614  #endif
615  }
616 
617  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
618  // TODO: This may not work with distributed matrices.
619  template < class Matrix >
620  void
621  MatrixAdapter<Matrix>::do_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
622  const Teuchos::ArrayView<global_ordinal_t> colind,
623  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
624  typename MatrixAdapter<Matrix>::global_size_t& nnz,
625  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
626  EDistribution distribution,
627  EStorage_Ordering ordering,
628  col_access ca) const
629  {
630  using Teuchos::Array;
631  // get the ccs and transpose
632 
633  Array<scalar_t> nzval_tmp(nzval.size(), 0);
634  Array<global_ordinal_t> rowind(colind.size(), 0);
635  Array<global_size_t> colptr(this->getGlobalNumCols() + 1);
636  this->getCcs(nzval_tmp(), rowind(), colptr(), nnz, rowmap, ordering, distribution);
637 
638  if( !nzval.is_null() && !colind.is_null() && !rowptr.is_null() )
639  Util::transpose(nzval_tmp(), rowind(), colptr(), nzval, colind, rowptr);
640  }
641 
642  template < class Matrix >
643  void
644  MatrixAdapter<Matrix>::help_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
645  const Teuchos::ArrayView<global_ordinal_t> rowind,
646  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
647  typename MatrixAdapter<Matrix>::global_size_t& nnz,
648  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
649  EDistribution distribution,
650  EStorage_Ordering ordering,
651  has_special_impl) const
652  {
653  static_cast<const adapter_t*>(this)->getCcs_spec(nzval, rowind, colptr,
654  nnz, colmap, ordering);
655  }
656  #endif
657 
658  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
659  template < class Matrix >
660  void
661  MatrixAdapter<Matrix>::help_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
662  const Teuchos::ArrayView<global_ordinal_t> rowind,
663  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
664  typename MatrixAdapter<Matrix>::global_size_t& nnz,
665  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
666  EDistribution distribution,
667  EStorage_Ordering ordering,
668  no_special_impl nsi) const
669  {
670  ((void)nsi);
671 
672  do_getCcs(nzval, rowind, colptr,
673  nnz, colmap, distribution, ordering,
674  typename adapter_t::major_access());
675  }
676  #endif
677 
678  template < class Matrix >
679  template<typename KV_S, typename KV_GO, typename KV_GS>
680  void
681  MatrixAdapter<Matrix>::help_getCcs_kokkos_view(KV_S & nzval,
682  KV_GO & rowind,
683  KV_GS & colptr,
684  typename MatrixAdapter<Matrix>::global_size_t& nnz,
685  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
686  EDistribution distribution,
687  EStorage_Ordering ordering,
688  no_special_impl nsi) const
689  {
690 
691  //Added void to remove parameter not used warning
692  ((void)nsi);
693  do_getCcs_kokkos_view(nzval, rowind, colptr,
694  nnz, colmap, distribution, ordering,
695  typename adapter_t::major_access());
696  }
697 
698  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
699  template < class Matrix >
700  void
701  MatrixAdapter<Matrix>::do_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
702  const Teuchos::ArrayView<global_ordinal_t> rowind,
703  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
704  typename MatrixAdapter<Matrix>::global_size_t& nnz,
705  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
706  EDistribution distribution,
707  EStorage_Ordering ordering,
708  row_access ra) const
709  {
710  using Teuchos::Array;
711  // get the crs and transpose
712 
713  ((void) ra);
714 
715  Array<scalar_t> nzval_tmp(nzval.size(), 0);
716  Array<global_ordinal_t> colind(rowind.size(), 0);
717  Array<global_size_t> rowptr(this->getGlobalNumRows() + 1);
718  this->getCrs(nzval_tmp(), colind(), rowptr(), nnz, colmap, ordering, distribution);
719 
720  if( !nzval.is_null() && !rowind.is_null() && !colptr.is_null() )
721  Util::transpose(nzval_tmp(), colind(), rowptr(), nzval, rowind, colptr);
722  }
723  #endif
724 
725  template < class Matrix >
726  template<typename KV_S, typename KV_GO, typename KV_GS>
727  void
728  MatrixAdapter<Matrix>::do_getCcs_kokkos_view(KV_S & nzval,
729  KV_GO & rowind,
730  KV_GS & colptr,
731  typename MatrixAdapter<Matrix>::global_size_t& nnz,
732  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
733  EDistribution distribution,
734  EStorage_Ordering ordering,
735  row_access ra) const
736  {
737  using Teuchos::ArrayView;
738  // get the crs and transpose
739 
740  ((void) ra);
741 
742  KV_S nzval_tmp(Kokkos::ViewAllocateWithoutInitializing("nzval_tmp"), nzval.size());
743  KV_GO colind(Kokkos::ViewAllocateWithoutInitializing("colind"), rowind.size());
744  KV_GS rowptr(Kokkos::ViewAllocateWithoutInitializing("rowptr"), this->getGlobalNumRows() + 1);
745 
746  this->getCrs_kokkos_view(nzval_tmp, colind, rowptr, nnz, colmap, ordering, distribution);
747 
748  if(nnz > 0) {
749  // This is currently just used by Cholmod in which case the views will be
750  // host, even if Cholmod is using GPU. Will need to upgrade this section
751  // to properly handle device when we have a solver that needs it.
752  ArrayView<typename KV_S::value_type> av_nzval_tmp(nzval_tmp.data(), nzval_tmp.size());
753  ArrayView<typename KV_GO::value_type> av_colind(colind.data(), colind.size());
754  ArrayView<typename KV_GS::value_type> av_rowptr(rowptr.data(), rowptr.size());
755  ArrayView<typename KV_S::value_type> av_nzval(nzval.data(), nzval.size());
756  ArrayView<typename KV_GO::value_type> av_rowind(rowind.data(), rowind.size());
757  ArrayView<typename KV_GS::value_type> av_colptr(colptr.data(), colptr.size());
758  Util::transpose(av_nzval_tmp, av_colind, av_rowptr, av_nzval, av_rowind, av_colptr);
759  }
760  }
761 
762  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
763  template < class Matrix >
764  void
765  MatrixAdapter<Matrix>::do_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
766  const Teuchos::ArrayView<global_ordinal_t> rowind,
767  const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
768  typename MatrixAdapter<Matrix>::global_size_t& nnz,
769  const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
770  EDistribution distribution,
771  EStorage_Ordering ordering,
772  col_access ca) const
773  {
774  using Teuchos::RCP;
775  using Teuchos::ArrayView;
776  using Teuchos::OrdinalTraits;
777 
778  RCP<const type> get_mat;
779  if( *colmap == *this->col_map_ && distribution != CONTIGUOUS_AND_ROOTED ){
780  // No need to redistribute
781  get_mat = rcp(this,false); // non-owning
782  } else {
783  get_mat = get(colmap, distribution);
784  }
785 
786  // If all is well and good, then colmap == cmap
787  RCP<const Tpetra::Map<scalar_t,local_ordinal_t,global_ordinal_t> > cmap = get_mat->getColMap();
788  TEUCHOS_ASSERT( *colmap == *cmap );
789 
790  ArrayView<global_ordinal_t> node_elements = cmap->getNodeElementList();
791  if( node_elements.size() == 0 ) return; // no more contribution
792 
793  typename ArrayView<global_ordinal_t>::iterator col_it, col_end;
794  col_end = node_elements.end();
795 
796  size_t colptr_ind = OrdinalTraits<size_t>::zero();
797  global_ordinal_t colInd = OrdinalTraits<global_ordinal_t>::zero();
798  for( col_it = node_elements.begin(); col_it != col_end; ++col_it ){
799  colptr[colptr_ind++] = colInd;
800  size_t colNNZ = getGlobalColNNZ(*col_it);
801  size_t nnzRet = 0;
802  ArrayView<global_ordinal_t> rowind_view = rowind.view(colInd,colNNZ);
803  ArrayView<scalar_t> nzval_view = nzval.view(colInd,colNNZ);
804  getGlobalColCopy(*col_it, rowind_view, nzval_view, nnzRet);
805 
806  // It was suggested that instead of sorting each row's indices
807  // individually, that we instead do a double-transpose at the
808  // end, which would also lead to the indices being sorted.
809  if( ordering == SORTED_INDICES ){
810  Tpetra::sort2(rowind_view.begin(), rowind_view.end(), nzval_view.begin());
811  }
812 
813  TEUCHOS_TEST_FOR_EXCEPTION( colNNZ != nnzRet,
814  std::runtime_error,
815  "Number of values returned different from "
816  "number of values reported");
817  colInd += colNNZ;
818  }
819  colptr[colptr_ind] = nnz = colInd;
820  }
821  #endif
822 
823 
824  // These will link to concrete implementations
825  template < class Matrix >
826  template<typename KV_GO, typename KV_S>
827  void
829  KV_GO & indices,
830  KV_S & vals,
831  size_t& nnz) const
832  {
833  static_cast<const adapter_t*>(this)->getGlobalRowCopy_kokkos_view_impl(row, indices, vals, nnz);
834  }
835 
836  #ifdef TPETRA_ENABLE_DEPRECATED_CODE
837  template < class Matrix >
838  void
839  MatrixAdapter<Matrix>::getGlobalRowCopy(global_ordinal_t row,
840  const Teuchos::ArrayView<global_ordinal_t>& indices,
841  const Teuchos::ArrayView<scalar_t>& vals,
842  size_t& nnz) const
843  {
844  static_cast<const adapter_t*>(this)->getGlobalRowCopy_impl(row, indices, vals, nnz);
845  }
846 
847  template < class Matrix >
848  void
849  MatrixAdapter<Matrix>::getGlobalColCopy(global_ordinal_t col,
850  const Teuchos::ArrayView<global_ordinal_t>& indices,
851  const Teuchos::ArrayView<scalar_t>& vals,
852  size_t& nnz) const
853  {
854  static_cast<const adapter_t*>(this)->getGlobalColCopy_impl(col, indices, vals, nnz);
855  }
856  #endif
857 
858  template < class Matrix >
859  size_t
861  {
862  return static_cast<const adapter_t*>(this)->getMaxRowNNZ_impl();
863  }
864 
865  template < class Matrix >
866  size_t
868  {
869  return static_cast<const adapter_t*>(this)->getMaxColNNZ_impl();
870  }
871 
872  template < class Matrix >
873  size_t
874  MatrixAdapter<Matrix>::getGlobalRowNNZ(global_ordinal_t row) const
875  {
876  return static_cast<const adapter_t*>(this)->getGlobalRowNNZ_impl(row);
877  }
878 
879  template < class Matrix >
880  size_t
881  MatrixAdapter<Matrix>::getLocalRowNNZ(local_ordinal_t row) const
882  {
883  return static_cast<const adapter_t*>(this)->getLocalRowNNZ_impl(row);
884  }
885 
886  template < class Matrix >
887  size_t
888  MatrixAdapter<Matrix>::getGlobalColNNZ(global_ordinal_t col) const
889  {
890  return static_cast<const adapter_t*>(this)->getGlobalColNNZ_impl(col);
891  }
892 
893  template < class Matrix >
894  size_t
895  MatrixAdapter<Matrix>::getLocalColNNZ(local_ordinal_t col) const
896  {
897  return static_cast<const adapter_t*>(this)->getLocalColNNZ_impl(col);
898  }
899 
900  template < class Matrix >
901  bool
902  MatrixAdapter<Matrix>::isLocallyIndexed() const
903  {
904  return static_cast<const adapter_t*>(this)->isLocallyIndexed_impl();
905  }
906 
907  template < class Matrix >
908  bool
909  MatrixAdapter<Matrix>::isGloballyIndexed() const
910  {
911  return static_cast<const adapter_t*>(this)->isGloballyIndexed_impl();
912  }
913 
914 
915  template < class Matrix >
916  Teuchos::RCP<const MatrixAdapter<Matrix> >
917  MatrixAdapter<Matrix>::get(const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > map, EDistribution distribution) const
918  {
919  return static_cast<const adapter_t*>(this)->get_impl(map, distribution);
920  }
921 
922 
923  template <class Matrix>
924  Teuchos::RCP<MatrixAdapter<Matrix> >
925  createMatrixAdapter(Teuchos::RCP<Matrix> m){
926  using Teuchos::rcp;
927  using Teuchos::rcp_const_cast;
928 
929  if(m.is_null()) return Teuchos::null;
930  return( rcp(new ConcreteMatrixAdapter<Matrix>(m)) );
931  }
932 
933  template <class Matrix>
934  Teuchos::RCP<const MatrixAdapter<Matrix> >
935  createConstMatrixAdapter(Teuchos::RCP<const Matrix> m){
936  using Teuchos::rcp;
937  using Teuchos::rcp_const_cast;
938 
939  if(m.is_null()) return Teuchos::null;
940  return( rcp(new ConcreteMatrixAdapter<Matrix>(rcp_const_cast<Matrix,const Matrix>(m))).getConst() );
941  }
942 
943 } // end namespace Amesos2
944 
945 #endif // AMESOS2_MATRIXADAPTER_DEF_HPP
EStorage_Ordering
Definition: Amesos2_TypeDecl.hpp:141
Indicates that the concrete class has a special implementation that should be called.
Definition: Amesos2_TypeDecl.hpp:82
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
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
std::string description() const
Returns a short description of this Solver.
Definition: Amesos2_MatrixAdapter_def.hpp:269
EDistribution
Definition: Amesos2_TypeDecl.hpp:123