Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Superlu_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 
00053 #ifndef AMESOS2_SUPERLU_DEF_HPP
00054 #define AMESOS2_SUPERLU_DEF_HPP
00055 
00056 #include <Teuchos_Tuple.hpp>
00057 #include <Teuchos_ParameterList.hpp>
00058 #include <Teuchos_StandardParameterEntryValidators.hpp>
00059 
00060 #include "Amesos2_SolverCore_def.hpp"
00061 #include "Amesos2_Superlu_decl.hpp"
00062 
00063 namespace Amesos2 {
00064 
00065 
00066 template <class Matrix, class Vector>
00067 Superlu<Matrix,Vector>::Superlu(
00068   Teuchos::RCP<const Matrix> A,
00069   Teuchos::RCP<Vector>       X,
00070   Teuchos::RCP<const Vector> B )
00071   : SolverCore<Amesos2::Superlu,Matrix,Vector>(A, X, B)
00072   , nzvals_()                   // initialize to empty arrays
00073   , rowind_()
00074   , colptr_()
00075 {
00076   // ilu_set_default_options is called later in set parameter list if required.
00077   // This is not the ideal way, but the other option to always call
00078   // ilu_set_default_options here and assuming it won't have any side effect
00079   // in the TPL is more dangerous. It is not a good idea to rely on external
00080   // libraries' internal "features".
00081   SLU::set_default_options(&(data_.options));
00082   // Override some default options
00083   data_.options.PrintStat = SLU::NO;
00084 
00085   SLU::StatInit(&(data_.stat));
00086 
00087   data_.perm_r.resize(this->globalNumRows_);
00088   data_.perm_c.resize(this->globalNumCols_);
00089   data_.etree.resize(this->globalNumCols_);
00090   data_.R.resize(this->globalNumRows_);
00091   data_.C.resize(this->globalNumCols_);
00092 
00093   data_.relax = SLU::sp_ienv(2); // Query optimal relax param from superlu
00094   data_.panel_size = SLU::sp_ienv(1); // Query optimal panel size
00095 
00096   data_.equed = 'N';            // No equilibration
00097   data_.A.Store = NULL;
00098   data_.L.Store = NULL;
00099   data_.U.Store = NULL;
00100   data_.X.Store = NULL;
00101   data_.B.Store = NULL;
00102   
00103   ILU_Flag_=false; // default: turn off ILU
00104 }
00105 
00106 
00107 template <class Matrix, class Vector>
00108 Superlu<Matrix,Vector>::~Superlu( )
00109 {
00110   /* Free Superlu data_types
00111    * - Matrices
00112    * - Vectors
00113    * - Stat object
00114    */
00115   SLU::StatFree( &(data_.stat) ) ;
00116 
00117   // Storage is initialized in numericFactorization_impl()
00118   if ( data_.A.Store != NULL ){
00119     SLU::Destroy_SuperMatrix_Store( &(data_.A) );
00120   }
00121 
00122   // only root allocated these SuperMatrices.
00123   if ( data_.L.Store != NULL ){ // will only be true for this->root_
00124     SLU::Destroy_SuperNode_Matrix( &(data_.L) );
00125     SLU::Destroy_CompCol_Matrix( &(data_.U) );
00126   }
00127 }
00128 
00129 template<class Matrix, class Vector>
00130 int
00131 Superlu<Matrix,Vector>::preOrdering_impl()
00132 {
00133   /*
00134    * Get column permutation vector perm_c[], according to permc_spec:
00135    *   permc_spec = NATURAL:  natural ordering
00136    *   permc_spec = MMD_AT_PLUS_A: minimum degree on structure of A'+A
00137    *   permc_spec = MMD_ATA:  minimum degree on structure of A'*A
00138    *   permc_spec = COLAMD:   approximate minimum degree column ordering
00139    *   permc_spec = MY_PERMC: the ordering already supplied in perm_c[]
00140    */
00141   int permc_spec = data_.options.ColPerm;
00142   if ( permc_spec != SLU::MY_PERMC && this->root_ ){
00143 #ifdef HAVE_AMESOS2_TIMERS
00144     Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
00145 #endif
00146 
00147     SLU::get_perm_c(permc_spec, &(data_.A), data_.perm_c.getRawPtr());
00148   }
00149 
00150   return(0);
00151 }
00152 
00153 
00154 template <class Matrix, class Vector>
00155 int
00156 Superlu<Matrix,Vector>::symbolicFactorization_impl()
00157 {
00158   /*
00159    * SuperLU performs symbolic factorization and numeric factorization
00160    * together, but does leave some options for reusing symbolic
00161    * structures that have been created on previous factorizations.  If
00162    * our Amesos2 user calls this function, that is an indication that
00163    * the symbolic structure of the matrix is no longer valid, and
00164    * SuperLU should do the factorization from scratch.
00165    *
00166    * This can be accomplished by setting the options.Fact flag to
00167    * DOFACT, as well as setting our own internal flag to false.
00168    */
00169   same_symbolic_ = false;
00170   data_.options.Fact = SLU::DOFACT;
00171 
00172   return(0);
00173 }
00174 
00175 
00176 template <class Matrix, class Vector>
00177 int
00178 Superlu<Matrix,Vector>::numericFactorization_impl()
00179 {
00180   using Teuchos::as;
00181 
00182   // Cleanup old L and U matrices if we are not reusing a symbolic
00183   // factorization.  Stores and other data will be allocated in gstrf.
00184   // Only rank 0 has valid pointers
00185   if ( !same_symbolic_ && data_.L.Store != NULL ){
00186     SLU::Destroy_SuperNode_Matrix( &(data_.L) );
00187     SLU::Destroy_CompCol_Matrix( &(data_.U) );
00188     data_.L.Store = NULL;
00189     data_.U.Store = NULL;
00190   }
00191 
00192   if( same_symbolic_ ) data_.options.Fact = SLU::SamePattern_SameRowPerm;
00193 
00194   int info = 0;
00195   if ( this->root_ ){
00196 
00197 #ifdef HAVE_AMESOS2_DEBUG
00198     TEUCHOS_TEST_FOR_EXCEPTION( data_.A.ncol != as<int>(this->globalNumCols_),
00199                         std::runtime_error,
00200                         "Error in converting to SuperLU SuperMatrix: wrong number of global columns." );
00201     TEUCHOS_TEST_FOR_EXCEPTION( data_.A.nrow != as<int>(this->globalNumRows_),
00202                         std::runtime_error,
00203                         "Error in converting to SuperLU SuperMatrix: wrong number of global rows." );
00204 #endif
00205 
00206     if( data_.options.Equil == SLU::YES ){
00207       magnitude_type rowcnd, colcnd, amax;
00208       int info2 = 0;
00209 
00210       // calculate row and column scalings
00211       function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
00212                           data_.C.getRawPtr(), &rowcnd, &colcnd,
00213                           &amax, &info2);
00214       TEUCHOS_TEST_FOR_EXCEPTION( info2 != 0,
00215                           std::runtime_error,
00216                           "SuperLU gsequ returned with status " << info2 );
00217 
00218       // apply row and column scalings if necessary
00219       function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
00220                           data_.C.getRawPtr(), rowcnd, colcnd,
00221                           amax, &(data_.equed));
00222 
00223       // // check what types of equilibration was actually done
00224       // data_.rowequ = (data_.equed == 'R') || (data_.equed == 'B');
00225       // data_.colequ = (data_.equed == 'C') || (data_.equed == 'B');
00226     }
00227 
00228     // Apply the column permutation computed in preOrdering.  Place the
00229     // column-permuted matrix in AC
00230     SLU::sp_preorder(&(data_.options), &(data_.A), data_.perm_c.getRawPtr(),
00231                      data_.etree.getRawPtr(), &(data_.AC));
00232 
00233     { // Do factorization
00234 #ifdef HAVE_AMESOS2_TIMERS
00235       Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
00236 #endif
00237 
00238 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
00239       std::cout << "Superlu:: Before numeric factorization" << std::endl;
00240       std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
00241       std::cout << "rowind_ : " << rowind_.toString() << std::endl;
00242       std::cout << "colptr_ : " << colptr_.toString() << std::endl;
00243 #endif
00244 
00245       if(ILU_Flag_==false) {
00246         function_map::gstrf(&(data_.options), &(data_.AC),
00247             data_.relax, data_.panel_size, data_.etree.getRawPtr(),
00248             NULL, 0, data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
00249             &(data_.L), &(data_.U), &(data_.stat), &info);
00250       }
00251       else {
00252         function_map::gsitrf(&(data_.options), &(data_.AC),
00253             data_.relax, data_.panel_size, data_.etree.getRawPtr(),
00254             NULL, 0, data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
00255             &(data_.L), &(data_.U), &(data_.stat), &info);
00256       }
00257 
00258     }
00259     // Cleanup. AC data will be alloc'd again for next factorization (if at all)
00260     SLU::Destroy_CompCol_Permuted( &(data_.AC) );
00261 
00262     // Set the number of non-zero values in the L and U factors
00263     this->setNnzLU( as<size_t>(((SLU::SCformat*)data_.L.Store)->nnz +
00264                                ((SLU::NCformat*)data_.U.Store)->nnz) );
00265   }
00266 
00267   /* All processes should have the same error code */
00268   Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
00269 
00270   global_size_type info_st = as<global_size_type>(info);
00271   TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
00272     std::runtime_error,
00273     "Factorization complete, but matrix is singular. Division by zero eminent");
00274   TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
00275     std::runtime_error,
00276     "Memory allocation failure in Superlu factorization");
00277 
00278   data_.options.Fact = SLU::FACTORED;
00279   same_symbolic_ = true;
00280 
00281   return(info);
00282 }
00283 
00284 
00285 template <class Matrix, class Vector>
00286 int
00287 Superlu<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> >       X,
00288                                    const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
00289 {
00290   using Teuchos::as;
00291 
00292   const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
00293   const size_t nrhs = X->getGlobalNumVectors();
00294 
00295   const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
00296   Teuchos::Array<slu_type> xValues(val_store_size);
00297   Teuchos::Array<slu_type> bValues(val_store_size);
00298 
00299   {                             // Get values from RHS B
00300 #ifdef HAVE_AMESOS2_TIMERS
00301     Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
00302     Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
00303 #endif
00304     Util::get_1d_copy_helper<MultiVecAdapter<Vector>,
00305                              slu_type>::do_get(B, bValues(),
00306                                                as<size_t>(ld_rhs),
00307                                                ROOTED, this->rowIndexBase_);
00308   }
00309 
00310   int ierr = 0; // returned error code
00311 
00312   magnitude_type rpg, rcond;
00313   if ( this->root_ ) {
00314     data_.ferr.resize(nrhs);
00315     data_.berr.resize(nrhs);
00316 
00317     {
00318 #ifdef HAVE_AMESOS2_TIMERS
00319       Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
00320 #endif
00321       SLU::Dtype_t dtype = type_map::dtype;
00322       int i_ld_rhs = as<int>(ld_rhs);
00323       function_map::create_Dense_Matrix(&(data_.B), i_ld_rhs, as<int>(nrhs),
00324                                         bValues.getRawPtr(), i_ld_rhs,
00325                                         SLU::SLU_DN, dtype, SLU::SLU_GE);
00326       function_map::create_Dense_Matrix(&(data_.X), i_ld_rhs, as<int>(nrhs),
00327                                         xValues.getRawPtr(), i_ld_rhs,
00328                                         SLU::SLU_DN, dtype, SLU::SLU_GE);
00329     }
00330 
00331     // Note: the values of B and X (after solution) are adjusted
00332     // appropriately within gssvx for row and column scalings.
00333 
00334     {                           // Do solve!
00335 #ifdef HAVE_AMESOS2_TIMERS
00336     Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
00337 #endif
00338 
00339     if(ILU_Flag_==false) {
00340       function_map::gssvx(&(data_.options), &(data_.A),
00341           data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
00342           data_.etree.getRawPtr(), &(data_.equed), data_.R.getRawPtr(),
00343           data_.C.getRawPtr(), &(data_.L), &(data_.U), NULL, 0, &(data_.B),
00344           &(data_.X), &rpg, &rcond, data_.ferr.getRawPtr(),
00345           data_.berr.getRawPtr(), &(data_.mem_usage), &(data_.stat), &ierr);
00346     }
00347     else {
00348       function_map::gsisx(&(data_.options), &(data_.A),
00349           data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
00350           data_.etree.getRawPtr(), &(data_.equed), data_.R.getRawPtr(),
00351           data_.C.getRawPtr(), &(data_.L), &(data_.U), NULL, 0, &(data_.B),
00352           &(data_.X), &rpg, &rcond, &(data_.mem_usage), &(data_.stat), &ierr);
00353     }
00354 
00355     }
00356 
00357     // Cleanup X and B stores
00358     SLU::Destroy_SuperMatrix_Store( &(data_.X) );
00359     SLU::Destroy_SuperMatrix_Store( &(data_.B) );
00360     data_.X.Store = NULL;
00361     data_.B.Store = NULL;
00362   }
00363 
00364   /* All processes should have the same error code */
00365   Teuchos::broadcast(*(this->getComm()), 0, &ierr);
00366 
00367   global_size_type ierr_st = as<global_size_type>(ierr);
00368   TEUCHOS_TEST_FOR_EXCEPTION( ierr < 0,
00369                       std::invalid_argument,
00370                       "Argument " << -ierr << " to SuperLU xgssvx had illegal value" );
00371   TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st <= this->globalNumCols_,
00372                       std::runtime_error,
00373                       "Factorization complete, but U is exactly singular" );
00374   TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st > this->globalNumCols_ + 1,
00375                       std::runtime_error,
00376                       "SuperLU allocated " << ierr - this->globalNumCols_ << " bytes of "
00377                       "memory before allocation failure occured." );
00378 
00379   /* Update X's global values */
00380   {
00381 #ifdef HAVE_AMESOS2_TIMERS
00382     Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
00383 #endif
00384 
00385     Util::put_1d_data_helper<
00386       MultiVecAdapter<Vector>,slu_type>::do_put(X, xValues(),
00387                                          as<size_t>(ld_rhs),
00388                                          ROOTED, this->rowIndexBase_);
00389   }
00390 
00391 
00392   return(ierr);
00393 }
00394 
00395 
00396 template <class Matrix, class Vector>
00397 bool
00398 Superlu<Matrix,Vector>::matrixShapeOK_impl() const
00399 {
00400   // The Superlu factorization routines can handle square as well as
00401   // rectangular matrices, but Superlu can only apply the solve routines to
00402   // square matrices, so we check the matrix for squareness.
00403   return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
00404 }
00405 
00406 
00407 template <class Matrix, class Vector>
00408 void
00409 Superlu<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
00410 {
00411   using Teuchos::RCP;
00412   using Teuchos::getIntegralValue;
00413   using Teuchos::ParameterEntryValidator;
00414 
00415   RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
00416 
00417   ILU_Flag_ = parameterList->get<bool>("ILU_Flag",false);
00418   if (ILU_Flag_) {
00419       SLU::ilu_set_default_options(&(data_.options));
00420       // Override some default options
00421       data_.options.PrintStat = SLU::NO;
00422   }
00423 
00424   data_.options.Trans = this->control_.useTranspose_ ? SLU::TRANS : SLU::NOTRANS;
00425   // The SuperLU transpose option can override the Amesos2 option
00426   if( parameterList->isParameter("Trans") ){
00427     RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("Trans").validator();
00428     parameterList->getEntry("Trans").setValidator(trans_validator);
00429 
00430     data_.options.Trans = getIntegralValue<SLU::trans_t>(*parameterList, "Trans");
00431   }
00432 
00433   if( parameterList->isParameter("IterRefine") ){
00434     RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry("IterRefine").validator();
00435     parameterList->getEntry("IterRefine").setValidator(refine_validator);
00436 
00437     data_.options.IterRefine = getIntegralValue<SLU::IterRefine_t>(*parameterList, "IterRefine");
00438   }
00439 
00440   if( parameterList->isParameter("ColPerm") ){
00441     RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
00442     parameterList->getEntry("ColPerm").setValidator(colperm_validator);
00443 
00444     data_.options.ColPerm = getIntegralValue<SLU::colperm_t>(*parameterList, "ColPerm");
00445   }
00446 
00447   data_.options.DiagPivotThresh = parameterList->get<double>("DiagPivotThresh", 1.0);
00448 
00449   bool equil = parameterList->get<bool>("Equil", true);
00450   data_.options.Equil = equil ? SLU::YES : SLU::NO;
00451 
00452   bool symmetric_mode = parameterList->get<bool>("SymmetricMode", false);
00453   data_.options.SymmetricMode = symmetric_mode ? SLU::YES : SLU::NO;
00454 
00455   // ILU parameters
00456   if( parameterList->isParameter("RowPerm") ){
00457     RCP<const ParameterEntryValidator> rowperm_validator = valid_params->getEntry("RowPerm").validator();
00458     parameterList->getEntry("RowPerm").setValidator(rowperm_validator);
00459     data_.options.RowPerm = getIntegralValue<SLU::rowperm_t>(*parameterList, "RowPerm");
00460   }
00461 
00462   /*if( parameterList->isParameter("ILU_DropRule") ) {
00463     RCP<const ParameterEntryValidator> droprule_validator = valid_params->getEntry("ILU_DropRule").validator();
00464     parameterList->getEntry("ILU_DropRule").setValidator(droprule_validator);
00465     data_.options.ILU_DropRule = getIntegralValue<SLU::rule_t>(*parameterList, "ILU_DropRule");
00466   }*/
00467 
00468   data_.options.ILU_DropTol = parameterList->get<double>("ILU_DropTol", 0.0001);
00469 
00470   data_.options.ILU_FillFactor = parameterList->get<double>("ILU_FillFactor", 10.0);
00471 
00472   if( parameterList->isParameter("ILU_Norm") ) {
00473     RCP<const ParameterEntryValidator> norm_validator = valid_params->getEntry("ILU_Norm").validator();
00474     parameterList->getEntry("ILU_Norm").setValidator(norm_validator);
00475     data_.options.ILU_Norm = getIntegralValue<SLU::norm_t>(*parameterList, "ILU_Norm");
00476   }
00477 
00478   if( parameterList->isParameter("ILU_MILU") ) {
00479     RCP<const ParameterEntryValidator> milu_validator = valid_params->getEntry("ILU_MILU").validator();
00480     parameterList->getEntry("ILU_MILU").setValidator(milu_validator);
00481     data_.options.ILU_MILU = getIntegralValue<SLU::milu_t>(*parameterList, "ILU_MILU");
00482   }
00483 
00484   data_.options.ILU_FillTol = parameterList->get<double>("ILU_FillTol", 0.01);
00485 
00486 }
00487 
00488 
00489 template <class Matrix, class Vector>
00490 Teuchos::RCP<const Teuchos::ParameterList>
00491 Superlu<Matrix,Vector>::getValidParameters_impl() const
00492 {
00493   using std::string;
00494   using Teuchos::tuple;
00495   using Teuchos::ParameterList;
00496   using Teuchos::EnhancedNumberValidator;
00497   using Teuchos::setStringToIntegralParameter;
00498   using Teuchos::stringToIntegralParameterEntryValidator;
00499 
00500   static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
00501 
00502   if( is_null(valid_params) ){
00503     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00504 
00505     setStringToIntegralParameter<SLU::trans_t>("Trans", "NOTRANS",
00506                                                "Solve for the transpose system or not",
00507                                                tuple<string>("TRANS","NOTRANS","CONJ"),
00508                                                tuple<string>("Solve with transpose",
00509                                                              "Do not solve with transpose",
00510                                                              "Solve with the conjugate transpose"),
00511                                                tuple<SLU::trans_t>(SLU::TRANS,
00512                                                                    SLU::NOTRANS,
00513                                                                    SLU::CONJ),
00514                                                pl.getRawPtr());
00515 
00516     setStringToIntegralParameter<SLU::IterRefine_t>("IterRefine", "NOREFINE",
00517                                                     "Type of iterative refinement to use",
00518                                                     tuple<string>("NOREFINE", "SLU_SINGLE", "SLU_DOUBLE"),
00519                                                     tuple<string>("Do not use iterative refinement",
00520                                                                   "Do single iterative refinement",
00521                                                                   "Do double iterative refinement"),
00522                                                     tuple<SLU::IterRefine_t>(SLU::NOREFINE,
00523                                                                              SLU::SLU_SINGLE,
00524                                                                              SLU::SLU_DOUBLE),
00525                                                     pl.getRawPtr());
00526 
00527     // Note: MY_PERMC not yet supported
00528     setStringToIntegralParameter<SLU::colperm_t>("ColPerm", "COLAMD",
00529                                                  "Specifies how to permute the columns of the "
00530                                                  "matrix for sparsity preservation",
00531                                                  tuple<string>("NATURAL", "MMD_AT_PLUS_A",
00532                                                                "MMD_ATA", "COLAMD"),
00533                                                  tuple<string>("Natural ordering",
00534                                                                "Minimum degree ordering on A^T + A",
00535                                                                "Minimum degree ordering on A^T A",
00536                                                                "Approximate minimum degree column ordering"),
00537                                                  tuple<SLU::colperm_t>(SLU::NATURAL,
00538                                                                        SLU::MMD_AT_PLUS_A,
00539                                                                        SLU::MMD_ATA,
00540                                                                        SLU::COLAMD),
00541                                                  pl.getRawPtr());
00542 
00543     Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
00544       = Teuchos::rcp( new EnhancedNumberValidator<double>(0.0, 1.0) );
00545     pl->set("DiagPivotThresh", 1.0,
00546             "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
00547             diag_pivot_thresh_validator); // partial pivoting
00548 
00549     pl->set("Equil", true, "Whether to equilibrate the system before solve");
00550 
00551     pl->set("SymmetricMode", false,
00552             "Specifies whether to use the symmetric mode. "
00553             "Gives preference to diagonal pivots and uses "
00554             "an (A^T + A)-based column permutation.");
00555 
00556     // ILU parameters
00557 
00558     setStringToIntegralParameter<SLU::rowperm_t>("RowPerm", "LargeDiag",
00559             "Type of row permutation strategy to use",
00560             tuple<string>("NOROWPERM","LargeDiag","MY_PERMR"),
00561             tuple<string>("Use natural ordering",
00562             "Use weighted bipartite matching algorithm",
00563             "Use the ordering given in perm_r input"),
00564             tuple<SLU::rowperm_t>(SLU::NOROWPERM,
00565             SLU::LargeDiag,
00566             SLU::MY_PERMR),
00567             pl.getRawPtr());
00568 
00569     /*setStringToIntegralParameter<SLU::rule_t>("ILU_DropRule", "DROP_BASIC",
00570             "Type of dropping strategy to use",
00571             tuple<string>("DROP_BASIC","DROP_PROWS",
00572             "DROP_COLUMN","DROP_AREA",
00573             "DROP_DYNAMIC","DROP_INTERP"),
00574             tuple<string>("ILUTP(t)","ILUTP(p,t)",
00575             "Variant of ILUTP(p,t) for j-th column",
00576             "Variant of ILUTP to control memory",
00577             "Dynamically adjust threshold",
00578             "Compute second dropping threshold by interpolation"),
00579             tuple<SLU::rule_t>(SLU::DROP_BASIC,SLU::DROP_PROWS,SLU::DROP_COLUMN,
00580             SLU::DROP_AREA,SLU::DROP_DYNAMIC,SLU::DROP_INTERP),
00581             pl.getRawPtr());*/
00582 
00583     pl->set("ILU_DropTol", 0.0001, "ILUT drop tolerance");
00584 
00585     pl->set("ILU_FillFactor", 10.0, "ILUT fill factor");
00586 
00587     setStringToIntegralParameter<SLU::norm_t>("ILU_Norm", "INF_NORM",
00588             "Type of norm to use",
00589             tuple<string>("ONE_NORM","TWO_NORM","INF_NORM"),
00590             tuple<string>("1-norm","2-norm","inf-norm"),
00591             tuple<SLU::norm_t>(SLU::ONE_NORM,SLU::TWO_NORM,SLU::INF_NORM),
00592             pl.getRawPtr());
00593 
00594     setStringToIntegralParameter<SLU::milu_t>("ILU_MILU", "SILU",
00595             "Type of modified ILU to use",
00596             tuple<string>("SILU","SMILU_1","SMILU_2","SMILU_3"),
00597             tuple<string>("Regular ILU","MILU 1","MILU 2","MILU 3"),
00598             tuple<SLU::milu_t>(SLU::SILU,SLU::SMILU_1,SLU::SMILU_2,
00599             SLU::SMILU_3),
00600             pl.getRawPtr());
00601 
00602     pl->set("ILU_FillTol", 0.01, "ILUT fill tolerance");
00603 
00604     pl->set("ILU_Flag", false, "ILU flag: if true, run ILU routines");
00605 
00606     valid_params = pl;
00607   }
00608 
00609   return valid_params;
00610 }
00611 
00612 
00613 template <class Matrix, class Vector>
00614 bool
00615 Superlu<Matrix,Vector>::loadA_impl(EPhase current_phase)
00616 {
00617   using Teuchos::as;
00618 
00619 #ifdef HAVE_AMESOS2_TIMERS
00620   Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
00621 #endif
00622 
00623   // SuperLU does not need the matrix at symbolicFactorization()
00624   if( current_phase == SYMBFACT ) return false;
00625 
00626   // Cleanup old store memory if it's non-NULL (should only ever be non-NULL at root_)
00627   if( data_.A.Store != NULL ){
00628     SLU::Destroy_SuperMatrix_Store( &(data_.A) );
00629     data_.A.Store = NULL;
00630   }
00631 
00632   // Only the root image needs storage allocated
00633   if( this->root_ ){
00634     nzvals_.resize(this->globalNumNonZeros_);
00635     rowind_.resize(this->globalNumNonZeros_);
00636     colptr_.resize(this->globalNumCols_ + 1);
00637   }
00638 
00639   int nnz_ret = 0;
00640   {
00641 #ifdef HAVE_AMESOS2_TIMERS
00642     Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
00643 #endif
00644 
00645     TEUCHOS_TEST_FOR_EXCEPTION( this->rowIndexBase_ != this->columnIndexBase_,
00646                         std::runtime_error,
00647                         "Row and column maps have different indexbase ");
00648     Util::get_ccs_helper<
00649     MatrixAdapter<Matrix>,slu_type,int,int>::do_get(this->matrixA_.ptr(),
00650                                                     nzvals_(), rowind_(),
00651                                                     colptr_(), nnz_ret, ROOTED,
00652                                                     ARBITRARY,
00653                                                     this->rowIndexBase_);
00654   }
00655 
00656   // Get the SLU data type for this type of matrix
00657   SLU::Dtype_t dtype = type_map::dtype;
00658 
00659   if( this->root_ ){
00660     TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
00661                         std::runtime_error,
00662                         "Did not get the expected number of non-zero vals");
00663 
00664     function_map::create_CompCol_Matrix( &(data_.A),
00665                                          this->globalNumRows_, this->globalNumCols_,
00666                                          nnz_ret,
00667                                          nzvals_.getRawPtr(),
00668                                          rowind_.getRawPtr(),
00669                                          colptr_.getRawPtr(),
00670                                          SLU::SLU_NC, dtype, SLU::SLU_GE);
00671   }
00672 
00673   return true;
00674 }
00675 
00676 
00677 template<class Matrix, class Vector>
00678 const char* Superlu<Matrix,Vector>::name = "SuperLU";
00679   
00680 
00681 } // end namespace Amesos2
00682 
00683 #endif  // AMESOS2_SUPERLU_DEF_HPP