Anasazi Version of the Day
TsqrAdaptor_Epetra_MultiVector.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef __TSQR_Trilinos_TsqrAdaptor_Epetra_MultiVector_hpp
00030 #define __TSQR_Trilinos_TsqrAdaptor_Epetra_MultiVector_hpp
00031 
00037 
00038 #include "TsqrTypeAdaptor_Epetra_Multivector.hpp"
00039 #include "TsqrCommFactory_Epetra.hpp"
00040 #include <stdexcept>
00041 #include <sstream>
00042 
00045 
00046 namespace TSQR {
00047   namespace Trilinos {
00048 
00049     class TsqrEpetraAdaptor : 
00050       public TsqrAdaptor< double, int, int, Epetra_MultiVector >
00051     {
00052     public:
00053       // C++ doesn't inherit typedefs from base to derived class, when
00054       // both the base and the derived classes are templated.  To
00055       // avoid code duplication, we pull in typedefs from base_type.
00056       typedef TsqrAdaptor< double, int, int, Epetra_MultiVector > base_type;
00057       typedef base_type::comm_ptr comm_ptr;
00058       typedef base_type::multivector_type multivector_type;
00059       typedef base_type::scalar_type scalar_type;
00060       typedef base_type::local_ordinal_type local_ordinal_type;
00061       typedef typename base_type::scalar_messenger_ptr scalar_messenger_ptr;
00062       typedef typename base_type::ordinal_messenger_ptr ordinal_messenger_ptr;
00063 
00077       TsqrEpetraAdaptor (const multivector_type& mv,
00078        const Teuchos::ParameterList& plist)
00079       {
00080   init (mv, plist);
00081       }
00082 
00083     private:
00084       virtual void
00085       fetchDims (const multivector_type& A,
00086      local_ordinal_type& nrowsLocal,
00087      local_ordinal_type& ncols,
00088      local_ordinal_type& LDA) const
00089       {
00090   nrowsLocal = A.MyLength();
00091   ncols = A.NumVectors();
00092   if (nrowsLocal < ncols)
00093     {
00094       std::ostringstream os;
00095       os << "TSQR: The local component of the input matrix has fewer row"
00096         "s (" << nrowsLocal << ") than columns (" << ncols << ").  TSQR "
00097         "does not support this case.";
00098       throw std::runtime_error (os.str());
00099     }
00100   if (! A.ConstantStride())
00101     {
00102       std::ostringstream os;
00103       os << "TSQR does not support Epetra_MultiVector inputs that do not"
00104         " have constant stride.";
00105       throw std::runtime_error (os.str());
00106     }
00107   LDA = A.Stride();
00108       }
00109 
00110       Teuchos::ArrayView< scalar_type >::size_type 
00111       fetchArrayLength (const multivector_type& A) const
00112       {
00113   typedef Teuchos::ArrayView< scalar_type >::size_type size_type;
00114 
00115   // Compute length (nelts) of the A.Values() array.  This only
00116   // makes sense if A has constant stride.  We convert from int
00117   // (the local_ordinal_type) to size_type, hoping that the
00118   // latter is no smaller than local_ordinal_type.  This gives
00119   // us a better chance that LDA*ncols fits in a size_type.  We
00120   // check for overflow both when converting each of LDA and
00121   // ncols from local_ordinal_type to size_type, and when
00122   // computing LDA*ncols (in size_type arithmetic).
00123 
00124   // Fetch dimensions of local part of A, as local_ordinal_type.
00125   local_ordinal_type nrowsLocal_LO, ncols_LO, LDA_LO;
00126   fetchDims (A, nrowsLocal_LO, ncols_LO, LDA_LO);
00127 
00128   // Convert dimensions from local_ordinal_type to size_type,
00129   // ensuring that the conversions don't overflow.
00130   const size_type LDA = static_cast< size_type > (LDA_LO);
00131   if (static_cast< local_ordinal_type > (LDA) != LDA_LO)
00132     throw std::runtime_error ("TSQR::Trilinos::TsqrEpetraAdaptor: "
00133             "Leading dimension of local part of "
00134             "Epetra_MultiVector overflowed on co"
00135             "nversion from local_ordinal_type to"
00136             " ArrayView::size_type");
00137   const size_type ncols = static_cast< size_type > (ncols_LO);
00138   if (static_cast< local_ordinal_type > (ncols) != ncols_LO)
00139     throw std::runtime_error ("TSQR::Trilinos::TsqrEpetraAdaptor: "
00140             "Number of columns of the given "
00141             "Epetra_MultiVector overflowed on co"
00142             "nversion from local_ordinal_type to"
00143             " ArrayView::size_type");
00144   // Make sure that the length of A.Values(), which is the
00145   // product of the leading dimension and the number of columns
00146   // of A, fits in a size_type.
00147   const size_type nelts = LDA * ncols;
00148   if (nelts / LDA != ncols)
00149     throw std::runtime_error ("TSQR::Trilinos::TsqrEpetraAdaptor: "
00150             " Length of A.Values() does not fit "
00151             "in an ArrayView::size_type");
00152   return nelts;
00153       }
00154 
00155       virtual Teuchos::ArrayRCP< scalar_type > 
00156       fetchNonConstView (multivector_type& A) const
00157       {
00158   using Teuchos::arcpFromArrayView;
00159   using Teuchos::arrayView;
00160   typedef Teuchos::ArrayView< scalar_type >::size_type size_type;
00161 
00162   const size_type nelts = fetchArrayLength (A);
00163   // The returned ArrayRCP does NOT own A.Values().
00164   return arcpFromArrayView (arrayView (A.Values(), nelts));
00165       }
00166 
00167       virtual Teuchos::ArrayRCP< const scalar_type > 
00168       fetchConstView (const multivector_type& A) const
00169       {
00170   using Teuchos::arcpFromArrayView;
00171   using Teuchos::arrayView;
00172   using Teuchos::ArrayView;
00173   typedef ArrayView< scalar_type >::size_type size_type;
00174 
00175   const size_type nelts = fetchArrayLength (A);
00176   const scalar_type* A_ptr = A.Values();
00177   ArrayView< const scalar_type > A_view = arrayView (A_ptr, nelts);
00178 
00179   // The returned ArrayRCP does NOT own A.Values().
00180   return arcpFromArrayView (A_view);
00181       }
00182 
00183       virtual void
00184       fetchMessengers (const multivector_type& mv,
00185            scalar_messenger_ptr& pScalarMessenger,
00186            ordinal_messenger_ptr& pOrdinalMessenger) const
00187       {
00188   typedef EpetraCommFactory comm_factory_type;
00189   // mv.Comm() returns a "const Epetra_Comm&".  rcpFromRef()
00190   // makes a non-owning (weak) RCP out of it.  This requires
00191   // that the mv's Epetra_Comm object not fall out of scope,
00192   // which it shouldn't as long as we are computing with
00193   // multivectors distributed according to that communicator.
00194   comm_ptr pComm = Teuchos::rcpFromRef (mv.Comm());
00195   comm_factory_type factory;
00196   factory.makeMessengers (pComm, pScalarMessenger, pOrdinalMessenger);
00197       }
00198     };
00199 
00200   } // namespace Trilinos
00201 } // namespace TSQR
00202 
00203 #endif // __TSQR_Trilinos_TsqrAdaptor_Epetra_MultiVector_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends