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