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
57namespace 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
208 MatrixAdapter<Matrix>::getGlobalNumRows() const
209 {
210 return static_cast<const adapter_t*>(this)->getGlobalNumRows_impl();
211 }
212
213 template < class Matrix >
214 typename MatrixAdapter<Matrix>::global_size_t
215 MatrixAdapter<Matrix>::getGlobalNumCols() const
216 {
217 return static_cast<const adapter_t*>(this)->getGlobalNumCols_impl();
218 }
219
220 template < class Matrix >
221 typename MatrixAdapter<Matrix>::global_size_t
222 MatrixAdapter<Matrix>::getRowIndexBase() const
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
231 MatrixAdapter<Matrix>::getColumnIndexBase() const
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
240 MatrixAdapter<Matrix>::getGlobalNNZ() const
241 {
242 return static_cast<const adapter_t*>(this)->getGlobalNNZ_impl();
243 }
244
245 template < class Matrix >
246 size_t
247 MatrixAdapter<Matrix>::getLocalNumRows() const
248 {
249 return row_map_->getNodeNumElements(); // TODO: verify
250 }
251
252 template < class Matrix >
253 size_t
254 MatrixAdapter<Matrix>::getLocalNumCols() const
255 {
256 return col_map_->getNodeNumElements();
257 }
258
259 template < class Matrix >
260 size_t
261 MatrixAdapter<Matrix>::getLocalNNZ() const
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
269 MatrixAdapter<Matrix>::description() const
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 >
285 typename MatrixAdapter<Matrix>::spmtx_ptr_t
286 MatrixAdapter<Matrix>::returnRowPtr() const
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 >
308 void MatrixAdapter<Matrix>::returnRowPtr_kokkos_view(KV & view) const
309 {
310 return static_cast<const adapter_t*>(this)->getSparseRowPtr_kokkos_view(view);
311 }
312
313 template < class Matrix >
314 template < class KV >
315 void MatrixAdapter<Matrix>::returnColInd_kokkos_view(KV & view) const
316 {
317 return static_cast<const adapter_t*>(this)->getSparseColInd_kokkos_view(view);
318 }
319
320 template < class Matrix >
321 template < class KV >
322 void MatrixAdapter<Matrix>::returnValues_kokkos_view(KV & view) const
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
828 MatrixAdapter<Matrix>::getGlobalRowCopy_kokkos_view(global_ordinal_t row,
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
860 MatrixAdapter<Matrix>::getMaxRowNNZ() const
861 {
862 return static_cast<const adapter_t*>(this)->getMaxRowNNZ_impl();
863 }
864
865 template < class Matrix >
866 size_t
867 MatrixAdapter<Matrix>::getMaxColNNZ() const
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
Amesos2::EDistribution
EDistribution
Definition Amesos2_TypeDecl.hpp:123
Amesos2::EStorage_Ordering
EStorage_Ordering
Definition Amesos2_TypeDecl.hpp:141
Amesos2::MatrixAdapter
A Matrix adapter interface for Amesos2.
Definition Amesos2_MatrixAdapter_decl.hpp:76
Amesos2::MatrixAdapter::description
std::string description() const
Returns a short description of this Solver.
Definition Amesos2_MatrixAdapter_def.hpp:269
Amesos2::has_special_impl
Indicates that the concrete class has a special implementation that should be called.
Definition Amesos2_TypeDecl.hpp:82
Amesos2::no_special_impl
Indicates that the concrete class can use the generic getC{c|r}s methods implemented in MatrixAdapter...
Definition Amesos2_TypeDecl.hpp:91
Amesos2::row_access
Indicates that the object of an adapter provides row access to its data.
Definition Amesos2_TypeDecl.hpp:100
Generated by doxygen 1.9.8