Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_MatrixAdapter_def.hpp
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //           Amesos2: Templated Direct Sparse Solver Package 
00006 //                  Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //
00042 // @HEADER
00043 
00044 
00045 #ifndef AMESOS2_MATRIXADAPTER_DEF_HPP
00046 #define AMESOS2_MATRIXADAPTER_DEF_HPP
00047 
00048 #include "Amesos2_MatrixAdapter_decl.hpp"
00049 #include "Amesos2_ConcreteMatrixAdapter_def.hpp"
00050 //#include "Amesos2_ConcreteMatrixAdapter.hpp"
00051 
00052 
00053 namespace Amesos2 {
00054 
00055   
00056   template < class Matrix >
00057   MatrixAdapter<Matrix>::MatrixAdapter(const Teuchos::RCP<Matrix> m)
00058     : mat_(m)
00059   {
00060     comm_ = static_cast<const adapter_t*>(this)->getComm_impl();
00061     col_map_ = static_cast<const adapter_t*>(this)->getColMap_impl();
00062     row_map_ = static_cast<const adapter_t*>(this)->getRowMap_impl();
00063   }
00064 
00065   template < class Matrix >
00066   void
00067   MatrixAdapter<Matrix>::getCrs(const Teuchos::ArrayView<scalar_t> nzval,
00068         const Teuchos::ArrayView<global_ordinal_t> colind,
00069         const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
00070         typename MatrixAdapter<Matrix>::global_size_t& nnz,
00071         const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t, global_ordinal_t, node_t> > rowmap,
00072         EStorage_Ordering ordering) const
00073   {
00074     help_getCrs(nzval, colind, rowptr,
00075     nnz, rowmap, ordering,
00076     typename adapter_t::get_crs_spec());
00077   }
00078 
00079   template < class Matrix >
00080   void
00081   MatrixAdapter<Matrix>::getCrs(const Teuchos::ArrayView<scalar_t> nzval,
00082         const Teuchos::ArrayView<global_ordinal_t> colind,
00083         const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
00084         typename MatrixAdapter<Matrix>::global_size_t& nnz,
00085         EDistribution distribution,
00086         EStorage_Ordering ordering) const
00087   {
00088     const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap
00089       = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
00090                       this->getGlobalNumRows(),
00091                       this->getComm());
00092     this->getCrs(nzval, colind, rowptr, nnz, Teuchos::ptrInArg(*rowmap), ordering);
00093   }
00094 
00095   template < class Matrix >
00096   void
00097   MatrixAdapter<Matrix>::getCcs(const Teuchos::ArrayView<scalar_t> nzval,
00098         const Teuchos::ArrayView<global_ordinal_t> rowind,
00099         const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
00100         typename MatrixAdapter<Matrix>::global_size_t& nnz,
00101         const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t, global_ordinal_t, node_t> > colmap,
00102         EStorage_Ordering ordering) const
00103   {
00104     help_getCcs(nzval, rowind, colptr,
00105     nnz, colmap, ordering,
00106     typename adapter_t::get_ccs_spec());
00107   }
00108 
00109   template < class Matrix >
00110   void
00111   MatrixAdapter<Matrix>::getCcs(const Teuchos::ArrayView<scalar_t> nzval,
00112         const Teuchos::ArrayView<global_ordinal_t> rowind,
00113         const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
00114         typename MatrixAdapter<Matrix>::global_size_t& nnz,
00115         EDistribution distribution,
00116         EStorage_Ordering ordering) const
00117   {
00118     const Teuchos::RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap
00119       = Util::getDistributionMap<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(distribution,
00120                       this->getGlobalNumCols(),
00121                       this->getComm());
00122     this->getCcs(nzval, rowind, colptr, nnz, Teuchos::ptrInArg(*colmap), ordering);
00123   }
00124 
00125   template < class Matrix >
00126   typename MatrixAdapter<Matrix>::global_size_t
00127   MatrixAdapter<Matrix>::getGlobalNumRows() const
00128   {
00129     return row_map_->getMaxAllGlobalIndex() + 1;
00130   }
00131 
00132   template < class Matrix >
00133   typename MatrixAdapter<Matrix>::global_size_t
00134   MatrixAdapter<Matrix>::getGlobalNumCols() const
00135   {
00136     return col_map_->getMaxAllGlobalIndex() + 1;
00137   }
00138   
00139   template < class Matrix >
00140   typename MatrixAdapter<Matrix>::global_size_t
00141   MatrixAdapter<Matrix>::getGlobalNNZ() const
00142   {
00143     return static_cast<const adapter_t*>(this)->getGlobalNNZ_impl();
00144   }
00145   
00146   template < class Matrix >
00147   size_t
00148   MatrixAdapter<Matrix>::getLocalNumRows() const
00149   {
00150     return row_map_->getNodeNumElements(); // TODO: verify
00151   }
00152 
00153   template < class Matrix >
00154   size_t
00155   MatrixAdapter<Matrix>::getLocalNumCols() const
00156   {
00157     return col_map_->getNodeNumElements();
00158   }
00159   
00160   template < class Matrix >
00161   size_t
00162   MatrixAdapter<Matrix>::getLocalNNZ() const
00163   {
00164     return static_cast<const adapter_t*>(this)->getLocalNNZ_impl();
00165   }
00166 
00167   template < class Matrix >
00168   std::string
00169   MatrixAdapter<Matrix>::description() const
00170   {
00171     std::ostringstream oss;
00172     oss << "Amesos2::MatrixAdapter wrapping: ";
00173     oss << mat_->description();
00174     return oss.str();
00175   }
00176   
00177   template < class Matrix >
00178   void
00179   MatrixAdapter<Matrix>::describe(Teuchos::FancyOStream &out,
00180           const Teuchos::EVerbosityLevel verbLevel) const
00181   {}
00182   
00183 
00184 
00185   /******************************
00186    * Private method definitions *
00187    ******************************/
00188 
00189   template < class Matrix >
00190   void
00191   MatrixAdapter<Matrix>::help_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
00192              const Teuchos::ArrayView<global_ordinal_t> colind,
00193              const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
00194              typename MatrixAdapter<Matrix>::global_size_t& nnz,
00195              const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
00196              EStorage_Ordering ordering,
00197              has_special_impl hsi) const
00198   {
00199     static_cast<const adapter_t*>(this)->getCrs_spec(nzval, colind, rowptr,
00200                  nnz, rowmap, ordering);
00201   }
00202 
00203   template < class Matrix >
00204   void
00205   MatrixAdapter<Matrix>::help_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
00206              const Teuchos::ArrayView<global_ordinal_t> colind,
00207              const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
00208              typename MatrixAdapter<Matrix>::global_size_t& nnz,
00209              const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
00210              EStorage_Ordering ordering,
00211              no_special_impl nsi) const
00212   {
00213     do_getCrs(nzval, colind, rowptr,
00214         nnz, rowmap, ordering,
00215         typename adapter_t::major_access());
00216   }
00217 
00218   template < class Matrix >
00219   void
00220   MatrixAdapter<Matrix>::do_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
00221            const Teuchos::ArrayView<global_ordinal_t> colind,
00222            const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
00223            typename MatrixAdapter<Matrix>::global_size_t& nnz,
00224            const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
00225            EStorage_Ordering ordering,
00226            row_access ra) const
00227   {
00228     using Teuchos::rcp;
00229     using Teuchos::RCP;
00230     using Teuchos::ArrayView;
00231     using Teuchos::OrdinalTraits;
00232     
00233     RCP<const type> get_mat;
00234     if( *rowmap == *this->row_map_ ){
00235       // No need to redistribute
00236       get_mat = rcp(this,false); // non-owning
00237     } else {
00238       get_mat = get(rowmap);
00239     }
00240     // RCP<const type> get_mat = get(rowmap);
00241 
00242     // rmap may not necessarily check same as rowmap because rmap may
00243     // have been constructued with Tpetra's "expert" constructor,
00244     // which assumes that the map points are non-contiguous.
00245     //
00246     // TODO: There may be some more checking between the row map
00247     // compatibility, but things are working fine now.
00248     RCP<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rmap = get_mat->getRowMap();
00249 
00250     ArrayView<const global_ordinal_t> node_elements = rmap->getNodeElementList();
00251     if( node_elements.size() == 0 ) return; // no more contribution
00252     
00253     typename ArrayView<const global_ordinal_t>::iterator row_it, row_end;
00254     row_end = node_elements.end();
00255 
00256     size_t rowptr_ind = OrdinalTraits<size_t>::zero();
00257     global_ordinal_t rowInd = OrdinalTraits<global_ordinal_t>::zero();
00258     for( row_it = node_elements.begin(); row_it != row_end; ++row_it ){
00259       rowptr[rowptr_ind++] = rowInd;
00260       size_t rowNNZ = get_mat->getGlobalRowNNZ(*row_it);
00261       size_t nnzRet = OrdinalTraits<size_t>::zero();
00262       ArrayView<global_ordinal_t> colind_view = colind.view(rowInd,rowNNZ);
00263       ArrayView<scalar_t> nzval_view = nzval.view(rowInd,rowNNZ);
00264       
00265       get_mat->getGlobalRowCopy(*row_it, colind_view, nzval_view, nnzRet);
00266 
00267       // It was suggested that instead of sorting each row's indices
00268       // individually, that we instead do a double-transpose at the
00269       // end, which would also lead to the indices being sorted.
00270       if( ordering == SORTED_INDICES ){
00271   Tpetra::sort2(colind_view.begin(), colind_view.end(), nzval_view.begin());
00272       }
00273       
00274       TEUCHOS_TEST_FOR_EXCEPTION( rowNNZ != nnzRet,
00275         std::runtime_error,
00276         "Number of values returned different from "
00277                           "number of values reported");
00278       rowInd += rowNNZ;
00279     }
00280     rowptr[rowptr_ind] = nnz = rowInd;
00281   }
00282 
00283   // TODO: This may not work with distributed matrices.
00284   template < class Matrix >
00285   void
00286   MatrixAdapter<Matrix>::do_getCrs(const Teuchos::ArrayView<scalar_t> nzval,
00287            const Teuchos::ArrayView<global_ordinal_t> colind,
00288            const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> rowptr,
00289            typename MatrixAdapter<Matrix>::global_size_t& nnz,
00290            const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > rowmap,
00291            EStorage_Ordering ordering,
00292            col_access ca) const
00293   {
00294     using Teuchos::Array;
00295     // get the ccs and transpose
00296 
00297     Array<scalar_t> nzval_tmp(nzval.size(), 0);
00298     Array<global_ordinal_t> rowind(colind.size(), 0);
00299     Array<global_size_t> colptr(this->getGlobalNumCols() + 1);
00300     this->getCcs(nzval_tmp(), rowind(), colptr(), nnz, rowmap, ordering);
00301     
00302     if( !nzval.is_null() && !colind.is_null() && !rowptr.is_null() )
00303       Util::transpose(nzval_tmp(), rowind(), colptr(), nzval, colind, rowptr);
00304   }
00305 
00306   template < class Matrix >
00307   void
00308   MatrixAdapter<Matrix>::help_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
00309              const Teuchos::ArrayView<global_ordinal_t> rowind,
00310              const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
00311              typename MatrixAdapter<Matrix>::global_size_t& nnz,
00312              const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
00313              EStorage_Ordering ordering,
00314              has_special_impl hsi) const
00315   {
00316     static_cast<const adapter_t*>(this)->getCcs_spec(nzval, rowind, colptr,
00317                  nnz, colmap, ordering);
00318   }
00319 
00320   template < class Matrix >
00321   void
00322   MatrixAdapter<Matrix>::help_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
00323              const Teuchos::ArrayView<global_ordinal_t> rowind,
00324              const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
00325              typename MatrixAdapter<Matrix>::global_size_t& nnz,
00326              const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
00327              EStorage_Ordering ordering,
00328              no_special_impl nsi) const
00329   {
00330     do_getCcs(nzval, rowind, colptr,
00331         nnz, colmap, ordering,
00332         typename adapter_t::major_access());
00333   }
00334 
00335   template < class Matrix >
00336   void
00337   MatrixAdapter<Matrix>::do_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
00338            const Teuchos::ArrayView<global_ordinal_t> rowind,
00339            const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
00340            typename MatrixAdapter<Matrix>::global_size_t& nnz,
00341            const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
00342            EStorage_Ordering ordering,
00343            row_access ra) const
00344   {
00345     using Teuchos::Array;
00346     // get the crs and transpose
00347 
00348     Array<scalar_t> nzval_tmp(nzval.size(), 0);
00349     Array<global_ordinal_t> colind(rowind.size(), 0);
00350     Array<global_size_t> rowptr(this->getGlobalNumRows() + 1);
00351     this->getCrs(nzval_tmp(), colind(), rowptr(), nnz, colmap, ordering);
00352 
00353     if( !nzval.is_null() && !rowind.is_null() && !colptr.is_null() )
00354       Util::transpose(nzval_tmp(), colind(), rowptr(), nzval, rowind, colptr);
00355   }
00356 
00357   template < class Matrix >
00358   void
00359   MatrixAdapter<Matrix>::do_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
00360            const Teuchos::ArrayView<global_ordinal_t> rowind,
00361            const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
00362            typename MatrixAdapter<Matrix>::global_size_t& nnz,
00363            const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
00364            EStorage_Ordering ordering,
00365            col_access ca) const
00366   {
00367     using Teuchos::RCP;
00368     using Teuchos::ArrayView;
00369     using Teuchos::OrdinalTraits;
00370     
00371     RCP<const type> get_mat;
00372     if( *colmap == *this->col_map_ ){
00373       // No need to redistribute
00374       get_mat = rcp(this,false); // non-owning
00375     } else {
00376       get_mat = get(colmap);
00377     }
00378 
00379     // If all is well and good, then colmap == cmap
00380     RCP<const Tpetra::Map<scalar_t,local_ordinal_t,global_ordinal_t> > cmap = get_mat->getColMap();
00381     TEUCHOS_ASSERT( *colmap == *cmap );
00382 
00383     ArrayView<global_ordinal_t> node_elements = cmap->getNodeElementList();
00384     if( node_elements.size() == 0 ) return; // no more contribution
00385     
00386     typename ArrayView<global_ordinal_t>::iterator col_it, col_end;
00387     col_end = node_elements.end();
00388 
00389     size_t colptr_ind = OrdinalTraits<size_t>::zero();
00390     global_ordinal_t colInd = OrdinalTraits<global_ordinal_t>::zero();
00391     for( col_it = node_elements.begin(); col_it != col_end; ++col_it ){
00392       colptr[colptr_ind++] = colInd;
00393       size_t colNNZ = getGlobalColNNZ(*col_it);
00394       size_t nnzRet = 0;
00395       ArrayView<global_ordinal_t> rowind_view = rowind.view(colInd,colNNZ);
00396       ArrayView<scalar_t> nzval_view = nzval.view(colInd,colNNZ);
00397       getGlobalColCopy(*col_it, rowind_view, nzval_view, nnzRet);
00398       
00399       // It was suggested that instead of sorting each row's indices
00400       // individually, that we instead do a double-transpose at the
00401       // end, which would also lead to the indices being sorted.
00402       if( ordering == SORTED_INDICES ){
00403   Tpetra::sort2(rowind_view.begin(), rowind_view.end(), nzval_view.begin());
00404       }
00405       
00406       TEUCHOS_TEST_FOR_EXCEPTION( colNNZ != nnzRet,
00407         std::runtime_error,
00408         "Number of values returned different from "
00409                           "number of values reported");
00410       colInd += colNNZ;
00411     }
00412     colptr[colptr_ind] = nnz = colInd;
00413   }
00414 
00415   
00416   // These will link to concrete implementations
00417   template < class Matrix >
00418   void
00419   MatrixAdapter<Matrix>::getGlobalRowCopy(global_ordinal_t row,
00420             const Teuchos::ArrayView<global_ordinal_t>& indices,
00421             const Teuchos::ArrayView<scalar_t>& vals,
00422             size_t& nnz) const
00423   {
00424     static_cast<const adapter_t*>(this)->getGlobalRowCopy_impl(row, indices, vals, nnz);
00425   }
00426   
00427   template < class Matrix >
00428   void
00429   MatrixAdapter<Matrix>::getGlobalColCopy(global_ordinal_t col,
00430             const Teuchos::ArrayView<global_ordinal_t>& indices,
00431             const Teuchos::ArrayView<scalar_t>& vals,
00432             size_t& nnz) const
00433   {
00434     static_cast<const adapter_t*>(this)->getGlobalColCopy_impl(col, indices, vals, nnz);
00435   }
00436 
00437   template < class Matrix >
00438   size_t
00439   MatrixAdapter<Matrix>::getMaxRowNNZ() const
00440   {
00441     return static_cast<const adapter_t*>(this)->getMaxRowNNZ_impl();
00442   }
00443 
00444   template < class Matrix >
00445   size_t
00446   MatrixAdapter<Matrix>::getMaxColNNZ() const
00447   {
00448     return static_cast<const adapter_t*>(this)->getMaxColNNZ_impl();
00449   }
00450     
00451   template < class Matrix >
00452   size_t
00453   MatrixAdapter<Matrix>::getGlobalRowNNZ(global_ordinal_t row) const
00454   {
00455     return static_cast<const adapter_t*>(this)->getGlobalRowNNZ_impl(row);
00456   }
00457 
00458   template < class Matrix >
00459   size_t
00460   MatrixAdapter<Matrix>::getLocalRowNNZ(local_ordinal_t row) const
00461   {
00462     return static_cast<const adapter_t*>(this)->getLocalRowNNZ_impl(row);
00463   }
00464 
00465   template < class Matrix >
00466   size_t
00467   MatrixAdapter<Matrix>::getGlobalColNNZ(global_ordinal_t col) const
00468   {
00469     return static_cast<const adapter_t*>(this)->getGlobalColNNZ_impl(col);
00470   }
00471 
00472   template < class Matrix >
00473   size_t
00474   MatrixAdapter<Matrix>::getLocalColNNZ(local_ordinal_t col) const
00475   {
00476     return static_cast<const adapter_t*>(this)->getLocalColNNZ_impl(col);
00477   }
00478 
00479   template < class Matrix >
00480   bool
00481   MatrixAdapter<Matrix>::isLocallyIndexed() const
00482   {
00483     return static_cast<const adapter_t*>(this)->isLocallyIndexed_impl();
00484   }
00485   
00486   template < class Matrix >
00487   bool
00488   MatrixAdapter<Matrix>::isGloballyIndexed() const
00489   {
00490     return static_cast<const adapter_t*>(this)->isGloballyIndexed_impl();
00491   }
00492   
00493   template < class Matrix >
00494   Teuchos::RCP<const MatrixAdapter<Matrix> >
00495   MatrixAdapter<Matrix>::get(const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > map) const
00496   {
00497     return static_cast<const adapter_t*>(this)->get_impl(map);
00498   }
00499 
00500 
00501   template <class Matrix>
00502   Teuchos::RCP<MatrixAdapter<Matrix> >
00503   createMatrixAdapter(Teuchos::RCP<Matrix> m){
00504     using Teuchos::rcp;
00505     using Teuchos::rcp_const_cast;
00506     
00507     if(m.is_null()) return Teuchos::null;
00508     return( rcp(new ConcreteMatrixAdapter<Matrix>(m)) );
00509   }
00510 
00511   template <class Matrix>
00512   Teuchos::RCP<const MatrixAdapter<Matrix> >
00513   createConstMatrixAdapter(Teuchos::RCP<const Matrix> m){
00514     using Teuchos::rcp;
00515     using Teuchos::rcp_const_cast;
00516     
00517     if(m.is_null()) return Teuchos::null;
00518     return( rcp(new ConcreteMatrixAdapter<Matrix>(rcp_const_cast<Matrix,const Matrix>(m))).getConst() );
00519   }
00520 
00521 } // end namespace Amesos2
00522 
00523 #endif  // AMESOS2_MATRIXADAPTER_DEF_HPP