Tpetra Matrix/Vector Services Version of the Day
Epetra_TsqrAdaptor.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright 2010 Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef __Epetra_TsqrAdaptor_hpp
00043 #define __Epetra_TsqrAdaptor_hpp
00044 
00058 
00059 #include <Tpetra_ConfigDefs.hpp>
00060 
00061 #if defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
00062 
00063 #include <Kokkos_DefaultNode.hpp> // Include minimal Kokkos Node types
00064 #include <Tsqr_NodeTsqrFactory.hpp> // create intranode TSQR object
00065 #include <Tsqr.hpp> // full (internode + intranode) TSQR
00066 #include <Tsqr_DistTsqr.hpp> // internode TSQR
00067 
00068 #include <Epetra_Comm.h>
00069 // Subclass of TSQR::MessengerBase, implemented using Teuchos
00070 // communicator template helper functions
00071 #include <Epetra_TsqrMessenger.hpp>
00072 #include <Epetra_MultiVector.h>
00073 
00074 #include <stdexcept>
00075 
00078 
00079 namespace Epetra {
00080 
00101   class TsqrAdaptor {
00102   public:
00103     typedef Epetra_MultiVector MV;
00104 
00111     typedef double scalar_type;
00112 
00119     typedef int ordinal_type;
00120 
00131     typedef Kokkos::SerialNode node_type;
00132 
00141     typedef Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > dense_matrix_type;
00142 
00148     typedef double magnitude_type;
00149 
00150   private:
00151     typedef TSQR::MatView< ordinal_type, scalar_type > matview_type;
00152     typedef TSQR::NodeTsqrFactory< node_type, scalar_type, ordinal_type > node_tsqr_factory_type;
00153     // Don't need a "typename" here, because there are no template
00154     // parameters involved in the type definition.
00155     typedef node_tsqr_factory_type::node_tsqr_type node_tsqr_type;
00156     typedef TSQR::DistTsqr< ordinal_type, scalar_type > dist_tsqr_type;
00157     typedef TSQR::Tsqr< ordinal_type, scalar_type, node_tsqr_type > tsqr_type;
00158 
00159   public:
00167     static Teuchos::RCP<const Teuchos::ParameterList>
00168     getDefaultParameters ()
00169     {
00170       // For now, only the intranode part of TSQR accepts parameters.
00171       return node_tsqr_factory_type::getDefaultParameters ();
00172     }
00173 
00193     TsqrAdaptor (const MV& mv,
00194      const Teuchos::RCP<const Teuchos::ParameterList>& plist) :
00195       pTsqr_ (new tsqr_type (makeNodeTsqr (plist), makeDistTsqr (mv)))
00196     {
00197       // This particular adaptor wants an instance of the specific
00198       // Kokkos Node type.  All of the Node constructors take
00199       // parameters.  SerialNode doesn't currently read any of them,
00200       // but we still have to pass in an empty parameter list.
00201       Teuchos::ParameterList emptyParams;
00202       pNode_ = Teuchos::rcp (new Kokkos::SerialNode (emptyParams));
00203     }
00204 
00207     void
00208     factorExplicit (MV& A,
00209         MV& Q,
00210         dense_matrix_type& R)
00211     {
00212       typedef Kokkos::MultiVector< scalar_type, node_type > KMV;
00213 
00214       // FIXME (mfh 25 Oct 2010) Check Epetra_Comm objects in A and Q
00215       // to make sure they are the same communicator as the one we are
00216       // using in our dist_tsqr_type implementation.
00217       //
00218       // mfh 11 Jan 2011: Actually we should also check whether the
00219       // Epetra_Map objects are compatible.  Fixing this at the TSQR
00220       // implementation level likely will require reimplementing TSQR
00221       // with Map awareness.  This will be troublesome if we want to
00222       // maintain compatibility with both Epetra and Tpetra.
00223       // Alternately, we can check for compatibility of maps at the
00224       // adaptor or even the (Mat)OrthoManager levels.
00225       KMV A_view = getNonConstView (A);
00226       KMV Q_view = getNonConstView (Q);
00227       pTsqr_->factorExplicit (A_view, Q_view, R, false);
00228     }
00229 
00261     int
00262     revealRank (MV& Q,
00263     dense_matrix_type& R,
00264     const magnitude_type& tol)
00265     {
00266       typedef Kokkos::MultiVector< scalar_type, node_type > KMV;
00267 
00268       // FIXME (mfh 25 Oct 2010) Check Epetra_Comm object in Q to make
00269       // sure it is the same communicator as the one we are using in
00270       // our dist_tsqr_type implementation.
00271 
00272       KMV Q_view = getNonConstView (Q);
00273       return pTsqr_->revealRank (Q_view, R, tol, false);
00274     }
00275 
00276   private:
00278     Teuchos::RCP< tsqr_type > pTsqr_;
00279 
00285     Teuchos::RCP< Kokkos::SerialNode > pNode_;
00286 
00304     Kokkos::MultiVector< scalar_type, node_type >
00305     getNonConstView (MV& A)
00306     {
00307       if (! A.ConstantStride())
00308   {
00309     // FIXME (mfh 25 Oct 2010) Storage of A uses nonconstant
00310     // stride internally, but that doesn't necessarily mean we
00311     // can't run TSQR.  We would have to copy and pack into a
00312     // matrix with constant stride, and then unpack on exit.
00313     std::ostringstream os;
00314     os << "TSQR does not currently support Epetra_MultiVector "
00315       "inputs that do not have constant stride.";
00316     throw std::runtime_error (os.str());
00317   }
00318 
00319       const int numRows = A.MyLength();
00320       const int numCols = A.NumVectors();
00321       const int stride  = A.Stride();
00322       // A_ptr does _not_ own the data.  TSQR only operates within the
00323       // scope of the multivector objects on which it operates, so it
00324       // doesn't need ownership of the data.
00325       Teuchos::ArrayRCP< double > A_ptr (A.Values(), 0, numRows*stride, false);
00326 
00327       typedef Kokkos::MultiVector< scalar_type, node_type > KMV;
00328       KMV A_kmv (pNode_);
00329       A_kmv.initializeValues (numRows, numCols, A_ptr, stride);
00330       return A_kmv;
00331     }
00332 
00334     static Teuchos::RCP< dist_tsqr_type > 
00335     makeDistTsqr (const MV& mv)
00336     {
00337       using Teuchos::RCP;
00338       typedef TSQR::MessengerBase< scalar_type > base_mess_type;
00339       using TSQR::Epetra::makeTsqrMessenger;
00340 
00341       // pComm is a nonowning pointer to the reference.
00342       RCP< const Epetra_Comm > pComm = Teuchos::rcpFromRef (mv.Comm());
00343       RCP< base_mess_type > pMessBase = makeTsqrMessenger< scalar_type >(pComm);
00344       RCP< dist_tsqr_type > pDistTsqr (new dist_tsqr_type (pMessBase));
00345       return pDistTsqr;
00346     }
00347 
00349     static Teuchos::RCP< node_tsqr_type >
00350     makeNodeTsqr (const Teuchos::RCP<const Teuchos::ParameterList>& plist)
00351     {
00352       return node_tsqr_factory_type::makeNodeTsqr (plist);
00353     }
00354   };
00355 
00356 } // namespace Epetra
00357 
00358 #endif // defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
00359 
00360 #endif // __Epetra_TsqrAdaptor_hpp
00361 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines