Anasazi Version of the Day
TsqrRandomizer_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_Randomizer_Epetra_MultiVector_hpp
00030 #define __TSQR_Trilinos_Randomizer_Epetra_MultiVector_hpp
00031 
00035 
00036 #include "TsqrTypeAdaptor_Epetra_Multivector.hpp"
00037 #include "TsqrCommFactory_Epetra.hpp"
00038 #include "Teuchos_ArrayView.hpp"
00039 #include <limits>
00040 #include <stdexcept>
00041 #include <sstream>
00042 
00045 
00046 namespace TSQR {
00047   namespace Trilinos {
00048 
00049     template< class Gen >
00050     class EpetraRandomizer :
00051       public Randomizer< double, int, int, Epetra_MultiVector, Gen >
00052     {
00053     public:
00054       // The C++ compiler needs help inheriting typedefs, when both
00055       // the base class and the derived class are templated.
00056       typedef Randomizer< double, int, int, Epetra_MultiVector, Gen > base_type;
00057       typedef typename base_type::multivector_type multivector_type;
00058       typedef typename base_type::local_ordinal_type local_ordinal_type;
00059       typedef typename base_type::scalar_type scalar_type;
00060 
00061       typedef typename base_type::normalgen_ptr normalgen_ptr;
00062       typedef typename base_type::ordinal_messenger_ptr ordinal_messenger_ptr;
00063       typedef typename base_type::scalar_messenger_ptr scalar_messenger_ptr;
00064 
00074       EpetraRandomizer (const multivector_type& mv,
00075       const normalgen_ptr& pGen)
00076       {
00077   init (mv, pGen);
00078       }
00079 
00080     private:
00081       // mfh 14 Jul 2010: For who knows what reason, this refuses to
00082       // compile if you replace "double" with "scalar_type".  And no,
00083       // adding "typename" doesn't work either.
00084       typedef Teuchos::ArrayView< double >::size_type size_type;
00085 
00086       virtual void
00087       fetchDims (const multivector_type& A,
00088      local_ordinal_type& nrowsLocal,
00089      local_ordinal_type& ncols,
00090      local_ordinal_type& LDA) const
00091       {
00092   nrowsLocal = A.MyLength();
00093   ncols = A.NumVectors();
00094   if (nrowsLocal < ncols)
00095     {
00096       std::ostringstream os;
00097       os << "TSQR: The local component of the input matrix has fewer row"
00098         "s (" << nrowsLocal << ") than columns (" << ncols << ").  TSQR "
00099         "does not support this case.";
00100       throw std::runtime_error (os.str());
00101     }
00102   if (! A.ConstantStride())
00103     {
00104       std::ostringstream os;
00105       os << "TSQR does not support Epetra_MultiVector inputs that do not"
00106         " have constant stride.";
00107       throw std::runtime_error (os.str());
00108     }
00109   LDA = A.Stride();
00110       }
00111 
00112       size_type
00113       fetchArrayLength (const multivector_type& A) const
00114       {
00115   //typedef Teuchos::ArrayView< scalar_type >::size_type size_type;
00116 
00117   // Compute length (nelts) of the A.Values() array.  This only
00118   // makes sense if A has constant stride.  We convert from int
00119   // (the local_ordinal_type) to size_type, hoping that the
00120   // latter is no smaller than local_ordinal_type.  This gives
00121   // us a better chance that LDA*ncols fits in a size_type.  We
00122   // check for overflow both when converting each of LDA and
00123   // ncols from local_ordinal_type to size_type, and when
00124   // computing LDA*ncols (in size_type arithmetic).
00125 
00126   // Fetch dimensions of local part of A, as local_ordinal_type.
00127   local_ordinal_type nrowsLocal_LO, ncols_LO, LDA_LO;
00128   fetchDims (A, nrowsLocal_LO, ncols_LO, LDA_LO);
00129 
00130   // Convert dimensions from local_ordinal_type to size_type,
00131   // ensuring that the conversions don't overflow.
00132   const size_type LDA = static_cast< size_type > (LDA_LO);
00133   if (static_cast< local_ordinal_type > (LDA) != LDA_LO)
00134     throw std::runtime_error ("TSQR::Trilinos::EpetraRandomizer: "
00135             "Leading dimension of local part of "
00136             "Epetra_MultiVector overflowed on co"
00137             "nversion from local_ordinal_type to"
00138             " ArrayView::size_type");
00139   const size_type ncols = static_cast< size_type > (ncols_LO);
00140   if (static_cast< local_ordinal_type > (ncols) != ncols_LO)
00141     throw std::runtime_error ("TSQR::Trilinos::EpetraRandomizer: "
00142             "Number of columns of the given "
00143             "Epetra_MultiVector overflowed on co"
00144             "nversion from local_ordinal_type to"
00145             " ArrayView::size_type");
00146   // Make sure that the length of A.Values(), which is the
00147   // product of the leading dimension and the number of columns
00148   // of A, fits in a size_type.
00149   const size_type nelts = LDA * ncols;
00150   if (nelts / LDA != ncols)
00151     throw std::runtime_error ("TSQR::Trilinos::EpetraRandomizer: "
00152             " Length of A.Values() does not fit "
00153             "in an ArrayView::size_type");
00154   return nelts;
00155       }
00156 
00157       virtual Teuchos::ArrayRCP< scalar_type > 
00158       fetchNonConstView (multivector_type& A) const
00159       {
00160   using Teuchos::arcpFromArrayView;
00161   using Teuchos::arrayView;
00162 
00163   const size_type nelts = fetchArrayLength (A);
00164   // The returned ArrayRCP does NOT own A.Values().
00165   return arcpFromArrayView (arrayView (A.Values(), nelts));
00166       }
00167 
00168       virtual void
00169       fetchMessengers (const multivector_type& mv,
00170            scalar_messenger_ptr& pScalarMessenger,
00171            ordinal_messenger_ptr& pOrdinalMessenger) const
00172       {
00173   typedef EpetraCommFactory comm_factory_type;
00174   // mv.Comm() returns a "const Epetra_Comm&".  rcpFromRef()
00175   // makes a non-owning (weak) RCP out of it.  This requires
00176   // that the mv's Epetra_Comm object not fall out of scope,
00177   // which it shouldn't as long as we are computing with
00178   // multivectors distributed according to that communicator.
00179   comm_ptr pComm = Teuchos::rcpFromRef (mv.Comm());
00180   comm_factory_type factory;
00181   factory.makeMessengers (pComm, pScalarMessenger, pOrdinalMessenger);
00182       }
00183     };
00184 
00185   } // namespace Trilinos
00186 } // namespace TSQR
00187 
00188 #endif // __TSQR_Trilinos_Randomizer_Epetra_MultiVector_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends