Tpetra Matrix/Vector Services Version of the Day
Tpetra_TsqrAdaptor.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) 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 __Tpetra_TsqrAdaptor_hpp
00043 #define __Tpetra_TsqrAdaptor_hpp
00044 
00048 
00049 #include <Tpetra_ConfigDefs.hpp> // HAVE_TPETRA_TSQR, etc.
00050 
00051 #ifdef HAVE_TPETRA_TSQR
00052 #  include <Tsqr_NodeTsqrFactory.hpp> // create intranode TSQR object
00053 #  include <Tsqr.hpp> // full (internode + intranode) TSQR
00054 #  include <Tsqr_DistTsqr.hpp> // internode TSQR
00055 // Subclass of TSQR::MessengerBase, implemented using Teuchos
00056 // communicator template helper functions
00057 #  include <Tsqr_TeuchosMessenger.hpp>
00058 #  include <Tpetra_MultiVector.hpp>
00059 #  include <Teuchos_ParameterListAcceptorDefaultBase.hpp>
00060 #  include <stdexcept>
00061 
00062 
00063 namespace Tpetra {
00064 
00087   template<class MV>
00088   class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
00089   public:
00090     typedef typename MV::scalar_type scalar_type;
00091     typedef typename MV::local_ordinal_type ordinal_type;
00092     typedef typename MV::node_type node_type;
00093     typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
00094     typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
00095 
00096   private:
00097     //typedef TSQR::MatView<ordinal_type, scalar_type> matview_type;
00098     typedef TSQR::NodeTsqrFactory<node_type, scalar_type, ordinal_type> node_tsqr_factory_type;
00099     typedef typename node_tsqr_factory_type::node_tsqr_type node_tsqr_type;
00100     typedef TSQR::DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
00101     typedef TSQR::Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;
00102 
00103   public:
00110     TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
00111       nodeTsqr_ (new node_tsqr_type),
00112       distTsqr_ (new dist_tsqr_type),
00113       tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
00114       ready_ (false)
00115     {
00116       setParameterList (plist);
00117     }
00118 
00120     TsqrAdaptor () :
00121       nodeTsqr_ (new node_tsqr_type),
00122       distTsqr_ (new dist_tsqr_type),
00123       tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
00124       ready_ (false)
00125     {
00126       setParameterList (Teuchos::null);
00127     }
00128 
00130     Teuchos::RCP<const Teuchos::ParameterList>
00131     getValidParameters () const
00132     {
00133       using Teuchos::RCP;
00134       using Teuchos::rcp;
00135       using Teuchos::ParameterList;
00136       using Teuchos::parameterList;
00137 
00138       if (defaultParams_.is_null()) {
00139         RCP<ParameterList> params = parameterList ("TSQR implementation");
00140         params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
00141         params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
00142         defaultParams_ = params;
00143       }
00144       return defaultParams_;
00145     }
00146 
00172     void
00173     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00174     {
00175       using Teuchos::ParameterList;
00176       using Teuchos::parameterList;
00177       using Teuchos::RCP;
00178       using Teuchos::sublist;
00179 
00180       RCP<ParameterList> params = plist.is_null() ?
00181         parameterList (*getValidParameters ()) : plist;
00182       nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
00183       distTsqr_->setParameterList (sublist (params, "DistTsqr"));
00184 
00185       this->setMyParamList (params);
00186     }
00187 
00209     void
00210     factorExplicit (MV& A,
00211                     MV& Q,
00212                     dense_matrix_type& R,
00213                     const bool forceNonnegativeDiagonal=false)
00214     {
00215       typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
00216 
00217       prepareTsqr (Q); // Finish initializing TSQR.
00218       KMV A_view = getNonConstView (A);
00219       KMV Q_view = getNonConstView (Q);
00220       tsqr_->factorExplicit (A_view, Q_view, R, false,
00221                              forceNonnegativeDiagonal);
00222     }
00223 
00254     int
00255     revealRank (MV& Q,
00256                 dense_matrix_type& R,
00257                 const magnitude_type& tol)
00258     {
00259       typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
00260 
00261       prepareTsqr (Q); // Finish initializing TSQR.
00262 
00263       // FIXME (mfh 18 Oct 2010) Check Teuchos::Comm<int> object in Q
00264       // to make sure it is the same communicator as the one we are
00265       // using in our dist_tsqr_type implementation.
00266       KMV Q_view = getNonConstView (Q);
00267       return tsqr_->revealRank (Q_view, R, tol, false);
00268     }
00269 
00270   private:
00272     Teuchos::RCP<node_tsqr_type> nodeTsqr_;
00273 
00275     Teuchos::RCP<dist_tsqr_type> distTsqr_;
00276 
00278     Teuchos::RCP<tsqr_type> tsqr_;
00279 
00281     mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00282 
00284     bool ready_;
00285 
00306     void
00307     prepareTsqr (const MV& mv)
00308     {
00309       if (! ready_) {
00310         prepareDistTsqr (mv);
00311         prepareNodeTsqr (mv);
00312         ready_ = true;
00313       }
00314     }
00315 
00319     void
00320     prepareNodeTsqr (const MV& mv)
00321     {
00322       node_tsqr_factory_type::prepareNodeTsqr (nodeTsqr_, mv.getMap()->getNode());
00323     }
00324 
00331     void
00332     prepareDistTsqr (const MV& mv)
00333     {
00334       using Teuchos::RCP;
00335       using Teuchos::rcp_implicit_cast;
00336       typedef TSQR::TeuchosMessenger<scalar_type> mess_type;
00337       typedef TSQR::MessengerBase<scalar_type> base_mess_type;
00338 
00339       RCP<const Teuchos::Comm<int> > comm = mv.getMap()->getComm();
00340       RCP<mess_type> mess (new mess_type (comm));
00341       RCP<base_mess_type> messBase = rcp_implicit_cast<base_mess_type> (mess);
00342       distTsqr_->init (messBase);
00343     }
00344 
00357     static KokkosClassic::MultiVector<scalar_type, node_type>
00358     getNonConstView (MV& A)
00359     {
00360       // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
00361       // storage of A uses nonconstant stride internally.  We would
00362       // have to copy and pack into a matrix with constant stride, and
00363       // then unpack on exit.  For now we choose just to raise an
00364       // exception.
00365       TEUCHOS_TEST_FOR_EXCEPTION(! A.isConstantStride(), std::invalid_argument,
00366                                  "TSQR does not currently support Tpetra::MultiVector "
00367                                  "inputs that do not have constant stride.");
00368       return A.getLocalMVNonConst();
00369     }
00370   };
00371 
00372 } // namespace Tpetra
00373 
00374 #endif // HAVE_TPETRA_TSQR
00375 
00376 #endif // __Tpetra_TsqrAdaptor_hpp
00377 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines