Ifpack2 Templated Preconditioning Package: Ifpack2_Details_TriDiSolver

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_TRIDISOLVER_DECL_HPP
00044 #define IFPACK2_DETAILS_TRIDISOLVER_DECL_HPP
00045 
00048 
00049 #include "Ifpack2_ConfigDefs.hpp"
00050 #include "Ifpack2_Preconditioner.hpp"
00051 #include "Ifpack2_Details_CanChangeMatrix.hpp"
00052 #include "Teuchos_SerialTriDiMatrix.hpp"
00053 
00054 
00055 namespace Ifpack2 {
00056 namespace Details {
00057 
00058 #ifndef IFPACK2_LAPACKSUPPORTSSCALAR
00059 #define IFPACK2_LAPACKSUPPORTSSCALAR
00060 
00068 
00069 template<class ScalarType>
00070 struct LapackSupportsScalar {
00072   static const bool value = false;
00073 };
00074 
00075 template<>
00076 struct LapackSupportsScalar<float> {
00077   static const bool value = true;
00078 };
00079 
00080 template<>
00081 struct LapackSupportsScalar<double> {
00082   static const bool value = true;
00083 };
00084 
00085 // FIXME (mfh 15 Nov 2013) Should we write IFPACK2_HAVE_COMPLEX ?
00086 #ifdef TEUCHOS_HAVE_COMPLEX
00087 template<>
00088 struct LapackSupportsScalar<std::complex<float> > {
00089   static const bool value = true;
00090 };
00091 
00092 template<>
00093 struct LapackSupportsScalar<std::complex<double> > {
00094   static const bool value = true;
00095 };
00096 #endif // TEUCHOS_HAVE_COMPLEX
00097 #endif // IFPACK2_LAPACKSUPPORTSSCALAR
00098 
00111 template<class MatrixType,
00112          const bool stub = ! LapackSupportsScalar<typename MatrixType::scalar_type>::value>
00113 class TriDiSolver :
00114     public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
00115                                    typename MatrixType::local_ordinal_type,
00116                                    typename MatrixType::global_ordinal_type,
00117                                    typename MatrixType::node_type>,
00118     virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
00119                                                                        typename MatrixType::local_ordinal_type,
00120                                                                        typename MatrixType::global_ordinal_type,
00121                                                                        typename MatrixType::node_type> >
00122 {};
00123 
00125 template<class MatrixType>
00126 class TriDiSolver<MatrixType, false> :
00127     public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
00128                                    typename MatrixType::local_ordinal_type,
00129                                    typename MatrixType::global_ordinal_type,
00130                                    typename MatrixType::node_type>,
00131     virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
00132                                                                        typename MatrixType::local_ordinal_type,
00133                                                                        typename MatrixType::global_ordinal_type,
00134                                                                        typename MatrixType::node_type> >
00135 {
00136 public:
00138 
00139 
00144   typedef MatrixType matrix_type;
00145 
00147   typedef typename MatrixType::scalar_type scalar_type;
00148 
00150   typedef typename MatrixType::local_ordinal_type local_ordinal_type;
00151 
00153   typedef typename MatrixType::global_ordinal_type global_ordinal_type;
00154 
00156   typedef typename MatrixType::node_type node_type;
00157 
00159   typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
00160 
00162   typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> row_matrix_type;
00163 
00165   typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00166 
00168 
00169 
00170 
00174   TriDiSolver (const Teuchos::RCP<const row_matrix_type>& matrix);
00175 
00177   virtual ~TriDiSolver ();
00178 
00180 
00181 
00182 
00183 
00188   Teuchos::RCP<const map_type> getDomainMap () const;
00189 
00194   Teuchos::RCP<const map_type> getRangeMap () const;
00195 
00201   void
00202   apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00203          Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00204          Teuchos::ETransp mode = Teuchos::NO_TRANS,
00205          scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
00206          scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero()) const;
00207 
00209 
00211   void setParameters (const Teuchos::ParameterList& params);
00212 
00222   void initialize ();
00223 
00225   bool isInitialized () const;
00226 
00236   void compute ();
00237 
00239   bool isComputed () const;
00240 
00250   virtual magnitude_type TEUCHOS_DEPRECATED
00251   computeCondEst (CondestType CT = Ifpack2::Cheap,
00252                   local_ordinal_type MaxIters = 1550,
00253                   magnitude_type Tol = 1e-9,
00254                   const Teuchos::Ptr<const row_matrix_type>& Matrix = Teuchos::null);
00255 
00259   virtual magnitude_type TEUCHOS_DEPRECATED getCondEst() const;
00260 
00262   Teuchos::RCP<const row_matrix_type> getMatrix () const;
00263 
00265   void setMatrix (const Teuchos::RCP<const row_matrix_type>& A);
00266 
00268   int getNumInitialize () const;
00269 
00271   int getNumCompute () const;
00272 
00274   int getNumApply() const;
00275 
00277   double getInitializeTime() const;
00278 
00280   double getComputeTime() const;
00281 
00283   double getApplyTime() const;
00284 
00286 
00287 
00288 
00290   std::string description () const;
00291 
00293   void
00294   describe (Teuchos::FancyOStream &out,
00295             const Teuchos::EVerbosityLevel verbLevel =
00296             Teuchos::Describable::verbLevel_default) const;
00298 private:
00299 
00301   void
00302   describeLocal (Teuchos::FancyOStream& out,
00303                  const Teuchos::EVerbosityLevel verbLevel) const;
00304 
00306   void reset ();
00307 
00315   static void
00316   extract (Teuchos::SerialTriDiMatrix<int, scalar_type>& A_local_tridi,
00317            const row_matrix_type& A_local);
00318 
00328   static void
00329   factor (Teuchos::SerialTriDiMatrix<int, scalar_type>& A,
00330           const Teuchos::ArrayView<int>& ipiv);
00331 
00333   typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
00334                               global_ordinal_type, node_type> MV;
00335 
00337   typedef Tpetra::Import<local_ordinal_type,
00338                          global_ordinal_type, node_type> import_type;
00339 
00341   typedef Tpetra::Export<local_ordinal_type,
00342                          global_ordinal_type, node_type> export_type;
00343 
00345   typedef Teuchos::ScalarTraits<scalar_type> STS;
00346 
00355   void
00356   applyImpl (const MV& X,
00357              MV& Y,
00358              const Teuchos::ETransp mode,
00359              const scalar_type alpha,
00360              const scalar_type beta) const;
00361 
00363   Teuchos::RCP<const row_matrix_type> A_;
00364 
00366   Teuchos::RCP<const row_matrix_type> A_local_;
00367 
00369   Teuchos::SerialTriDiMatrix<int, scalar_type> A_local_tridi_;
00370 
00372   Teuchos::Array<int> ipiv_;
00373 
00375   double initializeTime_;
00376 
00378   double computeTime_;
00379 
00381   mutable double applyTime_;
00382 
00384   int numInitialize_;
00385 
00387   int numCompute_;
00388 
00390   mutable int numApply_;
00391 
00393   bool isInitialized_;
00394 
00396   bool isComputed_;
00397 };
00398 
00399 
00401 template<class MatrixType>
00402 class TriDiSolver<MatrixType, true> :
00403     public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
00404                                    typename MatrixType::local_ordinal_type,
00405                                    typename MatrixType::global_ordinal_type,
00406                                    typename MatrixType::node_type>,
00407     virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
00408                                                                        typename MatrixType::local_ordinal_type,
00409                                                                        typename MatrixType::global_ordinal_type,
00410                                                                        typename MatrixType::node_type> >
00411 {
00412 public:
00414 
00415 
00420   typedef MatrixType matrix_type;
00421 
00423   typedef typename MatrixType::scalar_type scalar_type;
00424 
00426   typedef typename MatrixType::local_ordinal_type local_ordinal_type;
00427 
00429   typedef typename MatrixType::global_ordinal_type global_ordinal_type;
00430 
00432   typedef typename MatrixType::node_type node_type;
00433 
00435   typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
00436 
00438   typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> row_matrix_type;
00439 
00441   typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00442 
00444 
00445 
00446 
00450   TriDiSolver (const Teuchos::RCP<const row_matrix_type>& matrix);
00451 
00453   virtual ~TriDiSolver ();
00454 
00456 
00457 
00458 
00463   Teuchos::RCP<const map_type> getDomainMap () const;
00464 
00469   Teuchos::RCP<const map_type> getRangeMap () const;
00470 
00476   void
00477   apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00478          Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00479          Teuchos::ETransp mode = Teuchos::NO_TRANS,
00480          scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
00481          scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero()) const;
00482 
00484 
00486   void setParameters (const Teuchos::ParameterList& params);
00487 
00497   void initialize ();
00498 
00500   bool isInitialized () const;
00501 
00511   void compute ();
00512 
00514   bool isComputed () const;
00515 
00517   magnitude_type TEUCHOS_DEPRECATED
00518   computeCondEst (CondestType CT = Ifpack2::Cheap,
00519                   local_ordinal_type MaxIters = 1550,
00520                   magnitude_type Tol = 1e-9,
00521                   const Teuchos::Ptr<const row_matrix_type>& Matrix = Teuchos::null);
00522 
00524   magnitude_type TEUCHOS_DEPRECATED getCondEst () const;
00525 
00527   Teuchos::RCP<const row_matrix_type> getMatrix () const;
00528 
00530   void setMatrix (const Teuchos::RCP<const row_matrix_type>& A);
00531 
00533   int getNumInitialize () const;
00534 
00536   int getNumCompute () const;
00537 
00539   int getNumApply() const;
00540 
00542   double getInitializeTime() const;
00543 
00545   double getComputeTime() const;
00546 
00548   double getApplyTime() const;
00549 
00551 
00552 
00553 
00555   std::string description () const;
00556 
00558   void
00559   describe (Teuchos::FancyOStream &out,
00560             const Teuchos::EVerbosityLevel verbLevel =
00561             Teuchos::Describable::verbLevel_default) const;
00562   
00563   void 
00564   describeLocal (Teuchos::FancyOStream& out,
00565           const Teuchos::EVerbosityLevel verbLevel) const;
00566 
00568 private:
00570   typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
00571                               global_ordinal_type, node_type> MV;
00572 };
00573 
00574 } // namespace Details
00575 } // namespace Ifpack2
00576 
00577 #endif // IFPACK2_DETAILS_TRIDISOLVER_DECL_HPP