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   SLU::set_default_options(&(data_.options));
00077   // Override some default options
00078   data_.options.PrintStat = SLU::NO;
00079 
00080   SLU::StatInit(&(data_.stat));
00081 
00082   data_.perm_r.resize(this->globalNumRows_);
00083   data_.perm_c.resize(this->globalNumCols_);
00084   data_.etree.resize(this->globalNumCols_);
00085   data_.R.resize(this->globalNumRows_);
00086   data_.C.resize(this->globalNumCols_);
00087 
00088   data_.relax = SLU::sp_ienv(2); // Query optimal relax param from superlu
00089   data_.panel_size = SLU::sp_ienv(1); // Query optimal panel size
00090 
00091   data_.equed = 'N';            // No equilibration
00092   data_.A.Store = NULL;
00093   data_.L.Store = NULL;
00094   data_.U.Store = NULL;
00095   data_.X.Store = NULL;
00096   data_.B.Store = NULL;
00097 }
00098 
00099 
00100 template <class Matrix, class Vector>
00101 Superlu<Matrix,Vector>::~Superlu( )
00102 {
00103   /* Free Superlu data_types
00104    * - Matrices
00105    * - Vectors
00106    * - Stat object
00107    */
00108   SLU::StatFree( &(data_.stat) ) ;
00109 
00110   // Storage is initialized in numericFactorization_impl()
00111   if ( data_.A.Store != NULL ){
00112     SLU::Destroy_SuperMatrix_Store( &(data_.A) );
00113   }
00114 
00115   // only root allocated these SuperMatrices.
00116   if ( data_.L.Store != NULL ){ // will only be true for this->root_
00117     SLU::Destroy_SuperNode_Matrix( &(data_.L) );
00118     SLU::Destroy_CompCol_Matrix( &(data_.U) );
00119   }
00120 }
00121 
00122 template<class Matrix, class Vector>
00123 int
00124 Superlu<Matrix,Vector>::preOrdering_impl()
00125 {
00126   /*
00127    * Get column permutation vector perm_c[], according to permc_spec:
00128    *   permc_spec = NATURAL:  natural ordering
00129    *   permc_spec = MMD_AT_PLUS_A: minimum degree on structure of A'+A
00130    *   permc_spec = MMD_ATA:  minimum degree on structure of A'*A
00131    *   permc_spec = COLAMD:   approximate minimum degree column ordering
00132    *   permc_spec = MY_PERMC: the ordering already supplied in perm_c[]
00133    */
00134   int permc_spec = data_.options.ColPerm;
00135   if ( permc_spec != SLU::MY_PERMC && this->root_ ){
00136 #ifdef HAVE_AMESOS2_TIMERS
00137     Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
00138 #endif
00139 
00140     SLU::get_perm_c(permc_spec, &(data_.A), data_.perm_c.getRawPtr());
00141   }
00142 
00143   return(0);
00144 }
00145 
00146 
00147 template <class Matrix, class Vector>
00148 int
00149 Superlu<Matrix,Vector>::symbolicFactorization_impl()
00150 {
00151   /*
00152    * SuperLU performs symbolic factorization and numeric factorization
00153    * together, but does leave some options for reusing symbolic
00154    * structures that have been created on previous factorizations.  If
00155    * our Amesos2 user calls this function, that is an indication that
00156    * the symbolic structure of the matrix is no longer valid, and
00157    * SuperLU should do the factorization from scratch.
00158    *
00159    * This can be accomplished by setting the options.Fact flag to
00160    * DOFACT, as well as setting our own internal flag to false.
00161    */
00162   same_symbolic_ = false;
00163   data_.options.Fact = SLU::DOFACT;
00164 
00165   return(0);
00166 }
00167 
00168 
00169 template <class Matrix, class Vector>
00170 int
00171 Superlu<Matrix,Vector>::numericFactorization_impl()
00172 {
00173   using Teuchos::as;
00174 
00175   // Cleanup old L and U matrices if we are not reusing a symbolic
00176   // factorization.  Stores and other data will be allocated in gstrf.
00177   // Only rank 0 has valid pointers
00178   if ( !same_symbolic_ && data_.L.Store != NULL ){
00179     SLU::Destroy_SuperNode_Matrix( &(data_.L) );
00180     SLU::Destroy_CompCol_Matrix( &(data_.U) );
00181     data_.L.Store = NULL;
00182     data_.U.Store = NULL;
00183   }
00184 
00185   if( same_symbolic_ ) data_.options.Fact = SLU::SamePattern_SameRowPerm;
00186 
00187   int info = 0;
00188   if ( this->root_ ){
00189 
00190 #ifdef HAVE_AMESOS2_DEBUG
00191     TEUCHOS_TEST_FOR_EXCEPTION( data_.A.ncol != as<int>(this->globalNumCols_),
00192                         std::runtime_error,
00193                         "Error in converting to SuperLU SuperMatrix: wrong number of global columns." );
00194     TEUCHOS_TEST_FOR_EXCEPTION( data_.A.nrow != as<int>(this->globalNumRows_),
00195                         std::runtime_error,
00196                         "Error in converting to SuperLU SuperMatrix: wrong number of global rows." );
00197 #endif
00198 
00199     if( data_.options.Equil == SLU::YES ){
00200       magnitude_type rowcnd, colcnd, amax;
00201       int info2 = 0;
00202 
00203       // calculate row and column scalings
00204       function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
00205                           data_.C.getRawPtr(), &rowcnd, &colcnd,
00206                           &amax, &info2);
00207       TEUCHOS_TEST_FOR_EXCEPTION( info2 != 0,
00208                           std::runtime_error,
00209                           "SuperLU gsequ returned with status " << info2 );
00210 
00211       // apply row and column scalings if necessary
00212       function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
00213                           data_.C.getRawPtr(), rowcnd, colcnd,
00214                           amax, &(data_.equed));
00215 
00216       // // check what types of equilibration was actually done
00217       // data_.rowequ = (data_.equed == 'R') || (data_.equed == 'B');
00218       // data_.colequ = (data_.equed == 'C') || (data_.equed == 'B');
00219     }
00220 
00221     // Apply the column permutation computed in preOrdering.  Place the
00222     // column-permuted matrix in AC
00223     SLU::sp_preorder(&(data_.options), &(data_.A), data_.perm_c.getRawPtr(),
00224                      data_.etree.getRawPtr(), &(data_.AC));
00225 
00226     { // Do factorization
00227 #ifdef HAVE_AMESOS2_TIMERS
00228       Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
00229 #endif
00230 
00231 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
00232       std::cout << "Superlu:: Before numeric factorization" << std::endl;
00233       std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
00234       std::cout << "rowind_ : " << rowind_.toString() << std::endl;
00235       std::cout << "colptr_ : " << colptr_.toString() << std::endl;
00236 #endif
00237 
00238       function_map::gstrf(&(data_.options), &(data_.AC),
00239                           data_.relax, data_.panel_size, data_.etree.getRawPtr(),
00240                           NULL, 0, data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
00241                           &(data_.L), &(data_.U), &(data_.stat), &info);
00242     }
00243     // Cleanup. AC data will be alloc'd again for next factorization (if at all)
00244     SLU::Destroy_CompCol_Permuted( &(data_.AC) );
00245 
00246     // Set the number of non-zero values in the L and U factors
00247     this->setNnzLU( as<size_t>(((SLU::SCformat*)data_.L.Store)->nnz +
00248                                ((SLU::NCformat*)data_.U.Store)->nnz) );
00249   }
00250 
00251   /* All processes should have the same error code */
00252   Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
00253 
00254   global_size_type info_st = as<global_size_type>(info);
00255   TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
00256     std::runtime_error,
00257     "Factorization complete, but matrix is singular. Division by zero eminent");
00258   TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
00259     std::runtime_error,
00260     "Memory allocation failure in Superlu factorization");
00261 
00262   data_.options.Fact = SLU::FACTORED;
00263   same_symbolic_ = true;
00264 
00265   return(info);
00266 }
00267 
00268 
00269 template <class Matrix, class Vector>
00270 int
00271 Superlu<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> >       X,
00272                                    const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
00273 {
00274   using Teuchos::as;
00275 
00276   const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
00277   const size_t nrhs = X->getGlobalNumVectors();
00278 
00279   const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
00280   Teuchos::Array<slu_type> xValues(val_store_size);
00281   Teuchos::Array<slu_type> bValues(val_store_size);
00282 
00283   {                             // Get values from RHS B
00284 #ifdef HAVE_AMESOS2_TIMERS
00285     Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
00286     Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
00287 #endif
00288     Util::get_1d_copy_helper<MultiVecAdapter<Vector>,
00289                              slu_type>::do_get(B, bValues(),
00290                                                as<size_t>(ld_rhs),
00291                                                ROOTED);
00292   }
00293 
00294   int ierr = 0; // returned error code
00295 
00296   magnitude_type rpg, rcond;
00297   if ( this->root_ ) {
00298     data_.ferr.resize(nrhs);
00299     data_.berr.resize(nrhs);
00300 
00301     {
00302 #ifdef HAVE_AMESOS2_TIMERS
00303       Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
00304 #endif
00305       SLU::Dtype_t dtype = type_map::dtype;
00306       int i_ld_rhs = as<int>(ld_rhs);
00307       function_map::create_Dense_Matrix(&(data_.B), i_ld_rhs, as<int>(nrhs),
00308                                         bValues.getRawPtr(), i_ld_rhs,
00309                                         SLU::SLU_DN, dtype, SLU::SLU_GE);
00310       function_map::create_Dense_Matrix(&(data_.X), i_ld_rhs, as<int>(nrhs),
00311                                         xValues.getRawPtr(), i_ld_rhs,
00312                                         SLU::SLU_DN, dtype, SLU::SLU_GE);
00313     }
00314 
00315     // Note: the values of B and X (after solution) are adjusted
00316     // appropriately within gssvx for row and column scalings.
00317 
00318     {                           // Do solve!
00319 #ifdef HAVE_AMESOS2_TIMERS
00320     Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
00321 #endif
00322 
00323     function_map::gssvx(&(data_.options), &(data_.A),
00324       data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(), data_.etree.getRawPtr(),
00325       &(data_.equed), data_.R.getRawPtr(), data_.C.getRawPtr(), &(data_.L),
00326       &(data_.U), NULL, 0, &(data_.B), &(data_.X), &rpg, &rcond,
00327       data_.ferr.getRawPtr(), data_.berr.getRawPtr(), &(data_.mem_usage),
00328       &(data_.stat), &ierr);
00329     }
00330 
00331     // Cleanup X and B stores
00332     SLU::Destroy_SuperMatrix_Store( &(data_.X) );
00333     SLU::Destroy_SuperMatrix_Store( &(data_.B) );
00334     data_.X.Store = NULL;
00335     data_.B.Store = NULL;
00336   }
00337 
00338   /* All processes should have the same error code */
00339   Teuchos::broadcast(*(this->getComm()), 0, &ierr);
00340 
00341   global_size_type ierr_st = as<global_size_type>(ierr);
00342   TEUCHOS_TEST_FOR_EXCEPTION( ierr < 0,
00343                       std::invalid_argument,
00344                       "Argument " << -ierr << " to SuperLU xgssvx had illegal value" );
00345   TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st <= this->globalNumCols_,
00346                       std::runtime_error,
00347                       "Factorization complete, but U is exactly singular" );
00348   TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st > this->globalNumCols_ + 1,
00349                       std::runtime_error,
00350                       "SuperLU allocated " << ierr - this->globalNumCols_ << " bytes of "
00351                       "memory before allocation failure occured." );
00352 
00353   /* Update X's global values */
00354   {
00355 #ifdef HAVE_AMESOS2_TIMERS
00356     Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
00357 #endif
00358 
00359     Util::put_1d_data_helper<
00360       MultiVecAdapter<Vector>,slu_type>::do_put(X, xValues(),
00361                                          as<size_t>(ld_rhs),
00362                                          ROOTED);
00363   }
00364 
00365 
00366   return(ierr);
00367 }
00368 
00369 
00370 template <class Matrix, class Vector>
00371 bool
00372 Superlu<Matrix,Vector>::matrixShapeOK_impl() const
00373 {
00374   // The Superlu factorization routines can handle square as well as
00375   // rectangular matrices, but Superlu can only apply the solve routines to
00376   // square matrices, so we check the matrix for squareness.
00377   return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
00378 }
00379 
00380 
00381 template <class Matrix, class Vector>
00382 void
00383 Superlu<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
00384 {
00385   using Teuchos::RCP;
00386   using Teuchos::getIntegralValue;
00387   using Teuchos::ParameterEntryValidator;
00388 
00389   RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
00390 
00391   data_.options.Trans = this->control_.useTranspose_ ? SLU::TRANS : SLU::NOTRANS;
00392   // The SuperLU transpose option can override the Amesos2 option
00393   if( parameterList->isParameter("Trans") ){
00394     RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("Trans").validator();
00395     parameterList->getEntry("Trans").setValidator(trans_validator);
00396 
00397     data_.options.Trans = getIntegralValue<SLU::trans_t>(*parameterList, "Trans");
00398   }
00399 
00400   if( parameterList->isParameter("IterRefine") ){
00401     RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry("IterRefine").validator();
00402     parameterList->getEntry("IterRefine").setValidator(refine_validator);
00403 
00404     data_.options.IterRefine = getIntegralValue<SLU::IterRefine_t>(*parameterList, "IterRefine");
00405   }
00406 
00407   if( parameterList->isParameter("ColPerm") ){
00408     RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
00409     parameterList->getEntry("ColPerm").setValidator(colperm_validator);
00410 
00411     data_.options.ColPerm = getIntegralValue<SLU::colperm_t>(*parameterList, "ColPerm");
00412   }
00413 
00414   data_.options.DiagPivotThresh = parameterList->get<double>("DiagPivotThresh", 1.0);
00415 
00416   bool equil = parameterList->get<bool>("Equil", true);
00417   data_.options.Equil = equil ? SLU::YES : SLU::NO;
00418 
00419   bool symmetric_mode = parameterList->get<bool>("SymmetricMode", false);
00420   data_.options.SymmetricMode = symmetric_mode ? SLU::YES : SLU::NO;
00421 }
00422 
00423 
00424 template <class Matrix, class Vector>
00425 Teuchos::RCP<const Teuchos::ParameterList>
00426 Superlu<Matrix,Vector>::getValidParameters_impl() const
00427 {
00428   using std::string;
00429   using Teuchos::tuple;
00430   using Teuchos::ParameterList;
00431   using Teuchos::EnhancedNumberValidator;
00432   using Teuchos::setStringToIntegralParameter;
00433   using Teuchos::stringToIntegralParameterEntryValidator;
00434 
00435   static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
00436 
00437   if( is_null(valid_params) ){
00438     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00439 
00440     setStringToIntegralParameter<SLU::trans_t>("Trans", "NOTRANS",
00441                                                "Solve for the transpose system or not",
00442                                                tuple<string>("TRANS","NOTRANS","CONJ"),
00443                                                tuple<string>("Solve with transpose",
00444                                                              "Do not solve with transpose",
00445                                                              "Solve with the conjugate transpose"),
00446                                                tuple<SLU::trans_t>(SLU::TRANS,
00447                                                                    SLU::NOTRANS,
00448                                                                    SLU::CONJ),
00449                                                pl.getRawPtr());
00450 
00451     setStringToIntegralParameter<SLU::IterRefine_t>("IterRefine", "NOREFINE",
00452                                                     "Type of iterative refinement to use",
00453                                                     tuple<string>("NOREFINE", "SINGLE", "DOUBLE"),
00454                                                     tuple<string>("Do not use iterative refinement",
00455                                                                   "Do single iterative refinement",
00456                                                                   "Do double iterative refinement"),
00457                                                     tuple<SLU::IterRefine_t>(SLU::NOREFINE,
00458                                                                              SLU::SINGLE,
00459                                                                              SLU::DOUBLE),
00460                                                     pl.getRawPtr());
00461 
00462     // Note: MY_PERMC not yet supported
00463     setStringToIntegralParameter<SLU::colperm_t>("ColPerm", "COLAMD",
00464                                                  "Specifies how to permute the columns of the "
00465                                                  "matrix for sparsity preservation",
00466                                                  tuple<string>("NATURAL", "MMD_AT_PLUS_A",
00467                                                                "MMD_ATA", "COLAMD"),
00468                                                  tuple<string>("Natural ordering",
00469                                                                "Minimum degree ordering on A^T + A",
00470                                                                "Minimum degree ordering on A^T A",
00471                                                                "Approximate minimum degree column ordering"),
00472                                                  tuple<SLU::colperm_t>(SLU::NATURAL,
00473                                                                        SLU::MMD_AT_PLUS_A,
00474                                                                        SLU::MMD_ATA,
00475                                                                        SLU::COLAMD),
00476                                                  pl.getRawPtr());
00477 
00478     Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
00479       = Teuchos::rcp( new EnhancedNumberValidator<double>(0.0, 1.0) );
00480     pl->set("DiagPivotThresh", 1.0,
00481             "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
00482             diag_pivot_thresh_validator); // partial pivoting
00483 
00484     pl->set("Equil", true, "Whether to equilibrate the system before solve");
00485 
00486     pl->set("SymmetricMode", false,
00487             "Specifies whether to use the symmetric mode. "
00488             "Gives preference to diagonal pivots and uses "
00489             "an (A^T + A)-based column permutation.");
00490 
00491     valid_params = pl;
00492   }
00493 
00494   return valid_params;
00495 }
00496 
00497 
00498 template <class Matrix, class Vector>
00499 bool
00500 Superlu<Matrix,Vector>::loadA_impl(EPhase current_phase)
00501 {
00502   using Teuchos::as;
00503 
00504 #ifdef HAVE_AMESOS2_TIMERS
00505   Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
00506 #endif
00507 
00508   // SuperLU does not need the matrix at symbolicFactorization()
00509   if( current_phase == SYMBFACT ) return false;
00510 
00511   // Cleanup old store memory if it's non-NULL (should only ever be non-NULL at root_)
00512   if( data_.A.Store != NULL ){
00513     SLU::Destroy_SuperMatrix_Store( &(data_.A) );
00514     data_.A.Store = NULL;
00515   }
00516 
00517   // Only the root image needs storage allocated
00518   if( this->root_ ){
00519     nzvals_.resize(this->globalNumNonZeros_);
00520     rowind_.resize(this->globalNumNonZeros_);
00521     colptr_.resize(this->globalNumCols_ + 1);
00522   }
00523 
00524   int nnz_ret = 0;
00525   {
00526 #ifdef HAVE_AMESOS2_TIMERS
00527     Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
00528 #endif
00529 
00530     Util::get_ccs_helper<
00531     MatrixAdapter<Matrix>,slu_type,int,int>::do_get(this->matrixA_.ptr(),
00532                                                     nzvals_(), rowind_(), colptr_(),
00533                                                     nnz_ret, ROOTED, ARBITRARY);
00534   }
00535 
00536   // Get the SLU data type for this type of matrix
00537   SLU::Dtype_t dtype = type_map::dtype;
00538 
00539   if( this->root_ ){
00540     TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
00541                         std::runtime_error,
00542                         "Did not get the expected number of non-zero vals");
00543 
00544     function_map::create_CompCol_Matrix( &(data_.A),
00545                                          this->globalNumRows_, this->globalNumCols_,
00546                                          nnz_ret,
00547                                          nzvals_.getRawPtr(),
00548                                          rowind_.getRawPtr(),
00549                                          colptr_.getRawPtr(),
00550                                          SLU::SLU_NC, dtype, SLU::SLU_GE);
00551   }
00552 
00553   return true;
00554 }
00555 
00556 
00557 template<class Matrix, class Vector>
00558 const char* Superlu<Matrix,Vector>::name = "SuperLU";
00559   
00560 
00561 } // end namespace Amesos2
00562 
00563 #endif  // AMESOS2_SUPERLU_DEF_HPP