Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Superlu_decl.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //           Amesos2: Templated Direct Sparse Solver Package
00006 //                  Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //
00042 // @HEADER
00043 
00053 #ifndef AMESOS2_SUPERLU_DECL_HPP
00054 #define AMESOS2_SUPERLU_DECL_HPP
00055 
00056 #include "Amesos2_SolverTraits.hpp"
00057 #include "Amesos2_SolverCore.hpp"
00058 #include "Amesos2_Superlu_FunctionMap.hpp"
00059 
00060 
00061 namespace Amesos2 {
00062 
00063 
00071 template <class Matrix,
00072           class Vector>
00073 class Superlu : public SolverCore<Amesos2::Superlu, Matrix, Vector>
00074 {
00075   friend class SolverCore<Amesos2::Superlu,Matrix,Vector>; // Give our base access
00076                                                           // to our private
00077                                                           // implementation funcs
00078 public:
00079 
00081   static const char* name;      // declaration. Initialization outside.
00082 
00083   typedef Superlu<Matrix,Vector>                                       type;
00084   typedef SolverCore<Amesos2::Superlu,Matrix,Vector>             super_type;
00085 
00086   // Since typedef's are not inheritted, go grab them
00087   typedef typename super_type::scalar_type                      scalar_type;
00088   typedef typename super_type::local_ordinal_type        local_ordinal_type;
00089   typedef typename super_type::global_ordinal_type      global_ordinal_type;
00090   typedef typename super_type::global_size_type            global_size_type;
00091 
00092   typedef TypeMap<Amesos2::Superlu,scalar_type>                    type_map;
00093 
00094   /*
00095    * The SuperLU interface will need two other typedef's, which are:
00096    * - the superlu type that corresponds to scalar_type and
00097    * - the corresponding type to use for magnitude
00098    */
00099   typedef typename type_map::type                                  slu_type;
00100   typedef typename type_map::magnitude_type                  magnitude_type;
00101 
00102   typedef FunctionMap<Amesos2::Superlu,slu_type>               function_map;
00103 
00105 
00106 
00113   Superlu(Teuchos::RCP<const Matrix> A,
00114           Teuchos::RCP<Vector>       X,
00115           Teuchos::RCP<const Vector> B);
00116 
00117 
00119   ~Superlu( );
00120 
00122 
00123 private:
00124 
00130   int preOrdering_impl();
00131 
00132 
00140   int symbolicFactorization_impl();
00141 
00142 
00148   int numericFactorization_impl();
00149 
00150 
00162   int solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> >       X,
00163                  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const;
00164 
00165 
00169   bool matrixShapeOK_impl() const;
00170 
00171 
00199   void setParameters_impl(
00200     const Teuchos::RCP<Teuchos::ParameterList> & parameterList );
00201 
00202 
00209   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters_impl() const;
00210 
00211 
00220   bool loadA_impl(EPhase current_phase);
00221 
00222 
00223   // struct holds all data necessary to make a superlu factorization or solve call
00224   mutable struct SLUData {
00225     SLU::SuperMatrix A, B, X, L, U; // matrix A in NCformat
00226     SLU::SuperMatrix AC;        // permuted matrix A in NCPformat
00227 
00228     SLU::superlu_options_t options;
00229     SLU::mem_usage_t mem_usage;
00230     SLU::SuperLUStat_t stat;
00231 
00232     Teuchos::Array<magnitude_type> berr; 
00233     Teuchos::Array<magnitude_type> ferr; 
00234     Teuchos::Array<int> perm_r;
00235     Teuchos::Array<int> perm_c;
00236     Teuchos::Array<int> etree;
00237     Teuchos::Array<magnitude_type> R;
00238     Teuchos::Array<magnitude_type> C;
00239 
00240     char equed;
00241     bool rowequ, colequ;        // flags what type of equilibration
00242                                 // has been performed
00243 
00244     int relax;
00245     int panel_size;
00246   } data_;
00247 
00248   // The following Arrays are persisting storage arrays for A, X, and B
00250   Teuchos::Array<slu_type> nzvals_;
00252   Teuchos::Array<int> rowind_;
00254   Teuchos::Array<int> colptr_;
00255 
00257   Teuchos::Array<slu_type> xvals_;  int ldx_;
00259   Teuchos::Array<slu_type> bvals_;  int ldb_;
00260 
00261   /* Note: In the above, must use "Amesos2::Superlu" rather than
00262    * "Superlu" because otherwise the compiler references the
00263    * specialized type of the class, and not the templated type that is
00264    * required for Amesos2::TypeMap
00265    */
00266 
00267   /* SuperLU can accept input in either compressed-row or
00268    * compressed-column storage.  We will store and pass matrices in
00269    * *compressed-column* format.
00270    */
00271 
00272   /*
00273    * Internal flag that is used for the numericFactorization_impl
00274    * routine.  If true, then the superlu gstrf routine should have
00275    * SamePattern_SameRowPerm in its options.  Otherwise, it should
00276    * factor from scratch.
00277    *
00278    * This is somewhat of a kludge to get around the fact that the
00279    * superlu routines all expect something different from the options
00280    * struct.  The big issue is that we don't want gstrf doing the
00281    * symbolic factorization if it doesn't need to.  On the other hand,
00282    * we can't leave options.Fact set to SamePattern_SameRowPerm
00283    * because the solver driver needs it to be set at FACTORED.  But
00284    * having it set at FACTORED upon re-entrance into
00285    * numericFactorization prompts gstrf to redo the symbolic
00286    * factorization.
00287    */
00288   bool same_symbolic_;
00289 
00290 };                              // End class Superlu
00291 
00292 
00293 // Specialize solver_traits struct for SuperLU
00294 template <>
00295 struct solver_traits<Superlu> {
00296 #ifdef HAVE_TEUCHOS_COMPLEX
00297   typedef Meta::make_list6<float,
00298                            double,
00299                            std::complex<float>,
00300                            std::complex<double>,
00301                            SLU::C::complex,
00302                            SLU::Z::doublecomplex> supported_scalars;
00303 #else
00304   typedef Meta::make_list2<float, double> supported_scalars;
00305 #endif
00306 };
00307 
00308 } // end namespace Amesos2
00309 
00310 #endif  // AMESOS2_SUPERLU_DECL_HPP