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 
00224     int
00225     revealRank (MV& Q,
00226     dense_matrix_type& R,
00227     const magnitude_type& tol)
00228     {
00229       typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00230 
00231       prepareTsqr (Q); // Finish initializing TSQR.      
00232 
00233       // FIXME (mfh 18 Oct 2010) Check Teuchos::Comm<int> object in Q
00234       // to make sure it is the same communicator as the one we are
00235       // using in our dist_tsqr_type implementation.
00236       KMV Q_view = getNonConstView (Q);
00237       return tsqr_->revealRank (Q_view, R, tol, false);
00238     }
00239 
00240   private:
00242     Teuchos::RCP<node_tsqr_type> nodeTsqr_;
00243     
00245     Teuchos::RCP<dist_tsqr_type> distTsqr_;
00246 
00248     Teuchos::RCP<tsqr_type> tsqr_;
00249 
00251     mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00252 
00254     bool ready_;
00255 
00276     void 
00277     prepareTsqr (const MV& mv) 
00278     {
00279       if (! ready_) {
00280   prepareDistTsqr (mv);
00281   prepareNodeTsqr (mv);
00282   ready_ = true;
00283       }
00284     }
00285 
00289     void
00290     prepareNodeTsqr (const MV& mv)
00291     {
00292       node_tsqr_factory_type::prepareNodeTsqr (nodeTsqr_, mv.getMap()->getNode());
00293     }
00294 
00301     void
00302     prepareDistTsqr (const MV& mv)
00303     {
00304       using Teuchos::RCP;
00305       using Teuchos::rcp_implicit_cast;
00306       typedef TSQR::TeuchosMessenger<scalar_type> mess_type;
00307       typedef TSQR::MessengerBase<scalar_type> base_mess_type;
00308 
00309       RCP<const Teuchos::Comm<int> > comm = mv.getMap()->getComm();
00310       RCP<mess_type> mess (new mess_type (comm));
00311       RCP<base_mess_type> messBase = rcp_implicit_cast<base_mess_type> (mess);
00312       distTsqr_->init (messBase);
00313     }
00314 
00327     static Kokkos::MultiVector<scalar_type, node_type>
00328     getNonConstView (MV& A)
00329     {
00330       // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
00331       // storage of A uses nonconstant stride internally.  We would
00332       // have to copy and pack into a matrix with constant stride, and
00333       // then unpack on exit.  For now we choose just to raise an
00334       // exception.
00335       TEUCHOS_TEST_FOR_EXCEPTION(! A.isConstantStride(), std::invalid_argument,
00336          "TSQR does not currently support Tpetra::MultiVector "
00337          "inputs that do not have constant stride.");
00338       return A.getLocalMVNonConst();
00339     }
00340   };
00341 
00342 } // namespace Tpetra
00343 
00344 #endif // HAVE_TPETRA_TSQR
00345 
00346 #endif // __Tpetra_TsqrAdaptor_hpp
00347 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines