Amesos2_Factory.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//           Amesos2: Templated Direct Sparse Solver Package
//                  Copyright 2011 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

#ifndef AMESOS2_FACTORY_HPP
#define AMESOS2_FACTORY_HPP

#include "Amesos2_config.h"

#include "Amesos2_Solver.hpp"
#include "Amesos2_SolverTraits.hpp"

#include "Teuchos_ScalarTraits.hpp"
#include "Amesos2_MultiVecAdapter.hpp"
#include "Amesos2_MatrixTraits.hpp"
#include "Amesos2_ctassert.hpp"

#if defined(HAVE_AMESOS2_KLU2) & !(defined HAVE_TPETRA_INST_COMPLEX_FLOAT)
#include "Amesos2_KLU2.hpp"
#endif
#ifdef HAVE_AMESOS2_SUPERLUDIST // Distributed-memory SuperLU
#include "Amesos2_Superludist.hpp"
#endif
#ifdef HAVE_AMESOS2_SUPERLUMT   // Multi-threaded SuperLU
#include "Amesos2_Superlumt.hpp"
#endif
#ifdef HAVE_AMESOS2_SUPERLU     // Sequential SuperLU
#include "Amesos2_Superlu.hpp"
#endif
#ifdef HAVE_AMESOS2_PARDISO_MKL // MKL version of Pardiso
#include "Amesos2_PardisoMKL.hpp"
#endif
#ifdef HAVE_AMESOS2_LAPACK
#include "Amesos2_Lapack.hpp"
#endif
#if defined (HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
#include "Amesos2_Cholmod.hpp"
#endif

namespace Amesos2 {

  template <class,class> class Solver;

  /*
   * Utility function to transform a string into all lowercase
   */
  std::string tolower(const std::string& s);


  template < class Matrix,
             class Vector >
  Solver<Matrix,Vector>*
  create(const Matrix* A, Vector* X, const Vector* B);


  template < class Matrix,
             class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(Teuchos::RCP<const Matrix> A,
         Teuchos::RCP<Vector>       X,
         Teuchos::RCP<const Vector> B);


  template < class Matrix,
             class Vector >
  Solver<Matrix,Vector>*
  create(const char* solverName, const Matrix* A, Vector* X, const Vector* B);


  template < class Matrix,
             class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(const char* solverName,
         const Teuchos::RCP<const Matrix> A,
         const Teuchos::RCP<Vector>       X,
         const Teuchos::RCP<const Vector> B);


  template < class Matrix,
             class Vector >
  Solver<Matrix,Vector>*
  create(const std::string solverName, const Matrix* A, Vector* X, const Vector* B);


  template < class Matrix,
             class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(const std::string solverName,
         const Teuchos::RCP<const Matrix> A,
         const Teuchos::RCP<Vector>       X,
         const Teuchos::RCP<const Vector> B);


  template < class Matrix,
             class Vector >
  Solver<Matrix,Vector>*
  create(const std::string solverName, const Matrix* A);


  template < class Matrix,
             class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(const std::string solverName,
         const Teuchos::RCP<const Matrix> A);


  // Meta-functions to help with creation of solvers //

  template < template <class,class> class ConcreteSolver,
             class Matrix,
             class Vector >
  struct create_solver_with_supported_type {
    static Teuchos::RCP<Solver<Matrix,Vector> > apply(Teuchos::RCP<const Matrix> A,
                                                      Teuchos::RCP<Vector>       X,
                                                      Teuchos::RCP<const Vector> B )
    {
      ctassert<
        Meta::is_same<
          typename MatrixTraits<Matrix>::scalar_t,
          typename MultiVecAdapter<Vector>::scalar_t
        >::value
      > same_scalar_assertion;
      (void)same_scalar_assertion; // This stops the compiler from warning about unused declared variables

      // If our assertion did not fail, then create and return a new solver
      return rcp( new ConcreteSolver<Matrix,Vector>(A, X, B) );
    }
  };

template < template <class,class> class ConcreteSolver,
           class Matrix,
           class Vector >
struct throw_no_scalar_support_exception {
  static Teuchos::RCP<Solver<Matrix,Vector> > apply(Teuchos::RCP<const Matrix> A,
                                                    Teuchos::RCP<Vector>       X,
                                                    Teuchos::RCP<const Vector> B )
  {
    // typedef ConcreteSolver<Matrix,Vector> concretesolver_matrix_vector;
    typedef typename MatrixTraits<Matrix>::scalar_t scalar_t;
    TEUCHOS_TEST_FOR_EXCEPTION( true,
                        std::invalid_argument,
                        "The requested Amesos2 "
                        // << concretesolver_matrix_vector::name <<
                        " solver interface does not support the " <<
                        Teuchos::ScalarTraits<scalar_t>::name() <<
                        " scalar type." );
  }
};

  template < template <class,class> class ConcreteSolver,
             class Matrix,
             class Vector >
  struct handle_solver_type_support {
    static Teuchos::RCP<Solver<Matrix,Vector> > apply(Teuchos::RCP<const Matrix> A,
                                                      Teuchos::RCP<Vector>       X,
                                                      Teuchos::RCP<const Vector> B )
    {
      return Meta::if_then_else<
      solver_supports_scalar<ConcreteSolver, typename MatrixTraits<Matrix>::scalar_t>::value,
        create_solver_with_supported_type<ConcreteSolver,Matrix,Vector>,
        throw_no_scalar_support_exception<ConcreteSolver,Matrix,Vector> >::type::apply(A, X, B);
  }
    };


  // Query Functions //

  bool query(const char* solverName);


  bool query(const std::string solverName);


  // Definitions //

  template <class Matrix,
            class Vector >
  Solver<Matrix,Vector>*
  create(Matrix* A, Vector* X, Vector* B)
  {
    std::string solver = "Klu2";
    // Pass non-owning RCP objects to other factory method
    return( create(solver, rcp(A,false), rcp(X,false), rcp(B,false)).getRawPtr() );
  }


  template <class Matrix,
            class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(Teuchos::RCP<const Matrix> A,
         Teuchos::RCP<Vector>       X,
         Teuchos::RCP<const Vector> B)
  {
    std::string solver = "Klu2";
    return( create(solver, A, X, B) );
  }


  template <class Matrix,
            class Vector >
  Solver<Matrix,Vector>*
  create(const char* solverName, const Matrix* A, Vector* X, const Vector* B)
  {
    std::string solver = solverName;
    // Pass non-owning Teuchos::RCP objects to other factory method
    return( create(solver, rcp(A,false), rcp(X,false), rcp(B,false)).getRawPtr() );
  }


  template <class Matrix,
            class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(const char* solverName,
         const Teuchos::RCP<const Matrix> A,
         const Teuchos::RCP<Vector>       X,
         const Teuchos::RCP<const Vector> B)
  {
    std::string solver = solverName;
    return( create(solver, A, X, B) );
  }


  template <class Matrix,
            class Vector >
  Solver<Matrix,Vector>*
  create(const std::string solverName, const Matrix* A){
    return( create(solverName, rcp(A,false),
                   Teuchos::RCP<Vector>(),
                   Teuchos::RCP<const Vector>()).getRawPtr() );
  }


  template <class Matrix,
            class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(const std::string solverName, const Teuchos::RCP<const Matrix> A){
    return( create(solverName, A, Teuchos::RCP<Vector>(), Teuchos::RCP<const Vector>()) );
  }


  template <class Matrix,
            class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(const std::string solverName, const Matrix* A, Vector* X, const Vector* B)
  {
    // Pass non-owning Teuchos::RCP objects to other factory method
    return( create(solverName, rcp(A,false), rcp(X,false), rcp(B,false)) );
  }


  template <class Matrix,
            class Vector >
  Teuchos::RCP<Solver<Matrix,Vector> >
  create(const std::string solver_name,
         const Teuchos::RCP<const Matrix> A,
         const Teuchos::RCP<Vector>       X,
         const Teuchos::RCP<const Vector> B)
  {
    std::string solverName = tolower(solver_name); // for easy string checking

    // Check for our native solver first.  Treat KLU and KLU2 as equals.
    //
    // We use compiler guards in case a user does want to disable KLU2
#ifdef HAVE_AMESOS2_KLU2 

if((solverName == "amesos2_klu2") || (solverName == "klu2") ||
   (solverName == "amesos2_klu")  || (solverName == "klu")){
#ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
    std::string err_msg = solver_name + " is not is not supported with complex<float>";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, err_msg);
    return( Teuchos::null );
#else
      return handle_solver_type_support<KLU2,Matrix,Vector>::apply(A, X, B);
    }
#endif
#endif

#ifdef HAVE_AMESOS2_SUPERLUDIST
    if((solverName == "amesos2_superludist") ||
       (solverName == "superludist") ||
       (solverName == "amesos2_superlu_dist") ||
       (solverName == "superlu_dist")){
      return handle_solver_type_support<Superludist,Matrix,Vector>::apply(A, X, B);
    }
#endif

#ifdef HAVE_AMESOS2_SUPERLUMT
    if((solverName == "amesos2_superlumt") ||
       (solverName == "superlumt") ||
       (solverName == "amesos2_superlu_mt") ||
       (solverName == "superlu_mt")){
      return handle_solver_type_support<Superlumt,Matrix,Vector>::apply(A, X, B);
    }
#endif

#ifdef HAVE_AMESOS2_SUPERLU
    if((solverName == "amesos2_superlu") ||
       (solverName == "superlu")){
      return handle_solver_type_support<Superlu,Matrix,Vector>::apply(A, X, B);
    }
#endif

#ifdef HAVE_AMESOS2_PARDISO_MKL
    if((solverName == "amesos2_pardiso_mkl") ||
       (solverName == "pardiso_mkl") ||
       (solverName == "amesos2_pardisomkl")  ||
       (solverName == "pardisomkl")){
      return handle_solver_type_support<PardisoMKL,Matrix,Vector>::apply(A, X, B);
    }
#endif

#ifdef HAVE_AMESOS2_LAPACK
    if((solverName == "amesos2_lapack") ||
       (solverName == "lapack")){
      return handle_solver_type_support<Lapack,Matrix,Vector>::apply(A, X, B);
    }
#endif

#if defined (HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
    if(solverName == "amesos2_cholmod")
      return handle_solver_type_support<Cholmod,Matrix,Vector>::apply(A, X, B);
#endif

    /* If none of the above conditionals are satisfied, then the solver
     * requested is not yet supported.  We throw a runtime exception stating
     * this, and return null.
     */
    err_msg = solver_name + " is not enabled or is not supported";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, err_msg);
    return( Teuchos::null );
  }

} // end namespace Amesos2

#endif  // AMESOS2_FACTORY_HPP