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" 
   56   template < 
class Matrix >
 
   57   MatrixAdapter<Matrix>::MatrixAdapter(
const Teuchos::RCP<Matrix> m)
 
   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();
 
   65   template < 
class Matrix >
 
   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,
 
   75     help_getCrs(nzval, colind, rowptr,
 
   76     nnz, rowmap, distribution, ordering,
 
   77     typename adapter_t::get_crs_spec());
 
   80   template < 
class Matrix >
 
   82   MatrixAdapter<Matrix>::getCrs(
const Teuchos::ArrayView<scalar_t> nzval,
 
   83         const Teuchos::ArrayView<global_ordinal_t> colind,
 
   84         const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> rowptr,
 
   85         typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
   89     const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap
 
   90       = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
 
   91                       this->getGlobalNumRows(),
 
   93     this->getCrs(nzval, colind, rowptr, nnz, Teuchos::ptrInArg(*rowmap), ordering, distribution);
 
   96   template < 
class Matrix >
 
   98   MatrixAdapter<Matrix>::getCcs(
const Teuchos::ArrayView<scalar_t> nzval,
 
   99         const Teuchos::ArrayView<global_ordinal_t> rowind,
 
  100         const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> colptr,
 
  101         typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  102         const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t, global_ordinal_t, node_t> > colmap,
 
  106     help_getCcs(nzval, rowind, colptr,
 
  107     nnz, colmap, distribution, ordering,
 
  108     typename adapter_t::get_ccs_spec());
 
  111   template < 
class Matrix >
 
  113   MatrixAdapter<Matrix>::getCcs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  114         const Teuchos::ArrayView<global_ordinal_t> rowind,
 
  115         const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> colptr,
 
  116         typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  120     const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap
 
  121       = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
 
  122                       this->getGlobalNumCols(),
 
  124     this->getCcs(nzval, rowind, colptr, nnz, Teuchos::ptrInArg(*colmap), ordering, distribution);
 
  127   template < 
class Matrix >
 
  128   typename MatrixAdapter<Matrix>::global_size_t
 
  131     return static_cast<const adapter_t*
>(
this)->getGlobalNumRows_impl();
 
  134   template < 
class Matrix >
 
  135   typename MatrixAdapter<Matrix>::global_size_t
 
  138     return static_cast<const adapter_t*
>(
this)->getGlobalNumCols_impl();
 
  141   template < 
class Matrix >
 
  142   typename MatrixAdapter<Matrix>::global_size_t
 
  145     return row_map_->getIndexBase();
 
  148   template < 
class Matrix >
 
  149   typename MatrixAdapter<Matrix>::global_size_t
 
  152     return col_map_->getIndexBase();
 
  155   template < 
class Matrix >
 
  156   typename MatrixAdapter<Matrix>::global_size_t
 
  159     return static_cast<const adapter_t*
>(
this)->getGlobalNNZ_impl();
 
  162   template < 
class Matrix >
 
  166     return row_map_->getNodeNumElements(); 
 
  169   template < 
class Matrix >
 
  173     return col_map_->getNodeNumElements();
 
  176   template < 
class Matrix >
 
  180     return static_cast<const adapter_t*
>(
this)->getLocalNNZ_impl();
 
  184   template < 
class Matrix >
 
  188     std::ostringstream oss;
 
  189     oss << 
"Amesos2::MatrixAdapter wrapping: ";
 
  194   template < 
class Matrix >
 
  197           const Teuchos::EVerbosityLevel verbLevel)
 const 
  200   template < 
class Matrix >
 
  201   typename MatrixAdapter<Matrix>::spmtx_ptr_t
 
  204     return static_cast<const adapter_t*
>(
this)->getSparseRowPtr();
 
  207   template < 
class Matrix >
 
  208   typename MatrixAdapter<Matrix>::spmtx_idx_t
 
  211     return static_cast<const adapter_t*
>(
this)->getSparseColInd();
 
  214   template < 
class Matrix >
 
  215   typename MatrixAdapter<Matrix>::spmtx_vals_t
 
  218     return static_cast<const adapter_t*
>(
this)->getSparseValues();
 
  226   template < 
class Matrix >
 
  229              const Teuchos::ArrayView<global_ordinal_t> colind,
 
  230              const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> rowptr,
 
  231              typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  232              const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
 
  237     static_cast<const adapter_t*
>(
this)->getCrs_spec(nzval, colind, rowptr,
 
  238                  nnz, rowmap, ordering);
 
  241   template < 
class Matrix >
 
  243   MatrixAdapter<Matrix>::help_getCrs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  244              const Teuchos::ArrayView<global_ordinal_t> colind,
 
  245              const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> rowptr,
 
  246              typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  247              const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
 
  250              no_special_impl nsi)
 const 
  255     do_getCrs(nzval, colind, rowptr,
 
  256         nnz, rowmap, distribution, ordering,
 
  257         typename adapter_t::major_access());
 
  260   template < 
class Matrix >
 
  262   MatrixAdapter<Matrix>::do_getCrs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  263            const Teuchos::ArrayView<global_ordinal_t> colind,
 
  264            const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> rowptr,
 
  265            typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  266            const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
 
  273     using Teuchos::ArrayView;
 
  274     using Teuchos::OrdinalTraits;
 
  279     RCP<const type> get_mat;
 
  280     if( *rowmap == *this->row_map_ && distribution != CONTIGUOUS_AND_ROOTED ){
 
  282       get_mat = rcp(
this,
false); 
 
  284       get_mat = 
get(rowmap, distribution);
 
  295     RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rmap = get_mat->getRowMap();
 
  296     ArrayView<const global_ordinal_t> node_elements = rmap->getNodeElementList();
 
  297     if( node_elements.size() == 0 ) 
return; 
 
  299     typename ArrayView<const global_ordinal_t>::iterator row_it, row_end;
 
  300     row_end = node_elements.end();
 
  302     size_t rowptr_ind = OrdinalTraits<size_t>::zero();
 
  303     global_ordinal_t rowInd = OrdinalTraits<global_ordinal_t>::zero();
 
  304     for( row_it = node_elements.begin(); row_it != row_end; ++row_it ){
 
  305       rowptr[rowptr_ind++] = rowInd;
 
  306       size_t rowNNZ = get_mat->getGlobalRowNNZ(*row_it);
 
  307       size_t nnzRet = OrdinalTraits<size_t>::zero();
 
  308       ArrayView<global_ordinal_t> colind_view = colind.view(rowInd,rowNNZ); 
 
  309       ArrayView<scalar_t> nzval_view = nzval.view(rowInd,rowNNZ);
 
  311       get_mat->getGlobalRowCopy(*row_it, colind_view, nzval_view, nnzRet);
 
  312       for (
size_t rr = 0; rr < nnzRet ; rr++)
 
  314           colind_view[rr] = colind_view[rr] - rmap->getIndexBase();
 
  320       if( ordering == SORTED_INDICES ){
 
  321         Tpetra::sort2(colind_view.begin(), colind_view.end(), nzval_view.begin());
 
  324       TEUCHOS_TEST_FOR_EXCEPTION( rowNNZ != nnzRet,
 
  326         "Number of values returned different from " 
  327                           "number of values reported");
 
  330     rowptr[rowptr_ind] = nnz = rowInd;
 
  334   template < 
class Matrix >
 
  336   MatrixAdapter<Matrix>::do_getCrs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  337            const Teuchos::ArrayView<global_ordinal_t> colind,
 
  338            const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> rowptr,
 
  339            typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  340            const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
 
  345     using Teuchos::Array;
 
  348     Array<scalar_t> nzval_tmp(nzval.size(), 0);
 
  349     Array<global_ordinal_t> rowind(colind.size(), 0);
 
  350     Array<global_size_t> colptr(this->getGlobalNumCols() + 1);
 
  351     this->getCcs(nzval_tmp(), rowind(), colptr(), nnz, rowmap, ordering, distribution);
 
  353     if( !nzval.is_null() && !colind.is_null() && !rowptr.is_null() )
 
  354       Util::transpose(nzval_tmp(), rowind(), colptr(), nzval, colind, rowptr);
 
  357   template < 
class Matrix >
 
  359   MatrixAdapter<Matrix>::help_getCcs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  360              const Teuchos::ArrayView<global_ordinal_t> rowind,
 
  361              const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> colptr,
 
  362              typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  363              const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
 
  366              has_special_impl)
 const 
  368     static_cast<const adapter_t*
>(
this)->getCcs_spec(nzval, rowind, colptr,
 
  369                  nnz, colmap, ordering);
 
  372   template < 
class Matrix >
 
  374   MatrixAdapter<Matrix>::help_getCcs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  375              const Teuchos::ArrayView<global_ordinal_t> rowind,
 
  376              const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> colptr,
 
  377              typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  378              const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
 
  381              no_special_impl nsi)
 const 
  385     do_getCcs(nzval, rowind, colptr,
 
  386         nnz, colmap, distribution, ordering,
 
  387         typename adapter_t::major_access());
 
  390   template < 
class Matrix >
 
  392   MatrixAdapter<Matrix>::do_getCcs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  393            const Teuchos::ArrayView<global_ordinal_t> rowind,
 
  394            const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> colptr,
 
  395            typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  396            const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
 
  401     using Teuchos::Array;
 
  406     Array<scalar_t> nzval_tmp(nzval.size(), 0);
 
  407     Array<global_ordinal_t> colind(rowind.size(), 0);
 
  408     Array<global_size_t> rowptr(this->getGlobalNumRows() + 1);
 
  409     this->getCrs(nzval_tmp(), colind(), rowptr(), nnz, colmap, ordering, distribution);
 
  411     if( !nzval.is_null() && !rowind.is_null() && !colptr.is_null() )
 
  412       Util::transpose(nzval_tmp(), colind(), rowptr(), nzval, rowind, colptr);
 
  415   template < 
class Matrix >
 
  417   MatrixAdapter<Matrix>::do_getCcs(
const Teuchos::ArrayView<scalar_t> nzval,
 
  418            const Teuchos::ArrayView<global_ordinal_t> rowind,
 
  419            const Teuchos::ArrayView<
typename MatrixAdapter<Matrix>::global_size_t> colptr,
 
  420            typename MatrixAdapter<Matrix>::global_size_t& nnz,
 
  421            const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
 
  427     using Teuchos::ArrayView;
 
  428     using Teuchos::OrdinalTraits;
 
  430     RCP<const type> get_mat;
 
  431     if( *colmap == *this->col_map_ && distribution != CONTIGUOUS_AND_ROOTED ){
 
  433       get_mat = rcp(
this,
false); 
 
  435       get_mat = 
get(colmap, distribution);
 
  439     RCP<const Tpetra::Map<scalar_t,local_ordinal_t,global_ordinal_t> > cmap = get_mat->getColMap();
 
  440     TEUCHOS_ASSERT( *colmap == *cmap );
 
  442     ArrayView<global_ordinal_t> node_elements = cmap->getNodeElementList();
 
  443     if( node_elements.size() == 0 ) 
return; 
 
  445     typename ArrayView<global_ordinal_t>::iterator col_it, col_end;
 
  446     col_end = node_elements.end();
 
  448     size_t colptr_ind = OrdinalTraits<size_t>::zero();
 
  449     global_ordinal_t colInd = OrdinalTraits<global_ordinal_t>::zero();
 
  450     for( col_it = node_elements.begin(); col_it != col_end; ++col_it ){
 
  451       colptr[colptr_ind++] = colInd;
 
  452       size_t colNNZ = getGlobalColNNZ(*col_it);
 
  454       ArrayView<global_ordinal_t> rowind_view = rowind.view(colInd,colNNZ);
 
  455       ArrayView<scalar_t> nzval_view = nzval.view(colInd,colNNZ);
 
  456       getGlobalColCopy(*col_it, rowind_view, nzval_view, nnzRet);
 
  461       if( ordering == SORTED_INDICES ){
 
  462         Tpetra::sort2(rowind_view.begin(), rowind_view.end(), nzval_view.begin());
 
  465       TEUCHOS_TEST_FOR_EXCEPTION( colNNZ != nnzRet,
 
  467         "Number of values returned different from " 
  468                           "number of values reported");
 
  471     colptr[colptr_ind] = nnz = colInd;
 
  476   template < 
class Matrix >
 
  479             const Teuchos::ArrayView<global_ordinal_t>& indices,
 
  480             const Teuchos::ArrayView<scalar_t>& vals,
 
  483     static_cast<const adapter_t*
>(
this)->getGlobalRowCopy_impl(row, indices, vals, nnz);
 
  486   template < 
class Matrix >
 
  489             const Teuchos::ArrayView<global_ordinal_t>& indices,
 
  490             const Teuchos::ArrayView<scalar_t>& vals,
 
  493     static_cast<const adapter_t*
>(
this)->getGlobalColCopy_impl(col, indices, vals, nnz);
 
  496   template < 
class Matrix >
 
  500     return static_cast<const adapter_t*
>(
this)->getMaxRowNNZ_impl();
 
  503   template < 
class Matrix >
 
  505   MatrixAdapter<Matrix>::getMaxColNNZ()
 const 
  507     return static_cast<const adapter_t*
>(
this)->getMaxColNNZ_impl();
 
  510   template < 
class Matrix >
 
  512   MatrixAdapter<Matrix>::getGlobalRowNNZ(global_ordinal_t row)
 const 
  514     return static_cast<const adapter_t*
>(
this)->getGlobalRowNNZ_impl(row);
 
  517   template < 
class Matrix >
 
  519   MatrixAdapter<Matrix>::getLocalRowNNZ(local_ordinal_t row)
 const 
  521     return static_cast<const adapter_t*
>(
this)->getLocalRowNNZ_impl(row);
 
  524   template < 
class Matrix >
 
  526   MatrixAdapter<Matrix>::getGlobalColNNZ(global_ordinal_t col)
 const 
  528     return static_cast<const adapter_t*
>(
this)->getGlobalColNNZ_impl(col);
 
  531   template < 
class Matrix >
 
  533   MatrixAdapter<Matrix>::getLocalColNNZ(local_ordinal_t col)
 const 
  535     return static_cast<const adapter_t*
>(
this)->getLocalColNNZ_impl(col);
 
  538   template < 
class Matrix >
 
  540   MatrixAdapter<Matrix>::isLocallyIndexed()
 const 
  542     return static_cast<const adapter_t*
>(
this)->isLocallyIndexed_impl();
 
  545   template < 
class Matrix >
 
  547   MatrixAdapter<Matrix>::isGloballyIndexed()
 const 
  549     return static_cast<const adapter_t*
>(
this)->isGloballyIndexed_impl();
 
  553   template < 
class Matrix >
 
  554   Teuchos::RCP<const MatrixAdapter<Matrix> >
 
  555   MatrixAdapter<Matrix>::get(
const Teuchos::Ptr<
const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > map, 
EDistribution distribution)
 const 
  557     return static_cast<const adapter_t*
>(
this)->get_impl(map, distribution);
 
  561   template <
class Matrix>
 
  562   Teuchos::RCP<MatrixAdapter<Matrix> >
 
  563   createMatrixAdapter(Teuchos::RCP<Matrix> m){
 
  565     using Teuchos::rcp_const_cast;
 
  567     if(m.is_null()) 
return Teuchos::null;
 
  568     return( rcp(
new ConcreteMatrixAdapter<Matrix>(m)) );
 
  571   template <
class Matrix>
 
  572   Teuchos::RCP<const MatrixAdapter<Matrix> >
 
  573   createConstMatrixAdapter(Teuchos::RCP<const Matrix> m){
 
  575     using Teuchos::rcp_const_cast;
 
  577     if(m.is_null()) 
return Teuchos::null;
 
  578     return( rcp(
new ConcreteMatrixAdapter<Matrix>(rcp_const_cast<Matrix,const Matrix>(m))).getConst() );
 
  583 #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)
std::string description() const 
Returns a short description of this Solver. 
Definition: Amesos2_MatrixAdapter_def.hpp:186
A Matrix adapter interface for Amesos2. 
Definition: Amesos2_MatrixAdapter_decl.hpp:76
EDistribution
Definition: Amesos2_TypeDecl.hpp:123