Tpetra Matrix/Vector Services Version of the Day
Epetra_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 __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 #include <Epetra_Comm.h>
00068 // Subclass of TSQR::MessengerBase, implemented using Teuchos
00069 // communicator template helper functions
00070 #include <Epetra_TsqrMessenger.hpp>
00071 #include <Epetra_MultiVector.h>
00072 #include <Teuchos_ParameterListAcceptorDefaultBase.hpp>
00073 #include <stdexcept>
00074 
00075 
00076 namespace Epetra {
00077 
00102   class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
00103   public:
00104     typedef Epetra_MultiVector MV;
00105 
00112     typedef double scalar_type;
00113 
00120     typedef int ordinal_type;
00121 
00132     typedef Kokkos::SerialNode node_type;
00133 
00142     typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
00143 
00149     typedef double magnitude_type;
00150 
00151   private:
00152     typedef TSQR::MatView<ordinal_type, scalar_type> matview_type;
00153     typedef TSQR::NodeTsqrFactory<node_type, scalar_type, ordinal_type> node_tsqr_factory_type;
00154     // Don't need a "typename" here, because there are no template
00155     // parameters involved in the type definition.
00156     typedef node_tsqr_factory_type::node_tsqr_type node_tsqr_type;
00157     typedef TSQR::DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
00158     typedef TSQR::Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;
00159 
00160   public:
00167     TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
00168       nodeTsqr_ (new node_tsqr_type),
00169       distTsqr_ (new dist_tsqr_type),
00170       tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
00171       ready_ (false)
00172     {
00173       setParameterList (plist);
00174     }
00175 
00177     TsqrAdaptor () : 
00178       nodeTsqr_ (new node_tsqr_type),
00179       distTsqr_ (new dist_tsqr_type),
00180       tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
00181       ready_ (false)
00182     {
00183       setParameterList (Teuchos::null);
00184     }
00185 
00186     Teuchos::RCP<const Teuchos::ParameterList>
00187     getValidParameters () const
00188     {
00189       using Teuchos::RCP;
00190       using Teuchos::rcp;
00191       using Teuchos::ParameterList;
00192       using Teuchos::parameterList;
00193 
00194       if (defaultParams_.is_null()) {
00195   RCP<ParameterList> params = parameterList ("TSQR implementation");
00196   params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
00197   params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
00198   defaultParams_ = params;
00199       }
00200       return defaultParams_;
00201     }
00202 
00203     void 
00204     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00205     {
00206       using Teuchos::ParameterList;
00207       using Teuchos::parameterList;
00208       using Teuchos::RCP;
00209       using Teuchos::sublist;
00210 
00211       RCP<ParameterList> params = plist.is_null() ? 
00212   parameterList (*getValidParameters ()) : plist;
00213       nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
00214       distTsqr_->setParameterList (sublist (params, "DistTsqr"));
00215 
00216       this->setMyParamList (params);
00217     }
00218 
00240     void
00241     factorExplicit (MV& A,
00242         MV& Q,
00243         dense_matrix_type& R,
00244         const bool forceNonnegativeDiagonal=false)
00245     {
00246       typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00247 
00248       prepareTsqr (Q); // Finish initializing TSQR.
00249       KMV A_view = getNonConstView (A);
00250       KMV Q_view = getNonConstView (Q);
00251       tsqr_->factorExplicit (A_view, Q_view, R, false,
00252            forceNonnegativeDiagonal);
00253     }
00254 
00285     int
00286     revealRank (MV& Q,
00287     dense_matrix_type& R,
00288     const magnitude_type& tol)
00289     {
00290       typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00291 
00292       prepareTsqr (Q); // Finish initializing TSQR.      
00293 
00294       // FIXME (mfh 25 Oct 2010) Check Epetra_Comm object in Q to make
00295       // sure it is the same communicator as the one we are using in
00296       // our dist_tsqr_type implementation.
00297       KMV Q_view = getNonConstView (Q);
00298       return tsqr_->revealRank (Q_view, R, tol, false);
00299     }
00300 
00301   private:
00303     Teuchos::RCP<node_tsqr_type> nodeTsqr_;
00304     
00306     Teuchos::RCP<dist_tsqr_type> distTsqr_;
00307 
00309     Teuchos::RCP<tsqr_type> tsqr_;
00310 
00312     mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00313 
00315     bool ready_;
00316 
00335     void 
00336     prepareTsqr (const MV& mv) 
00337     {
00338       if (! ready_) {
00339   prepareDistTsqr (mv);
00340   prepareNodeTsqr (mv);
00341   ready_ = true;
00342       }
00343     }
00344 
00348     void
00349     prepareNodeTsqr (const MV& mv)
00350     {
00351       (void) mv; // Epetra objects don't have a Kokkos Node.
00352 
00353       // Kokkos::SerialNode wants an empty ParameterList.
00354       Teuchos::ParameterList plist;
00355       Teuchos::RCP<node_type> node (new node_type (plist));
00356       node_tsqr_factory_type::prepareNodeTsqr (nodeTsqr_, node);
00357     }
00358 
00365     void
00366     prepareDistTsqr (const MV& mv)
00367     {
00368       using Teuchos::RCP;
00369       using Teuchos::rcp;
00370       using TSQR::Epetra::makeTsqrMessenger;
00371       typedef TSQR::MessengerBase<scalar_type> base_mess_type;
00372 
00373       // If mv falls out of scope, its Epetra_Comm may become invalid.
00374       // Thus, we clone the input Epetra_Comm, so that the messenger
00375       // owns the object.
00376       RCP<const Epetra_Comm> comm = rcp (mv.Comm().Clone());
00377       RCP<base_mess_type> messBase = makeTsqrMessenger<scalar_type> (comm);
00378       distTsqr_->init (messBase);
00379     }
00380 
00393     static Kokkos::MultiVector<scalar_type, node_type>
00394     getNonConstView (MV& A)
00395     {
00396       // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
00397       // storage of A uses nonconstant stride internally.  We would
00398       // have to copy and pack into a matrix with constant stride, and
00399       // then unpack on exit.  For now we choose just to raise an
00400       // exception.
00401       TEUCHOS_TEST_FOR_EXCEPTION(! A.ConstantStride(), std::invalid_argument,
00402          "TSQR does not currently support Epetra_MultiVector "
00403          "inputs that do not have constant stride.");
00404       const int numRows = A.MyLength();
00405       const int numCols = A.NumVectors();
00406       const int stride  = A.Stride();
00407       // A_ptr does _not_ own the data.  TSQR only operates within the
00408       // scope of the multivector objects on which it operates, so it
00409       // doesn't need ownership of the data.
00410       Teuchos::ArrayRCP<double> A_ptr (A.Values(), 0, numRows*stride, false);
00411 
00412       typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00413       // KMV objects want a Kokkos Node instance.  Epetra objects
00414       // don't have a Kokkos Node, so we make a temporary node just
00415       // for the KMV.
00416       //
00417       // Kokkos::SerialNode wants an empty ParameterList.
00418       Teuchos::ParameterList plist;
00419       Teuchos::RCP<node_type> node (new node_type (plist));
00420       KMV A_kmv (node);
00421       A_kmv.initializeValues (numRows, numCols, A_ptr, stride);
00422       return A_kmv;
00423     }
00424   };
00425 
00426 } // namespace Epetra
00427 
00428 #endif // defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
00429 
00430 #endif // __Epetra_TsqrAdaptor_hpp
00431 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines