Ifpack2 Templated Preconditioning Package: Ifpack2_Details_DenseSolver

00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
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 
00043 #ifndef IFPACK2_DETAILS_DENSESOLVER_DECL_HPP
00044 #define IFPACK2_DETAILS_DENSESOLVER_DECL_HPP
00045 
00048 
00049 #include "Ifpack2_ConfigDefs.hpp"
00050 #include "Ifpack2_Preconditioner.hpp"
00051 #include "Ifpack2_Details_CanChangeMatrix.hpp"
00052 #include "Teuchos_SerialDenseMatrix.hpp"
00053 
00054 
00055 namespace Ifpack2 {
00056 namespace Details {
00057 
00058 #ifndef IFPACK2_LAPACKSUPPORTSSCALAR
00059 #define IFPACK2_LAPACKSUPPORTSSCALAR
00060 
00068 template<class ScalarType>
00069 struct LapackSupportsScalar {
00071   static const bool value = false;
00072 };
00073 
00074 template<>
00075 struct LapackSupportsScalar<float> {
00076   static const bool value = true;
00077 };
00078 
00079 template<>
00080 struct LapackSupportsScalar<double> {
00081   static const bool value = true;
00082 };
00083 
00084 // FIXME (mfh 15 Nov 2013) Should we write IFPACK2_HAVE_COMPLEX ?
00085 #ifdef TEUCHOS_HAVE_COMPLEX
00086 template<>
00087 struct LapackSupportsScalar<std::complex<float> > {
00088   static const bool value = true;
00089 };
00090 
00091 template<>
00092 struct LapackSupportsScalar<std::complex<double> > {
00093   static const bool value = true;
00094 };
00095 #endif // TEUCHOS_HAVE_COMPLEX
00096 #endif //  IFPACK2_LAPACKSUPPORTSSCALAR
00097 
00110 template<class MatrixType,
00111          const bool stub = ! LapackSupportsScalar<typename MatrixType::scalar_type>::value>
00112 class DenseSolver :
00113     public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
00114                                    typename MatrixType::local_ordinal_type,
00115                                    typename MatrixType::global_ordinal_type,
00116                                    typename MatrixType::node_type>,
00117     virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
00118                                                                        typename MatrixType::local_ordinal_type,
00119                                                                        typename MatrixType::global_ordinal_type,
00120                                                                        typename MatrixType::node_type> >
00121 {};
00122 
00124 template<class MatrixType>
00125 class DenseSolver<MatrixType, false> :
00126     public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
00127                                    typename MatrixType::local_ordinal_type,
00128                                    typename MatrixType::global_ordinal_type,
00129                                    typename MatrixType::node_type>,
00130     virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
00131                                                                        typename MatrixType::local_ordinal_type,
00132                                                                        typename MatrixType::global_ordinal_type,
00133                                                                        typename MatrixType::node_type> >
00134 {
00135 public:
00137 
00138 
00143   typedef MatrixType matrix_type;
00144 
00146   typedef typename MatrixType::scalar_type scalar_type;
00147 
00149   typedef typename MatrixType::local_ordinal_type local_ordinal_type;
00150 
00152   typedef typename MatrixType::global_ordinal_type global_ordinal_type;
00153 
00155   typedef typename MatrixType::node_type node_type;
00156 
00158   typedef typename ::Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
00159 
00161   typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> row_matrix_type;
00162 
00164   typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00165 
00167 
00168 
00169 
00173   DenseSolver (const ::Teuchos::RCP<const row_matrix_type>& matrix);
00174 
00176   virtual ~DenseSolver ();
00177 
00179 
00180 
00181 
00182 
00187   ::Teuchos::RCP<const map_type> getDomainMap () const;
00188 
00193   ::Teuchos::RCP<const map_type> getRangeMap () const;
00194 
00200   void
00201   apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00202          Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00203          ::Teuchos::ETransp mode = ::Teuchos::NO_TRANS,
00204          scalar_type alpha = ::Teuchos::ScalarTraits<scalar_type>::one(),
00205          scalar_type beta = ::Teuchos::ScalarTraits<scalar_type>::zero()) const;
00206 
00208 
00210   void setParameters (const ::Teuchos::ParameterList& params);
00211 
00221   void initialize ();
00222 
00224   bool isInitialized () const;
00225 
00235   void compute ();
00236 
00238   bool isComputed () const;
00239 
00249   virtual magnitude_type TEUCHOS_DEPRECATED
00250   computeCondEst (CondestType CT = Ifpack2::Cheap,
00251                   local_ordinal_type MaxIters = 1550,
00252                   magnitude_type Tol = 1e-9,
00253                   const ::Teuchos::Ptr<const row_matrix_type>& Matrix = ::Teuchos::null);
00254 
00258   virtual magnitude_type TEUCHOS_DEPRECATED getCondEst() const;
00259 
00261   ::Teuchos::RCP<const row_matrix_type> getMatrix () const;
00262 
00264   void setMatrix (const ::Teuchos::RCP<const row_matrix_type>& A);
00265 
00267   int getNumInitialize () const;
00268 
00270   int getNumCompute () const;
00271 
00273   int getNumApply() const;
00274 
00276   double getInitializeTime() const;
00277 
00279   double getComputeTime() const;
00280 
00282   double getApplyTime() const;
00283 
00285 
00286 
00287 
00289   std::string description () const;
00290 
00292   void
00293   describe (::Teuchos::FancyOStream &out,
00294             const ::Teuchos::EVerbosityLevel verbLevel =
00295             ::Teuchos::Describable::verbLevel_default) const;
00297 private:
00298 
00300   void
00301   describeLocal (::Teuchos::FancyOStream& out,
00302                  const ::Teuchos::EVerbosityLevel verbLevel) const;
00303 
00305   void reset ();
00306 
00314   static void
00315   extract (::Teuchos::SerialDenseMatrix<int, scalar_type>& A_local_dense,
00316            const row_matrix_type& A_local);
00317 
00327   static void
00328   factor (::Teuchos::SerialDenseMatrix<int, scalar_type>& A,
00329           const ::Teuchos::ArrayView<int>& ipiv);
00330 
00332   typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
00333                               global_ordinal_type, node_type> MV;
00334 
00336   typedef Tpetra::Import<local_ordinal_type,
00337                          global_ordinal_type, node_type> import_type;
00338 
00340   typedef Tpetra::Export<local_ordinal_type,
00341                          global_ordinal_type, node_type> export_type;
00342 
00344   typedef ::Teuchos::ScalarTraits<scalar_type> STS;
00345 
00354   void
00355   applyImpl (const MV& X,
00356              MV& Y,
00357              const ::Teuchos::ETransp mode,
00358              const scalar_type alpha,
00359              const scalar_type beta) const;
00360 
00362   ::Teuchos::RCP<const row_matrix_type> A_;
00363 
00365   ::Teuchos::RCP<const row_matrix_type> A_local_;
00366 
00368   ::Teuchos::SerialDenseMatrix<int, scalar_type> A_local_dense_;
00369 
00371   ::Teuchos::Array<int> ipiv_;
00372 
00374   double initializeTime_;
00375 
00377   double computeTime_;
00378 
00380   mutable double applyTime_;
00381 
00383   int numInitialize_;
00384 
00386   int numCompute_;
00387 
00389   mutable int numApply_;
00390 
00392   bool isInitialized_;
00393 
00395   bool isComputed_;
00396 };
00397 
00398 
00400 template<class MatrixType>
00401 class DenseSolver<MatrixType, true> :
00402     public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
00403                                    typename MatrixType::local_ordinal_type,
00404                                    typename MatrixType::global_ordinal_type,
00405                                    typename MatrixType::node_type>,
00406     virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
00407                                                                        typename MatrixType::local_ordinal_type,
00408                                                                        typename MatrixType::global_ordinal_type,
00409                                                                        typename MatrixType::node_type> >
00410 {
00411 public:
00413 
00414 
00419   typedef MatrixType matrix_type;
00420 
00422   typedef typename MatrixType::scalar_type scalar_type;
00423 
00425   typedef typename MatrixType::local_ordinal_type local_ordinal_type;
00426 
00428   typedef typename MatrixType::global_ordinal_type global_ordinal_type;
00429 
00431   typedef typename MatrixType::node_type node_type;
00432 
00434   typedef typename ::Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
00435 
00437   typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> row_matrix_type;
00438 
00440   typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00441 
00443 
00444 
00445 
00449   DenseSolver (const ::Teuchos::RCP<const row_matrix_type>& matrix);
00450 
00452   virtual ~DenseSolver ();
00453 
00455 
00456 
00457 
00462   ::Teuchos::RCP<const map_type> getDomainMap () const;
00463 
00468   ::Teuchos::RCP<const map_type> getRangeMap () const;
00469 
00475   void
00476   apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00477          Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00478          ::Teuchos::ETransp mode = ::Teuchos::NO_TRANS,
00479          scalar_type alpha = ::Teuchos::ScalarTraits<scalar_type>::one(),
00480          scalar_type beta = ::Teuchos::ScalarTraits<scalar_type>::zero()) const;
00481 
00483 
00485   void setParameters (const ::Teuchos::ParameterList& params);
00486 
00496   void initialize ();
00497 
00499   bool isInitialized () const;
00500 
00510   void compute ();
00511 
00513   bool isComputed () const;
00514 
00516   magnitude_type TEUCHOS_DEPRECATED
00517   computeCondEst (CondestType CT = Ifpack2::Cheap,
00518                   local_ordinal_type MaxIters = 1550,
00519                   magnitude_type Tol = 1e-9,
00520                   const ::Teuchos::Ptr<const row_matrix_type>& Matrix = ::Teuchos::null);
00521 
00523   magnitude_type TEUCHOS_DEPRECATED getCondEst () const;
00524 
00526   ::Teuchos::RCP<const row_matrix_type> getMatrix () const;
00527 
00529   void setMatrix (const ::Teuchos::RCP<const row_matrix_type>& A);
00530 
00532   int getNumInitialize () const;
00533 
00535   int getNumCompute () const;
00536 
00538   int getNumApply() const;
00539 
00541   double getInitializeTime() const;
00542 
00544   double getComputeTime() const;
00545 
00547   double getApplyTime() const;
00548 
00550 
00551 
00552 
00554   std::string description () const;
00555 
00557   void
00558   describe (::Teuchos::FancyOStream &out,
00559             const ::Teuchos::EVerbosityLevel verbLevel =
00560             ::Teuchos::Describable::verbLevel_default) const;
00562 private:
00564   typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
00565                               global_ordinal_type, node_type> MV;
00566 };
00567 
00568 } // namespace Details
00569 } // namespace Ifpack2
00570 
00571 #endif // IFPACK2_DETAILS_DENSESOLVER_DECL_HPP