Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Superlumt_def.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 
00052 #ifndef AMESOS2_SUPERLUMT_DEF_HPP
00053 #define AMESOS2_SUPERLUMT_DEF_HPP
00054 
00055 #include <Teuchos_Tuple.hpp>
00056 #include <Teuchos_ParameterList.hpp>
00057 #include <Teuchos_StandardParameterEntryValidators.hpp>
00058 
00059 #include "Amesos2_SolverCore_def.hpp"
00060 #include "Amesos2_Util.hpp"
00061 
00062 
00063 namespace SLUMT {
00064   /*
00065    * We have to declare this extern function because of a bug in the
00066    * SuperLU_MT header files.  In each header is declared a function
00067    * "extern int xsp_colorder(...)" where x is in {'s','d','c','z'},
00068    * but this function is never defined anywhere, so if you use the
00069    * function as defined, it will compile fine, but will choke during
00070    * linking.  No other code in SuperLU_MT actually calls these
00071    * functions. Instead, all other SuperLU_MT function just call
00072    * "sp_colorder", which is defined within the SuperLU_MT library,
00073    * but not declared.
00074    */
00075   extern "C" {
00076     int sp_colorder(SuperMatrix*,int*,superlumt_options_t*,SuperMatrix*);
00077   }
00078 } // end namespace SLUMT
00079 
00080 
00081 namespace Amesos2 {
00082 
00083   /*
00084    * Note: Although many of the type definitions for SuperLU_MT are
00085    * identical to those of SuperLU, we do not mix the definitions so
00086    * that we do not introduce messy coupling between the two
00087    * interfaces.  Also, there exist enough differences between the two
00088    * packages to merit dedicated utility functions.
00089    *
00090    * We have also taken a different approach to interfacing with
00091    * SuperLU_MT than with SuperLU which I believe leads to a better
00092    * seperation of the 4 solution steps.  We may in the future adopt a
00093    * similar strategy for SuperLU.
00094    */
00095 
00096   template <class Matrix, class Vector>
00097   Superlumt<Matrix,Vector>::Superlumt(Teuchos::RCP<const Matrix> A,
00098                                       Teuchos::RCP<Vector>       X,
00099                                       Teuchos::RCP<const Vector> B )
00100     : SolverCore<Amesos2::Superlumt,Matrix,Vector>(A, X, B)
00101     , nzvals_()                 // initialize to empty arrays
00102     , rowind_()
00103     , colptr_()
00104   {
00105     Teuchos::RCP<Teuchos::ParameterList> default_params
00106       = Teuchos::parameterList( *(this->getValidParameters()) );
00107     this->setParameters(default_params);
00108 
00109     data_.options.lwork = 0;    // Use system memory for factorization
00110 
00111     data_.perm_r.resize(this->globalNumRows_);
00112     data_.perm_c.resize(this->globalNumCols_);
00113     data_.options.perm_r = data_.perm_r.getRawPtr();
00114     data_.options.perm_c = data_.perm_c.getRawPtr();
00115 
00116     data_.R.resize(this->globalNumRows_);
00117     data_.C.resize(this->globalNumCols_);
00118 
00119     data_.options.refact = SLUMT::NO; // initially we are not refactoring
00120     data_.equed = SLUMT::NOEQUIL;       // No equilibration has yet been performed
00121     data_.rowequ = false;
00122     data_.colequ = false;
00123 
00124     data_.A.Store  = NULL;
00125     data_.AC.Store = NULL;
00126     data_.BX.Store = NULL;
00127     data_.L.Store  = NULL;
00128     data_.U.Store  = NULL;
00129 
00130     data_.stat.ops = NULL;
00131   }
00132 
00133 
00134   template <class Matrix, class Vector>
00135   Superlumt<Matrix,Vector>::~Superlumt( )
00136   {
00137     /* Free SuperLU_MT data_types
00138      * - Matrices
00139      * - Vectors
00140      * - Stat object
00141      */
00142 
00143     // Storage is initialized in numericFactorization_impl()
00144     if ( data_.A.Store != NULL ){
00145       // Our Teuchos::Array's will destroy rowind, colptr, and nzval for us
00146       SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
00147     }
00148 
00149     // Cleanup memory allocated during last call to sp_colorder if needed
00150     if( data_.AC.Store != NULL ){
00151       SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
00152     }
00153 
00154     if ( data_.L.Store != NULL ){ // will only ever be true for this->root_
00155       SLUMT::D::Destroy_SuperNode_SCP( &(data_.L) );
00156       SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
00157 
00158       // memory alloc'd in sp_colorder
00159       free( data_.options.etree );
00160       free( data_.options.colcnt_h );
00161       free( data_.options.part_super_h );
00162     }
00163 
00164 
00165     // Storage is initialized in solve_impl()
00166     if ( data_.BX.Store != NULL ){
00167       /* Cannot use SLU::Destroy_Dense_Matrix routine here, since it attempts to
00168        * free the array of non-zero values, but that array has already been
00169        * deallocated by the MultiVector object.  So we release just the Store
00170        * instead.
00171        */
00172       SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
00173     }
00174     
00175     if ( data_.stat.ops != NULL )
00176       SLUMT::D::StatFree( &(data_.stat) );
00177   }
00178 
00179   template<class Matrix, class Vector>
00180   int
00181   Superlumt<Matrix,Vector>::preOrdering_impl()
00182   {
00183     // Use either the column-ordering found in the users perm_c or the requested computed ordering
00184     int perm_spec = data_.options.ColPerm;
00185     if( perm_spec != SLUMT::MY_PERMC && this->root_ ){
00186 #ifdef HAVE_AMESOS2_TIMERS
00187       Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
00188 #endif
00189 
00190       SLUMT::S::get_perm_c(perm_spec, &(data_.A), data_.perm_c.getRawPtr());
00191     }
00192     // Ordering will be applied directly before numeric factorization
00193 
00194     return(0);
00195   }
00196 
00197 
00198 
00199   template <class Matrix, class Vector>
00200   int
00201   Superlumt<Matrix,Vector>::symbolicFactorization_impl()
00202   {
00203     // We assume that a call to symbolicFactorization is indicating that
00204     // the structure of the matrix has changed.  SuperLU doesn't allow
00205     // us to do a symbolic factorization directly, but we can leave a
00206     // flag that tells it it needs to symbolically factor later.
00207     data_.options.refact = SLUMT::NO;
00208 
00209     if( this->status_.numericFactorizationDone() && this->root_ ){
00210       // If we've done a numeric factorization already, then we need to
00211       // cleanup the old L and U. Stores and other data will be
00212       // allocated during numeric factorization.  Only rank 0 has valid
00213       // pointers
00214       SLUMT::D::Destroy_SuperNode_Matrix( &(data_.L) );
00215       SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
00216       data_.L.Store = NULL;
00217       data_.U.Store = NULL;
00218     }
00219 
00220     return(0);
00221   }
00222 
00223 
00224   template <class Matrix, class Vector>
00225   int
00226   Superlumt<Matrix,Vector>::numericFactorization_impl()
00227   {
00228     using Teuchos::as;
00229 
00230 #ifdef HAVE_AMESOS2_DEBUG
00231     const int nprocs = data_.options.nprocs;
00232     TEUCHOS_TEST_FOR_EXCEPTION( nprocs <= 0,
00233                         std::invalid_argument,
00234                         "The number of threads to spawn should be greater than 0." );
00235 #endif
00236 
00237     int info = 0;
00238 
00239     if ( this->root_ ) {
00240 
00241       if( data_.options.fact == SLUMT::EQUILIBRATE ){
00242         magnitude_type rowcnd, colcnd, amax;
00243         int info;
00244 
00245         function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
00246                             data_.C.getRawPtr(), &rowcnd, &colcnd,
00247                             &amax, &info);
00248         TEUCHOS_TEST_FOR_EXCEPTION( info != 0,
00249                             std::runtime_error,
00250                             "SuperLU_MT gsequ returned with status " << info );
00251 
00252         function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
00253                             data_.C.getRawPtr(), rowcnd, colcnd,
00254                             amax, &(data_.equed));
00255 
00256         data_.rowequ = (data_.equed == SLUMT::ROW) || (data_.equed == SLUMT::BOTH);
00257         data_.colequ = (data_.equed == SLUMT::COL) || (data_.equed == SLUMT::BOTH);
00258 
00259         data_.options.fact = SLUMT::DOFACT;
00260       }
00261 
00262       // Cleanup memory allocated during last call to sp_colorder if needed
00263       if( data_.AC.Store != NULL ){
00264         SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
00265         if( data_.options.refact == SLUMT::NO ){                 // then we recompute etree; free the old one
00266           free( data_.options.etree );
00267           free( data_.options.colcnt_h );
00268           free( data_.options.part_super_h );
00269         }
00270         data_.AC.Store = NULL;
00271       }
00272 
00273       // Apply the column ordering, so that AC is the column-permuted A, and compute etree
00274       SLUMT::sp_colorder(&(data_.A), data_.perm_c.getRawPtr(),
00275                          &(data_.options), &(data_.AC));
00276 
00277 
00278       // Allocate and initialize status variable
00279       const int n = as<int>(this->globalNumCols_); // n is the number of columns in A
00280       if( data_.stat.ops != NULL ){ SLUMT::D::StatFree( &(data_.stat) ); data_.stat.ops = NULL; }
00281       SLUMT::D::StatAlloc(n, data_.options.nprocs, data_.options.panel_size,
00282                           data_.options.relax, &(data_.stat));
00283       SLUMT::D::StatInit(n, data_.options.nprocs, &(data_.stat));
00284 
00285 
00286       { // Do factorization
00287 #ifdef HAVE_AMESOS2_TIMERS
00288         Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
00289 #endif
00290 
00291 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
00292         std::cout << "SuperLU_MT:: Before numeric factorization" << std::endl;
00293         std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
00294         std::cout << "rowind_ : " << rowind_.toString() << std::endl;
00295         std::cout << "colptr_ : " << colptr_.toString() << std::endl;
00296 #endif
00297 
00298         function_map::gstrf(&(data_.options), &(data_.AC),
00299                             data_.perm_r.getRawPtr(), &(data_.L), &(data_.U),
00300                             &(data_.stat), &info);
00301       }
00302 
00303       // Set the number of non-zero values in the L and U factors
00304       this->setNnzLU( as<size_t>(((SLUMT::SCformat*)data_.L.Store)->nnz +
00305                                  ((SLUMT::NCformat*)data_.U.Store)->nnz) );
00306     }
00307 
00308     // Check output
00309     const global_size_type info_st = as<global_size_type>(info);
00310     TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
00311                         std::runtime_error,
00312                         "Factorization complete, but matrix is singular. Division by zero eminent");
00313     TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
00314                         std::runtime_error,
00315                         "Memory allocation failure in SuperLU_MT factorization");
00316     // The other option, that info_st < 0 denotes invalid parameters to
00317     // the function, but we'll assume for now that that won't happen.
00318 
00319     data_.options.fact = SLUMT::FACTORED;
00320     data_.options.refact = SLUMT::YES;
00321 
00322     /* All processes should return the same error code */
00323     Teuchos::broadcast(*(this->getComm()),0,&info);
00324     return(info);
00325   }
00326 
00327 
00328   template <class Matrix, class Vector>
00329   int
00330   Superlumt<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
00331                                        const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
00332   {
00333     using Teuchos::as;
00334 
00335     const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
00336     const size_t nrhs = X->getGlobalNumVectors();
00337 
00338     Teuchos::Array<slu_type> bxvals_(ld_rhs * nrhs);
00339     size_t ldbx_ = as<size_t>(ld_rhs);
00340 
00341     {                           // Get values from B
00342 #ifdef HAVE_AMESOS2_TIMERS
00343       Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
00344       Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
00345 #endif
00346 
00347       Util::get_1d_copy_helper<MultiVecAdapter<Vector>,
00348         slu_type>::do_get(B, bxvals_(),
00349                           ldbx_,
00350                           ROOTED);
00351     }
00352 
00353     int info = 0; // returned error code (0 = success)
00354     if ( this->root_ ) {
00355       // Clean up old B stores if it has already been created
00356       if( data_.BX.Store != NULL ){
00357         SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
00358         data_.BX.Store = NULL;
00359       }
00360 
00361       {
00362 #ifdef HAVE_AMESOS2_TIMERS
00363         Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
00364 #endif
00365         SLUMT::Dtype_t dtype = type_map::dtype;
00366         function_map::create_Dense_Matrix(&(data_.BX), as<int>(ld_rhs), as<int>(nrhs),
00367                                           bxvals_.getRawPtr(), as<int>(ldbx_),
00368                                           SLUMT::SLU_DN, dtype, SLUMT::SLU_GE);
00369       }
00370 
00371       if( data_.options.trans == SLUMT::NOTRANS ){
00372         if( data_.rowequ ){             // row equilibration has been done on AC
00373           // scale bxvals_ by diag(R)
00374           Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
00375                       SLUMT::slu_mt_mult<slu_type,magnitude_type>());
00376         }
00377       } else if( data_.colequ ){        // column equilibration has been done on AC
00378         // scale bxvals_ by diag(C)
00379         Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
00380                     SLUMT::slu_mt_mult<slu_type,magnitude_type>());
00381       }
00382 
00383 
00384       {
00385 #ifdef HAVE_AMESOS2_TIMERS
00386         Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
00387 #endif
00388 
00389         function_map::gstrs(data_.options.trans, &(data_.L),
00390                             &(data_.U), data_.perm_r.getRawPtr(),
00391                             data_.perm_c.getRawPtr(), &(data_.BX),
00392                             &(data_.stat), &info);
00393       }
00394     } // end block for solve time
00395 
00396     /* All processes should have the same error code */
00397     Teuchos::broadcast(*(this->getComm()),0,&info);
00398 
00399     TEUCHOS_TEST_FOR_EXCEPTION( info < 0,
00400                         std::runtime_error,
00401                         "Argument " << -info << " to gstrs had an illegal value" );
00402 
00403     // "Un-scale" the solution so that it is a solution of the original system
00404     if( data_.options.trans == SLUMT::NOTRANS ){
00405       if( data_.colequ ){       // column equilibration has been done on AC
00406         // scale bxvals_ by diag(C)
00407         Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
00408                     SLUMT::slu_mt_mult<slu_type,magnitude_type>());
00409       }
00410     } else if( data_.rowequ ){          // row equilibration has been done on AC
00411       // scale bxvals_ by diag(R)
00412       Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
00413                   SLUMT::slu_mt_mult<slu_type,magnitude_type>());
00414     }
00415 
00416     /* Update X's global values */
00417     {
00418 #ifdef HAVE_AMESOS2_TIMERS
00419       Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
00420 #endif
00421 
00422       Util::put_1d_data_helper<
00423         MultiVecAdapter<Vector>, slu_type>::do_put(X, bxvals_(), ldbx_, ROOTED);
00424     }
00425 
00426     return(info);
00427   }
00428 
00429 
00430   template <class Matrix, class Vector>
00431   bool
00432   Superlumt<Matrix,Vector>::matrixShapeOK_impl() const
00433   {
00434     // The SuperLU_MT factorization routines can handle square as well as
00435     // rectangular matrices, but SuperLU_MT can only apply the solve routines to
00436     // square matrices, so we check the matrix for squareness.
00437     return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
00438   }
00439 
00440 
00441   template <class Matrix, class Vector>
00442   void
00443   Superlumt<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
00444   {
00445     using Teuchos::as;
00446     using Teuchos::RCP;
00447     using Teuchos::getIntegralValue;
00448     using Teuchos::ParameterEntryValidator;
00449 
00450     RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
00451 
00452 
00453     data_.options.nprocs = parameterList->get<int>("nprocs", 1);
00454 
00455     data_.options.trans = this->control_.useTranspose_ ? SLUMT::TRANS : SLUMT::NOTRANS;
00456     // SuperLU_MT "trans" option can override the Amesos2 option
00457     if( parameterList->isParameter("trans") ){
00458       RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("trans").validator();
00459       parameterList->getEntry("trans").setValidator(trans_validator);
00460 
00461       data_.options.trans = getIntegralValue<SLUMT::trans_t>(*parameterList, "trans");
00462     }
00463 
00464     data_.options.panel_size = parameterList->get<int>("panel_size", SLUMT::D::sp_ienv(1));
00465 
00466     data_.options.relax = parameterList->get<int>("relax", SLUMT::D::sp_ienv(2));
00467 
00468     const bool equil = parameterList->get<bool>("Equil", true);
00469     data_.options.fact = equil ? SLUMT::EQUILIBRATE : SLUMT::DOFACT;
00470 
00471     const bool symmetric_mode = parameterList->get<bool>("SymmetricMode", false);
00472     data_.options.SymmetricMode = symmetric_mode ? SLUMT::YES : SLUMT::NO;
00473 
00474     const bool print_stat = parameterList->get<bool>("PrintStat", false);
00475     data_.options.PrintStat = print_stat ? SLUMT::YES : SLUMT::NO;
00476 
00477     data_.options.diag_pivot_thresh = parameterList->get<double>("diag_pivot_thresh", 1.0);
00478 
00479     if( parameterList->isParameter("ColPerm") ){
00480       RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
00481       parameterList->getEntry("ColPerm").setValidator(colperm_validator);
00482 
00483       data_.options.ColPerm = getIntegralValue<SLUMT::colperm_t>(*parameterList, "ColPerm");
00484     }
00485 
00486     // TODO: until we support retrieving col/row permutations from the user
00487     data_.options.usepr = SLUMT::NO;
00488   }
00489 
00490 
00491   template <class Matrix, class Vector>
00492   Teuchos::RCP<const Teuchos::ParameterList>
00493   Superlumt<Matrix,Vector>::getValidParameters_impl() const
00494   {
00495     using std::string;
00496     using Teuchos::tuple;
00497     using Teuchos::ParameterList;
00498     using Teuchos::EnhancedNumberValidator;
00499     using Teuchos::setStringToIntegralParameter;
00500     using Teuchos::stringToIntegralParameterEntryValidator;
00501 
00502     static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
00503 
00504     if( is_null(valid_params) ){
00505       Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00506 
00507       Teuchos::RCP<EnhancedNumberValidator<int> > nprocs_validator
00508         = Teuchos::rcp( new EnhancedNumberValidator<int>() );
00509       nprocs_validator->setMin(1);
00510       pl->set("nprocs", 1, "The number of processors to factorize with", nprocs_validator);
00511 
00512       setStringToIntegralParameter<SLUMT::trans_t>("trans", "NOTRANS",
00513                                                    "Solve for the transpose system or not",
00514                                                    tuple<string>("TRANS","NOTRANS","CONJ"),
00515                                                    tuple<string>("Solve with transpose",
00516                                                                  "Do not solve with transpose",
00517                                                                  "Solve with the conjugate transpose"),
00518                                                    tuple<SLUMT::trans_t>(SLUMT::TRANS,
00519                                                                          SLUMT::NOTRANS,
00520                                                                          SLUMT::CONJ),
00521                                                    pl.getRawPtr());
00522 
00523       Teuchos::RCP<EnhancedNumberValidator<int> > panel_size_validator
00524         = Teuchos::rcp( new EnhancedNumberValidator<int>() );
00525       panel_size_validator->setMin(1);
00526       pl->set("panel_size", SLUMT::D::sp_ienv(1),
00527               "Specifies the number of consecutive columns to be treated as a unit of task.",
00528               panel_size_validator);
00529 
00530       Teuchos::RCP<EnhancedNumberValidator<int> > relax_validator
00531         = Teuchos::rcp( new EnhancedNumberValidator<int>() );
00532       relax_validator->setMin(1);
00533       pl->set("relax", SLUMT::D::sp_ienv(2),
00534               "Specifies the number of columns to be grouped as a relaxed supernode",
00535               relax_validator);
00536 
00537       pl->set("Equil", true, "Whether to equilibrate the system before solve");
00538 
00539       Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
00540         = Teuchos::rcp( new EnhancedNumberValidator<double>(0.0, 1.0) );
00541       pl->set("diag_pivot_thresh", 1.0,
00542               "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
00543               diag_pivot_thresh_validator); // partial pivoting
00544 
00545       // Note: MY_PERMC not yet supported
00546       setStringToIntegralParameter<SLUMT::colperm_t>("ColPerm", "COLAMD",
00547                                                      "Specifies how to permute the columns of the "
00548                                                      "matrix for sparsity preservation",
00549                                                      tuple<string>("NATURAL",
00550                                                                    "MMD_AT_PLUS_A",
00551                                                                    "MMD_ATA",
00552                                                                    "COLAMD"),
00553                                                      tuple<string>("Natural ordering",
00554                                                                        "Minimum degree ordering on A^T + A",
00555                                                                    "Minimum degree ordering on A^T A",
00556                                                                    "Approximate minimum degree column ordering"),
00557                                                      tuple<SLUMT::colperm_t>(SLUMT::NATURAL,
00558                                                                              SLUMT::MMD_AT_PLUS_A,
00559                                                                              SLUMT::MMD_ATA,
00560                                                                              SLUMT::COLAMD),
00561                                                      pl.getRawPtr());
00562 
00563       pl->set("SymmetricMode", false, "Specifies whether to use the symmetric mode");
00564 
00565       // TODO: until we support getting row/col permutations from user
00566       //      pl->set("usepr", false);
00567 
00568       pl->set("PrintStat", false, "Specifies whether to print the solver's statistics");
00569 
00570       valid_params = pl;
00571     }
00572 
00573     return valid_params;
00574   }
00575 
00576 
00577   template <class Matrix, class Vector>
00578   bool
00579   Superlumt<Matrix,Vector>::loadA_impl(EPhase current_phase)
00580   {
00581     using Teuchos::as;
00582 
00583 #ifdef HAVE_AMESOS2_TIMERS
00584     Teuchos::TimeMonitor convTimer( this->timers_.mtxConvTime_ );
00585 #endif
00586 
00587     if( current_phase == SYMBFACT ) return false;
00588 
00589     // Store is allocated on create_CompCol_Matrix
00590     if( data_.A.Store != NULL ){
00591       SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
00592       data_.A.Store = NULL;
00593     }
00594 
00595     if( this->root_ ){
00596       nzvals_.resize(this->globalNumNonZeros_);
00597       rowind_.resize(this->globalNumNonZeros_);
00598       colptr_.resize(this->globalNumCols_ + 1);
00599     }
00600 
00601     int nnz_ret = 0;
00602     {
00603 #ifdef HAVE_AMESOS2_TIMERS
00604       Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
00605 #endif
00606 
00607       // Use compressed-column store for SuperLU_MT
00608       Util::get_ccs_helper<
00609       MatrixAdapter<Matrix>,slu_type,int,int>::do_get(this->matrixA_.ptr(),
00610                                                       nzvals_, rowind_, colptr_,
00611                                                       nnz_ret, ROOTED, ARBITRARY);
00612     }
00613 
00614     if( this->root_ ){
00615       TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
00616                           std::runtime_error,
00617                           "rank=0 failed to get all nonzero values in getCcs()");
00618 
00619       SLUMT::Dtype_t dtype = type_map::dtype;
00620       function_map::create_CompCol_Matrix(&(data_.A),
00621                                           as<int>(this->globalNumRows_),
00622                                           as<int>(this->globalNumCols_),
00623                                           nnz_ret,
00624                                           nzvals_.getRawPtr(),
00625                                           rowind_.getRawPtr(),
00626                                           colptr_.getRawPtr(),
00627                                           SLUMT::SLU_NC,
00628                                           dtype, SLUMT::SLU_GE);
00629 
00630       TEUCHOS_TEST_FOR_EXCEPTION( data_.A.Store == NULL,
00631                           std::runtime_error,
00632                           "Failed to allocate matrix store" );
00633     }
00634 
00635     return true;
00636   }
00637 
00638 
00639   template<class Matrix, class Vector>
00640   const char* Superlumt<Matrix,Vector>::name = "SuperLU_MT";
00641 
00642 
00643 } // end namespace Amesos2
00644 
00645 #endif  // AMESOS2_SUPERLUMT_DEF_HPP