Amesos2_PardisoMKL_def.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//           Amesos2: Templated Direct Sparse Solver Package
//
// Copyright 2011 NTESS and the Amesos2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
 
#ifndef AMESOS2_PARDISOMKL_DEF_HPP
#define AMESOS2_PARDISOMKL_DEF_HPP
 
#include <Teuchos_Tuple.hpp>
#include <Teuchos_toString.hpp>
#include <Teuchos_StandardParameterEntryValidators.hpp>
 
#include "Amesos2_SolverCore_def.hpp"
#include "Amesos2_PardisoMKL_decl.hpp"
 
 
namespace Amesos2 {
 
  namespace PMKL {
#   include <mkl.h>
#   include <mkl_pardiso.h>
  }
 
  template <class Matrix, class Vector>
  PardisoMKL<Matrix,Vector>::PardisoMKL(Teuchos::RCP<const Matrix> A,
                                        Teuchos::RCP<Vector>       X,
                                        Teuchos::RCP<const Vector> B)
    : SolverCore<Amesos2::PardisoMKL,Matrix,Vector>(A, X, B) // instantiate superclass
    , n_(Teuchos::as<int_t>(this->globalNumRows_))
    , perm_(this->globalNumRows_)
    , nrhs_(0)
    , partial_facto_ (0)
    , schur_size_ (0)
    , schur_out_ptr_(nullptr)
    , pardiso_initialized_(false)
    , only_forward_solve_(false)
    , only_backward_solve_(false)
    , is_contiguous_(true)
    , msglvl_(0)
    , debug_level_(0)
  {
    // set the default matrix type
    set_pardiso_mkl_matrix_type();
 
    PMKL::_INTEGER_t iparm_temp[64];
    PMKL::_INTEGER_t mtype_temp = mtype_;
    PMKL::pardisoinit(pt_, &mtype_temp, iparm_temp);
 
    for( int i = 0; i < 64; ++i ){
      iparm_[i] = iparm_temp[i];
    }
    iparm_[0] = 1; // do not use solver defaults
 
    // set single or double precision
    if constexpr ( std::is_same_v<solver_magnitude_type, PMKL::_REAL_t> ) {
      iparm_[27] = 1;           // single-precision
    } else {
      iparm_[27] = 0;           // double-precision
    }
 
    // Reset some of the default parameters
    iparm_[34] = 1;             // Use zero-based indexing
#ifdef HAVE_AMESOS2_DEBUG
    iparm_[26] = 1;             // turn the Pardiso matrix checker on
#endif
  }
 
 
  template <class Matrix, class Vector>
  PardisoMKL<Matrix,Vector>::~PardisoMKL( )
  {
    /*
     * Free any memory allocated by the PardisoMKL library functions
     */
    int_t error = 0;
    void *bdummy, *xdummy;
    if( this->root_ && pardiso_initialized_){
      int_t phase = -1;         // release all internal solver memory
      function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
                             const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
                             nzvals_view_.data(), rowptr_view_.data(),
                             colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
                             const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &error );
      pardiso_initialized_ = false;
    }
 
    check_pardiso_mkl_error(Amesos2::CLEAN, error);
  }
 
 
  template<class Matrix, class Vector>
  int
  PardisoMKL<Matrix,Vector>::preOrdering_impl()
  {
    // preOrdering done in PardisoMKL during "Analysis" (aka symbolic
    // factorization) phase
 
    return(0);
  }
 
 
  template <class Matrix, class Vector>
  int
  PardisoMKL<Matrix,Vector>::symbolicFactorization_impl()
  {
    using Teuchos::as;
 
    int_t error = 0;
    if( this->root_ ){
#ifdef HAVE_AMESOS2_TIMERS
      Teuchos::TimeMonitor symbFactTimer( this->timers_.symFactTime_ );
#endif
 
      solver_scalar_type bdummy, xdummy; // to be cased to void*
      if( pardiso_initialized_){
        int_t phase = -1;         // release all internal solver memory
        function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
                               const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
                               nzvals_view_.data(), rowptr_view_.data(),
                               colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
                               const_cast<int_t*>(&msglvl_), as<void*>(&bdummy), as<void*>(&xdummy), &error );
        if (msglvl_ > 0 && error != 0) {
          if (error != 0) {
            std::cout << " PardisoMKL::symbolicFactorization: clean-up failed with " << error << std::endl;
          } else {
            std::cout << " PardisoMKL::symbolicFactorization: cleaned-up before calling symbolic" << error;
          }
        }
        pardiso_initialized_ = false;
      }
      int_t phase = 11;
      function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
                             const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
                             nzvals_view_.data(), rowptr_view_.data(),
                             colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
                             const_cast<int_t*>(&msglvl_), as<void*>(&bdummy), as<void*>(&xdummy), &error );
      pardiso_initialized_ = true;
    }
    check_pardiso_mkl_error(Amesos2::SYMBFACT, error);
 
    if (msglvl_ > 0 && this->root_) {
      std::cout << " PardisoMKL::symbolicFactorization done:" << std::endl;
#ifdef HAVE_AMESOS2_TIMERS
      std::cout << " * Time : " << this->timers_.symFactTime_.totalElapsedTime() << std::endl;
#else
      std::cout << " * Time : not enabled" << std::endl;
#endif
    }
 
    // Pardiso only lets you retrieve the total number of factor
    // non-zeros, not for each individually.  We should document how
    // such a situation is reported.
    this->setNnzLU(iparm_[17]);
 
    return(0);
  }
 
 
  template <class Matrix, class Vector>
  int
  PardisoMKL<Matrix,Vector>::numericFactorization_impl()
  {
    using Teuchos::as;
 
    int_t error = 0;
    if( this->root_ ){
#ifdef HAVE_AMESOS2_TIMERS
      Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
#endif
 
      int_t phase = 22;
      solver_scalar_type bdummy; // to be casted to void*
      solver_scalar_type *xdummy = schur_out_.data();
      function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
                             const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
                             nzvals_view_.data(), rowptr_view_.data(),
                             colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
                             const_cast<int_t*>(&msglvl_), as<void*>(&bdummy), as<void*>(xdummy), &error );
 
      if (error == 0 && partial_facto_ != 0) {
        if (schur_out_ptr_ != nullptr) {
          // copy schur out if the output pointer is valid (assuming enough space has been allocated)
          // Pardiso returns the Schur complement in row-major
          // So, we transpose it to store in column-major
          for (size_t i = 0; i < schur_size_; i++) {
            for (size_t j = 0; j < schur_size_; j++) {
              schur_out_ptr_[i+j*schur_size_] = as<scalar_type>(schur_out_[j+i*schur_size_]);
            }
          }
        }
      }
    }
    check_pardiso_mkl_error(Amesos2::NUMFACT, error);
 
    if (msglvl_ > 0 && this->root_) {
      std::cout << " PardisoMKL::numericFactorization done:" << std::endl;
#ifdef HAVE_AMESOS2_TIMERS
      std::cout << " * Time : " << this->timers_.numFactTime_.totalElapsedTime() << std::endl;
#else
      std::cout << " * Time : not enabled" << std::endl;
#endif
    }
 
    return( 0 );
  }
 
 
  template <class Matrix, class Vector>
  int
  PardisoMKL<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> >       X,
                                        const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
  {
    using Teuchos::as;
 
    int_t error = 0;
 
    // Get B data
    const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
    nrhs_ = as<int_t>(X->getGlobalNumVectors());
    if (debug_level_ > 0) {
      if (this->root_) std::cout << "\n == Amesos2_PardisoMKL::solve_impl ==" << std::endl;
      if (debug_level_ == 1) {
        B->description();
      } else {
        Teuchos::RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
        if (!is_null(B->getMap())) B->getMap()->describe(*fancy, Teuchos::VERB_EXTREME);
        std::cout << std::endl;
        B->describe(*fancy, Teuchos::VERB_EXTREME);
      }
    }
 
    {                             // Get values from RHS B
#ifdef HAVE_AMESOS2_TIMERS
      Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
      Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
#endif
 
      const bool initialize_data = true;
      const bool do_not_initialize_data = false;
      Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
        host_solver_scalar_view>::do_get(initialize_data, B, bvals_,
          as<size_t>(ld_rhs),
          (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
          this->rowIndexBase_);
      Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
        host_solver_scalar_view>::do_get(do_not_initialize_data, X, xvals_,
          as<size_t>(ld_rhs),
          (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
          this->rowIndexBase_);
    }
 
    if( this->root_ ){
#ifdef HAVE_AMESOS2_TIMERS
      Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
#endif
 
      int_t phase = 33; // forward & backward solve
      if (only_forward_solve_) {
        // forward solve
        phase = 331;
        if (wvals_.extent(0) != n_ || wvals_.extent(1) != nrhs_) {
          Kokkos::resize(wvals_, n_, nrhs_);
        }
      } else if (only_backward_solve_) {
        // backward solve
        phase = 333;
        if (wvals_.extent(0) != n_ || wvals_.extent(1) != nrhs_) {
          Kokkos::resize(wvals_, n_, nrhs_);
        }
        // for consistent interface with other solvers (e.g., ShyLU-Basker)
        //  scatter RHS vectors (based on schur-part)
        //  input RHS for backward-solve has interior part followed by schur complement
        size_t n1 = 0;
        size_t n2 = this->globalNumCols_-schur_size_;
        for (global_size_type i=0; i<n_; i++) {
          if (schur_part_(i) == 1) {
            for (size_t j=0; j<nrhs_; j++) {
              wvals_(i,j) = bvals_(n2,j);
            }
            n2 ++;
          } else {
            for (size_t j=0; j<nrhs_; j++) {
              wvals_(i,j) = bvals_(n1,j);
            }
            n1 ++;
          }
        }
      }
      solver_scalar_type *b_in  = (only_backward_solve_ ? wvals_.data() : bvals_.data());
      solver_scalar_type *x_out = (only_forward_solve_  ? wvals_.data() : xvals_.data());
      function_map::pardiso( pt_,
                             const_cast<int_t*>(&maxfct_),
                             const_cast<int_t*>(&mnum_),
                             const_cast<int_t*>(&mtype_),
                             const_cast<int_t*>(&phase),
                             const_cast<int_t*>(&n_),
                             const_cast<solver_scalar_type*>(nzvals_view_.data()),
                             const_cast<int_t*>(rowptr_view_.data()),
                             const_cast<int_t*>(colind_view_.data()),
                             const_cast<int_t*>(perm_.getRawPtr()),
                             &nrhs_,
                             const_cast<int_t*>(iparm_),
                             const_cast<int_t*>(&msglvl_),
                             as<void*>(b_in),
                             as<void*>(x_out), &error );
      if (only_forward_solve_) {
        // for consistent interface with other solvers (e.g., ShyLU-Basker)
        //  gather solution vectors (based on schur-part)
        //  output vector should have interior part followed by schur complement
        size_t n1 = 0;
        size_t n2 = this->globalNumCols_-schur_size_;
        for (global_size_type i=0; i<n_; i++) {
          if (schur_part_(i) == 1) {
            for (size_t j=0; j<nrhs_; j++) {
              xvals_(n2,j) = wvals_(i,j);
            }
            n2 ++;
          } else {
            for (size_t j=0; j<nrhs_; j++) {
              xvals_(n1,j) = wvals_(i,j);
            }
            n1 ++;
          }
        }
      }
    }
    check_pardiso_mkl_error(Amesos2::SOLVE, error);
 
    /* Export X from root to the global space */
    {
#ifdef HAVE_AMESOS2_TIMERS
      Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
#endif
 
      Util::put_1d_data_helper_kokkos_view<
        MultiVecAdapter<Vector>,host_solver_scalar_view>::do_put(X, xvals_,
          as<size_t>(ld_rhs),
          (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
          this->rowIndexBase_);
    }
    if (debug_level_ > 0) {
      if (debug_level_ == 1) {
        X->description();
      } else {
        Teuchos::RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
        if (!is_null(X->getMap())) X->getMap()->describe(*fancy, Teuchos::VERB_EXTREME);
        std::cout << std::endl;
        X->describe(*fancy, Teuchos::VERB_EXTREME);
      }
    }
    if (msglvl_ > 0 && this->root_) {
      std::cout << " PardisoMKL::solve done:" << std::endl;
#ifdef HAVE_AMESOS2_TIMERS
      std::cout << " * Time : " << this->timers_.vecRedistTime_.totalElapsedTime()
                << " + " << this->timers_.solveTime_.totalElapsedTime() << std::endl;
#else
      std::cout << " * Time : not enabled" << std::endl;
#endif
    }
 
    return( 0 );
}
 
 
  template <class Matrix, class Vector>
  bool
  PardisoMKL<Matrix,Vector>::matrixShapeOK_impl() const
  {
    // PardisoMKL supports square matrices
    return( this->globalNumRows_ == this->globalNumCols_ );
  }
 
 
  template <class Matrix, class Vector>
  void
  PardisoMKL<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
  {
    using Teuchos::RCP;
    using Teuchos::getIntegralValue;
    using Teuchos::ParameterEntryValidator;
 
    RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
 
    // Fill-in reordering: 0 = minimum degree, 2 = METIS 4.0.1 (default), 3 = METIS 5.1, 4 = AMD,
    if( parameterList->isParameter("IPARM(2)") )
    {
      RCP<const ParameterEntryValidator> fillin_validator = valid_params->getEntry("IPARM(2)").validator();
      parameterList->getEntry("IPARM(2)").setValidator(fillin_validator);
      iparm_[1] = getIntegralValue<int>(*parameterList, "IPARM(2)");
    }
 
    // Iterative-direct algorithm
    if( parameterList->isParameter("IPARM(4)") )
    {
      RCP<const ParameterEntryValidator> prec_validator = valid_params->getEntry("IPARM(4)").validator();
      parameterList->getEntry("IPARM(4)").setValidator(prec_validator);
      iparm_[3] = getIntegralValue<int>(*parameterList, "IPARM(4)");
    }
   
    // Max numbers of iterative refinement steps
    if( parameterList->isParameter("IPARM(8)") )
    {
      RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry("IPARM(8)").validator();
      parameterList->getEntry("IPARM(8)").setValidator(refine_validator);
      iparm_[7] = getIntegralValue<int>(*parameterList, "IPARM(8)");
    }
 
    // Perturb the pivot elements
    if( parameterList->isParameter("IPARM(10)") )
    {
      RCP<const ParameterEntryValidator> pivot_perturb_validator = valid_params->getEntry("IPARM(10)").validator();
      parameterList->getEntry("IPARM(10)").setValidator(pivot_perturb_validator);
      iparm_[9] = getIntegralValue<int>(*parameterList, "IPARM(10)");
    }
 
    // Scale vector for stability
    if( parameterList->isParameter("IPARM(11)") )
    {
      RCP<const ParameterEntryValidator> mwm_validator = valid_params->getEntry("IPARM(11)").validator();
      parameterList->getEntry("IPARM(11)").setValidator(mwm_validator);
      iparm_[10] = getIntegralValue<int>(*parameterList, "IPARM(11)");
    }
 
    // First check if the control object requests a transpose solve.
    // Then solver specific options can override this.
    iparm_[11] = this->control_.useTranspose_ ? 2 : 0;
    
    // Normal solve (0), or a transpose solve (1)
    if( parameterList->isParameter("IPARM(12)") )
    {
      RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("IPARM(12)").validator();
      parameterList->getEntry("IPARM(12)").setValidator(trans_validator);
      iparm_[11] = getIntegralValue<int>(*parameterList, "IPARM(12)");
    }
 
    // (Non-)symmetric matchings : detault 1 for nonsymmetric and 0 for symmetric matrix (default is nonsymmetric)
    if( parameterList->isParameter("IPARM(13)") )
    {
      RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("IPARM(13)").validator();
      parameterList->getEntry("IPARM(13)").setValidator(trans_validator);
      iparm_[12] = getIntegralValue<int>(*parameterList, "IPARM(13)");
    }
 
    // Output: Number of nonzeros in the factor LU
    if( parameterList->isParameter("IPARM(18)") )
    {
      RCP<const ParameterEntryValidator> report_validator = valid_params->getEntry("IPARM(18)").validator();
      parameterList->getEntry("IPARM(18)").setValidator(report_validator);
      iparm_[17] = getIntegralValue<int>(*parameterList, "IPARM(18)");
    }
   
    // Scheduling method for the parallel numerical factorization
    if( parameterList->isParameter("IPARM(24)") )
    {
      RCP<const ParameterEntryValidator> par_fact_validator = valid_params->getEntry("IPARM(24)").validator();
      parameterList->getEntry("IPARM(24)").setValidator(par_fact_validator);
      iparm_[23] = getIntegralValue<int>(*parameterList, "IPARM(24)");
    }
  
    // Parallelization scheme for the forward and backward solve
    if( parameterList->isParameter("IPARM(25)") )
    {
      RCP<const ParameterEntryValidator> par_fbsolve_validator = valid_params->getEntry("IPARM(25)").validator();
      parameterList->getEntry("IPARM(25)").setValidator(par_fbsolve_validator);
      iparm_[24] = getIntegralValue<int>(*parameterList, "IPARM(25)");
    } 
 
    // Check matrix
    if( parameterList->isParameter("IPARM(27)") )
    {
      RCP<const ParameterEntryValidator> check_validator = valid_params->getEntry("IPARM(27)").validator();
      parameterList->getEntry("IPARM(27)").setValidator(check_validator);
      iparm_[26] = getIntegralValue<int>(*parameterList, "IPARM(27)");
    }
 
    // Graph compression scheme for METIS.
    if( parameterList->isParameter("IPARM(60)") )
    {
      RCP<const ParameterEntryValidator> ooc_validator = valid_params->getEntry("IPARM(60)").validator();
      parameterList->getEntry("IPARM(60)").setValidator(ooc_validator);
      iparm_[59] = getIntegralValue<int>(*parameterList, "IPARM(60)");
    } 
    
    // Partial factorization
    // Note: partial_facto_ == 1 means iparam_[35] = 2 (compute and factor Schur)
    // while partial_facto_ == 2 means iparam_[35] = 1 (only compute Schur)
    if(parameterList->isParameter("PartialFacto")) {
      partial_facto_ = parameterList->get<int>("PartialFacto");
      if (partial_facto_ == 1)
        iparm_[35] = 2; // calculate and factor Schur complement
      else if (partial_facto_ == 2)
        iparm_[35] = 1; // calculate Schur complement
    }
    if(parameterList->isParameter("SchurPart")) {
      // copy schur-part to the internal view (in case the user-pointer go out of scope?)
      auto schur_part_ptr = parameterList->get<const local_ordinal_type*>("SchurPart");
      Kokkos::resize(schur_part_, this->globalNumCols_);
 
      schur_size_ = 0;
      for (global_size_type i=0; i<this->globalNumCols_; i++) {
        schur_part_(i) = schur_part_ptr[i]; // keep track or schur part
        perm_[i] = schur_part_ptr[i];       // input for pardiso
        if (perm_[i] == 1) schur_size_ ++;
      }
      // allocate internal storage to store the schur complement (user may not want it and may not provide a valid pointer?)
      size_t schur_size_2 = schur_size_*schur_size_;
      Kokkos::resize(schur_out_, (schur_size_2 > 0 ? schur_size_2 : 1));
 
    }
    if(parameterList->isParameter("SchurOut")) {
      // store schur-part to the internal view (if user wants the output, then the pointer should stay, so no need for internal view?)
      schur_out_ptr_ = parameterList->get<scalar_type*>("SchurOut");
    }
    if(parameterList->isParameter("OnlyForwardSolve")) {
      only_forward_solve_ = parameterList->get<bool>("OnlyForwardSolve");
    }
    if(parameterList->isParameter("OnlyBackwardSolve")) {
      only_backward_solve_ = parameterList->get<bool>("OnlyBackwardSolve");
    }
 
 
    if( parameterList->isParameter("IsContiguous") ){
      is_contiguous_ = parameterList->get<bool>("IsContiguous");
    }
    if( parameterList->isParameter("MessageLevel") ){
      msglvl_ = parameterList->get<int>("MessageLevel");
    }
    if(parameterList->isParameter("verbose")){
      bool verbose = parameterList->get<bool>("verbose");
      if (verbose) msglvl_ = 1;
    }
    if( parameterList->isParameter("DebugLevel") ){
      debug_level_ = parameterList->get<int>("DebugLevel");
    }
  }
 
 
/*
 * TODO: It would be nice if the parameters could be expressed as
 * either all string or as all integers.  I see no way of doing this
 * at present with the standard validators.  However, we could create
 * our own validators or kindly ask the Teuchos team to add some
 * features for use.
 *
 * The issue is that with the current validators we cannot specify
 * arbitrary sets of numbers that are the only allowed parameters.
 * For example the IPARM(2) parameter can take only the values 0, 2,
 * and 3.  The EnhancedNumberValidator can take a min value, and max
 * value, and a step size, but with those options there is no way to
 * specify the needed set.
 *
 * Another missing feature is the ability to give docstrings for such
 * numbers.  For example IPARM(25) can take on the values 0 and 1.
 * This would be easy enough to accomplish with just a number
 * validator, but then have no way to document the effect of each
 * value.
 */
template <class Matrix, class Vector>
Teuchos::RCP<const Teuchos::ParameterList>
PardisoMKL<Matrix,Vector>::getValidParameters_impl() const
{
  using std::string;
  using Teuchos::as;
  using Teuchos::RCP;
  using Teuchos::tuple;
  using Teuchos::toString;
  using Teuchos::EnhancedNumberValidator;
  using Teuchos::setStringToIntegralParameter;
  using Teuchos::anyNumberParameterEntryValidator;
 
  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
 
  if( is_null(valid_params) ){
    Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
 
    // Use pardisoinit to get some default values;
    void *pt_dummy[64];
    PMKL::_INTEGER_t mtype_temp = mtype_;
    PMKL::_INTEGER_t iparm_temp[64];
    PMKL::pardisoinit(pt_dummy,
                      const_cast<PMKL::_INTEGER_t*>(&mtype_temp),
                      const_cast<PMKL::_INTEGER_t*>(iparm_temp));
 
    setStringToIntegralParameter<int>("IPARM(2)", toString(iparm_temp[1]),
                                      "Fill-in reducing ordering for the input matrix",
                                      tuple<string>("0", "2", "3"),
                                      tuple<string>("The minimum degree algorithm",
                                      "Nested dissection algorithm from METIS",
                                      "OpenMP parallel nested dissection algorithm"),
                                      tuple<int>(0, 2, 3),
                                      pl.getRawPtr());
    
    Teuchos::RCP<EnhancedNumberValidator<int> > iparm_4_validator
      = Teuchos::rcp( new EnhancedNumberValidator<int>() );
    iparm_4_validator->setMin(0);
    pl->set("IPARM(4)" , as<int>(iparm_temp[3]) , "Preconditioned CGS/CG",
            iparm_4_validator);
 
    setStringToIntegralParameter<int>("IPARM(12)", toString(iparm_temp[11]),
                                      "Solve with transposed or conjugate transposed matrix A",
                                      tuple<string>("0", "1", "2"),
                                      tuple<string>("Non-transposed",
                                      "Conjugate-transposed",
                                      "Transposed"),
                                      tuple<int>(0, 1, 2),
                                      pl.getRawPtr());
 
    setStringToIntegralParameter<int>("IPARM(13)", toString(iparm_temp[12]),
                                      "Use weighted matching",
                                      tuple<string>("0", "1"),
                                      tuple<string>("No matching", "Use matching"),
                                      tuple<int>(0, 1),
                                      pl.getRawPtr());
 
    setStringToIntegralParameter<int>("IPARM(24)", toString(iparm_temp[23]),
                                      "Parallel factorization control",
                                      tuple<string>("0", "1"),
                                      tuple<string>("PARDISO uses the previous algorithm for factorization",
                                      "PARDISO uses the new two-level factorization algorithm"),
                                      tuple<int>(0, 1),
                                      pl.getRawPtr());
 
    setStringToIntegralParameter<int>("IPARM(25)", toString(iparm_temp[24]),
                                      "Parallel forward/backward solve control",
                                      tuple<string>("0", "1"),
                                      tuple<string>("PARDISO uses the parallel algorithm for the solve step",
                                      "PARDISO uses the sequential forward and backward solve"),
                                      tuple<int>(0, 1),
                                      pl.getRawPtr());
 
    setStringToIntegralParameter<int>("IPARM(60)", toString(iparm_temp[59]),
                                      "PARDISO mode (OOC mode)",
                                      tuple<string>("0", "2"),
                                      tuple<string>("In-core PARDISO",
                                      "Out-of-core PARDISO.  The OOC PARDISO can solve very "
                                      "large problems by holding the matrix factors in files "
                                      "on the disk. Hence the amount of RAM required by OOC "
                                      "PARDISO is significantly reduced."),
                                      tuple<int>(0, 2),
                                      pl.getRawPtr());
 
    Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
      Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
 
    Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int( false );
    accept_int.allowInt( true );
 
    pl->set("IPARM(8)" , as<int>(iparm_temp[7]) , "Iterative refinement step",
            anyNumberParameterEntryValidator(preferred_int, accept_int));
 
    pl->set("IPARM(10)", as<int>(iparm_temp[9]) , "Pivoting perturbation",
            anyNumberParameterEntryValidator(preferred_int, accept_int));
    pl->set("IPARM(11)", as<int>(iparm_temp[10]) , "Scaling vectors",
            anyNumberParameterEntryValidator(preferred_int, accept_int));
 
    pl->set("IPARM(18)", as<int>(iparm_temp[17]), "Report the number of non-zero elements in the factors",
            anyNumberParameterEntryValidator(preferred_int, accept_int));
 
    pl->set("IPARM(27)", as<int>(iparm_temp[26]) , "Check input matrix",
            anyNumberParameterEntryValidator(preferred_int, accept_int));
 
    scalar_type *dummy_scalar_ptr;
    const local_ordinal_type *dummy_ordinal_ptr;
    pl->set("PartialFacto", 0,
            "Perform partial factorization to extract dense Schur complement (0: no, 1: form + factor Schur, 2: ony form");
    pl->set("SchurPart", dummy_ordinal_ptr,
            "Specify rows/columns belonging to Schur complement for partial factorization");
    pl->set("SchurOut", dummy_scalar_ptr,
            "Store output Schur complement from partial factorization");
    pl->set("OnlyForwardSolve", false,
            "Perform only the forward substitution");
    pl->set("OnlyBackwardSolve", false,
            "Perform only the backward substitution");
 
    pl->set("IsContiguous", true, "Whether GIDs contiguous");
    pl->set("MessageLevel", 0, "PardisoMKL message level (0 to turn off message, and 1 to turn on message");
    pl->set("verbose", false, "Set PardisoMKL message level to be 1");
    pl->set("DebugLevel", 0, "Debug message level (0 for no message, and >0 for more message");
 
    valid_params = pl;
  }
 
  return valid_params;
}
 
 
 
template <class Matrix, class Vector>
bool
PardisoMKL<Matrix,Vector>::loadA_impl(EPhase current_phase)
{
#ifdef HAVE_AMESOS2_TIMERS
  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
#endif
  if (debug_level_ > 0) {
    if (this->root_) {
      std::cout << "\n == Amesos2_PardisoMKL::loadA_impl";
      if (current_phase == PREORDERING) std::cout << "(PreOrder)";
      if (current_phase == SYMBFACT) std::cout << "(SymFact)";
      if (current_phase == NUMFACT)  std::cout << "(NumFact)";
      std::cout << " ==" << std::endl;
    }
    Teuchos::RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    this->matrixA_->describe(*fancy, (debug_level_ == 1 ? Teuchos::VERB_LOW : Teuchos::VERB_EXTREME));
  }
 
  // PardisoMKL does not need matrix data in the pre-ordering phase
  if( current_phase == PREORDERING ) return( false );
 
  if( this->root_ ){
    Kokkos::resize(nzvals_view_, this->globalNumNonZeros_);
    Kokkos::resize(colind_view_, this->globalNumNonZeros_);
    Kokkos::resize(rowptr_view_, this->globalNumRows_ + 1);
  }
  {
#ifdef HAVE_AMESOS2_TIMERS
    Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
#endif
 
    int_t nnz_ret = 0;
    Util::get_crs_helper_kokkos_view<
      MatrixAdapter<Matrix>,
      host_value_type_array, host_ordinal_type_array, host_size_type_array>::do_get(
        this->matrixA_.ptr(),
        nzvals_view_, colind_view_, rowptr_view_, nnz_ret, 
        (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
        SORTED_INDICES,
        this->rowIndexBase_);
  }
 
  return( true );
}
 
 
template <class Matrix, class Vector>
void
PardisoMKL<Matrix,Vector>::describe_impl(Teuchos::FancyOStream &out,
                                         const Teuchos::EVerbosityLevel verbLevel) const
{
  out << " PardisoMKL current parameters:" << std::endl;
  out << "  > IPARM(2)  = " << iparm_[1]  << std::endl;
  out << "  > IPARM(4)  = " << iparm_[3]  << std::endl;
  out << "  > IPARM(8)  = " << iparm_[7]  << std::endl;
  out << "  > IPARM(10) = " << iparm_[9]  << std::endl;
  out << "  > IPARM(11) = " << iparm_[10] << std::endl;
  out << "  > IPARM(12) = " << iparm_[11] << std::endl;
  out << "  > IPARM(13) = " << iparm_[12] << std::endl;
  out << "  > IPARM(18) = " << iparm_[17] << std::endl;
  out << "  > IPARM(24) = " << iparm_[23] << std::endl;
  out << "  > IPARM(25) = " << iparm_[24] << std::endl;
  out << "  > IPARM(27) = " << iparm_[26] << std::endl;
  out << "  > IPARM(60) = " << iparm_[59] << std::endl;
  out << "  > PartialFacto = " << partial_facto_;
  if (partial_facto_ == 0)
    out << " (no partial factorization)" << std::endl;
  if (partial_facto_ == 1)
    out << " (compute and factor schur complement : iparam_[35] = 2)" << std::endl;
  if (partial_facto_ == 2)
    out << " (only compute schur complement : iparam_[35] = 1)" << std::endl;
  out << "  > IsContiguous = " << (is_contiguous_ ? "YES" : "NO") << std::endl;
  out << "  > MessageLevel = " << msglvl_ << std::endl;
  out << "  > DebugLevel   = " << debug_level_ << std::endl;
  out << std::endl;
}
 
 
template <class Matrix, class Vector>
void
PardisoMKL<Matrix,Vector>::check_pardiso_mkl_error(EPhase phase,
                                                   int_t error) const
{
  int error_i = error;
  Teuchos::broadcast(*(this->getComm()), 0, &error_i); // We only care about root's value
 
  if( error_i == 0 ) return;      // No error
 
  std::string errmsg = "Other error";
  switch( error_i ){
  case -1:
    errmsg = "PardisoMKL reported error: 'Input inconsistent'";
    break;
  case -2:
    errmsg = "PardisoMKL reported error: 'Not enough memory'";
    break;
  case -3:
    errmsg = "PardisoMKL reported error: 'Reordering problem'";
    break;
  case -4:
    errmsg =
      "PardisoMKL reported error: 'Zero pivot, numerical "
      "factorization or iterative refinement problem'";
    break;
  case -5:
    errmsg = "PardisoMKL reported error: 'Unclassified (internal) error'";
    break;
  case -6:
    errmsg = "PardisoMKL reported error: 'Reordering failed'";
    break;
  case -7:
    errmsg = "PardisoMKL reported error: 'Diagonal matrix is singular'";
    break;
  case -8:
    errmsg = "PardisoMKL reported error: '32-bit integer overflow problem'";
    break;
  case -9:
    errmsg = "PardisoMKL reported error: 'Not enough memory for OOC'";
    break;
  case -10:
    errmsg = "PardisoMKL reported error: 'Problems with opening OOC temporary files'";
    break;
  case -11:
    errmsg = "PardisoMKL reported error: 'Read/write problem with OOC data file'";
    break;
  }
 
  TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error, errmsg );
}
 
 
template <class Matrix, class Vector>
void
PardisoMKL<Matrix,Vector>::set_pardiso_mkl_matrix_type(int_t mtype)
{
  if( mtype == 0 ){
    if( complex_ ){
      mtype_ = 13;              // complex, unsymmetric
    } else {
      mtype_ = 11;              // real, unsymmetric
    }
  } else {
    switch( mtype ){
    case 11:
      TEUCHOS_TEST_FOR_EXCEPTION( complex_,
                          std::invalid_argument,
                          "Cannot set a real Pardiso matrix type with scalar type complex" );
      mtype_ = 11; break;
    case 13:
      TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
                          std::invalid_argument,
                          "Cannot set a complex Pardiso matrix type with non-complex scalars" );
      mtype_ = 13; break;
    default:
      TEUCHOS_TEST_FOR_EXCEPTION( true,
                          std::invalid_argument,
                          "Symmetric matrices are not yet supported by the Amesos2 interface" );
    }
  }
}
 
 
template <class Matrix, class Vector>
const char* PardisoMKL<Matrix,Vector>::name = "PARDISOMKL";
 
template <class Matrix, class Vector>
const typename PardisoMKL<Matrix,Vector>::int_t
PardisoMKL<Matrix,Vector>::maxfct_ = 1;
 
template <class Matrix, class Vector>
const typename PardisoMKL<Matrix,Vector>::int_t
PardisoMKL<Matrix,Vector>::mnum_ = 1;
 
} // end namespace Amesos
 
#endif  // AMESOS2_PARDISOMKL_DEF_HPP