Tpetra Matrix/Vector Services Version of the Day
Tpetra_TsqrAdaptor.hpp
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 
00045 #include <Tpetra_ConfigDefs.hpp> // HAVE_TPETRA_TSQR, etc.
00046 
00047 #ifdef HAVE_TPETRA_TSQR
00048 #  include <Tsqr_NodeTsqrFactory.hpp> // create intranode TSQR object
00049 #  include <Tsqr.hpp> // full (internode + intranode) TSQR
00050 #  include <Tsqr_DistTsqr.hpp> // internode TSQR
00051 // Subclass of TSQR::MessengerBase, implemented using Teuchos
00052 // communicator template helper functions
00053 #  include <Tsqr_TeuchosMessenger.hpp> 
00054 #  include <Tpetra_MultiVector.hpp>
00055 #  include <Teuchos_ParameterListAcceptorDefaultBase.hpp>
00056 #  include <stdexcept>
00057 
00058 
00059 namespace Tpetra {
00060 
00083   template<class MV>
00084   class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
00085   public:
00086     typedef typename MV::scalar_type scalar_type;
00087     typedef typename MV::local_ordinal_type ordinal_type;
00088     typedef typename MV::node_type node_type;
00089     typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
00090     typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
00091 
00092   private:
00093     typedef TSQR::MatView<ordinal_type, scalar_type> matview_type;
00094     typedef TSQR::NodeTsqrFactory<node_type, scalar_type, ordinal_type> node_tsqr_factory_type;
00095     typedef typename node_tsqr_factory_type::node_tsqr_type node_tsqr_type;
00096     typedef TSQR::DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
00097     typedef TSQR::Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;
00098 
00099   public:
00106     TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
00107       nodeTsqr_ (new node_tsqr_type),
00108       distTsqr_ (new dist_tsqr_type),
00109       tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
00110       ready_ (false)
00111     {
00112       setParameterList (plist);
00113     }
00114 
00116     TsqrAdaptor () : 
00117       nodeTsqr_ (new node_tsqr_type),
00118       distTsqr_ (new dist_tsqr_type),
00119       tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
00120       ready_ (false)
00121     {
00122       setParameterList (Teuchos::null);
00123     }
00124 
00125     Teuchos::RCP<const Teuchos::ParameterList>
00126     getValidParameters () const
00127     {
00128       using Teuchos::RCP;
00129       using Teuchos::rcp;
00130       using Teuchos::ParameterList;
00131       using Teuchos::parameterList;
00132 
00133       if (defaultParams_.is_null()) {
00134         RCP<ParameterList> params = parameterList ("TSQR implementation");
00135         params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
00136         params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
00137         defaultParams_ = params;
00138       }
00139       return defaultParams_;
00140     }
00141 
00142     void 
00143     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00144     {
00145       using Teuchos::ParameterList;
00146       using Teuchos::parameterList;
00147       using Teuchos::RCP;
00148       using Teuchos::sublist;
00149 
00150       RCP<ParameterList> params = plist.is_null() ? 
00151         parameterList (*getValidParameters ()) : plist;
00152       nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
00153       distTsqr_->setParameterList (sublist (params, "DistTsqr"));
00154 
00155       this->setMyParamList (params);
00156     }
00157 
00179     void
00180     factorExplicit (MV& A,
00181                     MV& Q,
00182                     dense_matrix_type& R,
00183                     const bool forceNonnegativeDiagonal=false)
00184     {
00185       typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00186 
00187       prepareTsqr (Q); // Finish initializing TSQR.
00188       KMV A_view = getNonConstView (A);
00189       KMV Q_view = getNonConstView (Q);
00190       tsqr_->factorExplicit (A_view, Q_view, R, false, 
00191                              forceNonnegativeDiagonal);
00192     }
00193 
00225     int
00226     revealRank (MV& Q,
00227                 dense_matrix_type& R,
00228                 const magnitude_type& tol)
00229     {
00230       typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00231 
00232       prepareTsqr (Q); // Finish initializing TSQR.      
00233 
00234       // FIXME (mfh 18 Oct 2010) Check Teuchos::Comm<int> object in Q
00235       // to make sure it is the same communicator as the one we are
00236       // using in our dist_tsqr_type implementation.
00237       KMV Q_view = getNonConstView (Q);
00238       return tsqr_->revealRank (Q_view, R, tol, false);
00239     }
00240 
00241   private:
00243     Teuchos::RCP<node_tsqr_type> nodeTsqr_;
00244     
00246     Teuchos::RCP<dist_tsqr_type> distTsqr_;
00247 
00249     Teuchos::RCP<tsqr_type> tsqr_;
00250 
00252     mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00253 
00255     bool ready_;
00256 
00277     void 
00278     prepareTsqr (const MV& mv) 
00279     {
00280       if (! ready_) {
00281         prepareDistTsqr (mv);
00282         prepareNodeTsqr (mv);
00283         ready_ = true;
00284       }
00285     }
00286 
00290     void
00291     prepareNodeTsqr (const MV& mv)
00292     {
00293       node_tsqr_factory_type::prepareNodeTsqr (nodeTsqr_, mv.getMap()->getNode());
00294     }
00295 
00302     void
00303     prepareDistTsqr (const MV& mv)
00304     {
00305       using Teuchos::RCP;
00306       using Teuchos::rcp_implicit_cast;
00307       typedef TSQR::TeuchosMessenger<scalar_type> mess_type;
00308       typedef TSQR::MessengerBase<scalar_type> base_mess_type;
00309 
00310       RCP<const Teuchos::Comm<int> > comm = mv.getMap()->getComm();
00311       RCP<mess_type> mess (new mess_type (comm));
00312       RCP<base_mess_type> messBase = rcp_implicit_cast<base_mess_type> (mess);
00313       distTsqr_->init (messBase);
00314     }
00315 
00328     static Kokkos::MultiVector<scalar_type, node_type>
00329     getNonConstView (MV& A)
00330     {
00331       // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
00332       // storage of A uses nonconstant stride internally.  We would
00333       // have to copy and pack into a matrix with constant stride, and
00334       // then unpack on exit.  For now we choose just to raise an
00335       // exception.
00336       TEUCHOS_TEST_FOR_EXCEPTION(! A.isConstantStride(), std::invalid_argument,
00337                                  "TSQR does not currently support Tpetra::MultiVector "
00338                                  "inputs that do not have constant stride.");
00339       return A.getLocalMVNonConst();
00340     }
00341   };
00342 
00343 } // namespace Tpetra
00344 
00345 #endif // HAVE_TPETRA_TSQR
00346 
00347 #endif // __Tpetra_TsqrAdaptor_hpp
00348 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines