Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Superlumt_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Amesos2: Templated Direct Sparse Solver Package
6 // Copyright 2011 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //
42 // @HEADER
43 
52 #ifndef AMESOS2_SUPERLUMT_DEF_HPP
53 #define AMESOS2_SUPERLUMT_DEF_HPP
54 
55 #include <Teuchos_Tuple.hpp>
56 #include <Teuchos_ParameterList.hpp>
57 #include <Teuchos_StandardParameterEntryValidators.hpp>
58 
60 #include "Amesos2_Util.hpp"
61 
62 
63 namespace SLUMT {
64  /*
65  * We have to declare this extern function because of a bug in the
66  * SuperLU_MT header files. In each header is declared a function
67  * "extern int xsp_colorder(...)" where x is in {'s','d','c','z'},
68  * but this function is never defined anywhere, so if you use the
69  * function as defined, it will compile fine, but will choke during
70  * linking. No other code in SuperLU_MT actually calls these
71  * functions. Instead, all other SuperLU_MT function just call
72  * "sp_colorder", which is defined within the SuperLU_MT library,
73  * but not declared.
74  */
75  extern "C" {
76  int sp_colorder(SuperMatrix*,int*,superlumt_options_t*,SuperMatrix*);
77  }
78 } // end namespace SLUMT
79 
80 
81 namespace Amesos2 {
82 
83  /*
84  * Note: Although many of the type definitions for SuperLU_MT are
85  * identical to those of SuperLU, we do not mix the definitions so
86  * that we do not introduce messy coupling between the two
87  * interfaces. Also, there exist enough differences between the two
88  * packages to merit dedicated utility functions.
89  *
90  * We have also taken a different approach to interfacing with
91  * SuperLU_MT than with SuperLU which I believe leads to a better
92  * seperation of the 4 solution steps. We may in the future adopt a
93  * similar strategy for SuperLU.
94  */
95 
96  template <class Matrix, class Vector>
97  Superlumt<Matrix,Vector>::Superlumt(Teuchos::RCP<const Matrix> A,
98  Teuchos::RCP<Vector> X,
99  Teuchos::RCP<const Vector> B )
100  : SolverCore<Amesos2::Superlumt,Matrix,Vector>(A, X, B)
101  , nzvals_() // initialize to empty arrays
102  , rowind_()
103  , colptr_()
104  , is_contiguous_(true)
105  {
106  Teuchos::RCP<Teuchos::ParameterList> default_params
107  = Teuchos::parameterList( *(this->getValidParameters()) );
108  this->setParameters(default_params);
109 
110  data_.options.lwork = 0; // Use system memory for factorization
111 
112  data_.perm_r.resize(this->globalNumRows_);
113  data_.perm_c.resize(this->globalNumCols_);
114  data_.options.perm_r = data_.perm_r.getRawPtr();
115  data_.options.perm_c = data_.perm_c.getRawPtr();
116 
117  data_.R.resize(this->globalNumRows_);
118  data_.C.resize(this->globalNumCols_);
119 
120  data_.options.refact = SLUMT::NO; // initially we are not refactoring
121  data_.equed = SLUMT::NOEQUIL; // No equilibration has yet been performed
122  data_.rowequ = false;
123  data_.colequ = false;
124 
125  data_.A.Store = NULL;
126  data_.AC.Store = NULL;
127  data_.BX.Store = NULL;
128  data_.L.Store = NULL;
129  data_.U.Store = NULL;
130 
131  data_.stat.ops = NULL;
132  }
133 
134 
135  template <class Matrix, class Vector>
137  {
138  /* Free SuperLU_MT data_types
139  * - Matrices
140  * - Vectors
141  * - Stat object
142  */
143 
144  // Storage is initialized in numericFactorization_impl()
145  if ( data_.A.Store != NULL ){
146  // Our Teuchos::Array's will destroy rowind, colptr, and nzval for us
147  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
148  }
149 
150  // Cleanup memory allocated during last call to sp_colorder if needed
151  if( data_.AC.Store != NULL ){
152  SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
153  }
154 
155  if ( data_.L.Store != NULL ){ // will only ever be true for this->root_
156  SLUMT::D::Destroy_SuperNode_SCP( &(data_.L) );
157  SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
158 
159  // memory alloc'd in sp_colorder
160  free( data_.options.etree );
161  free( data_.options.colcnt_h );
162  free( data_.options.part_super_h );
163  }
164 
165 
166  // Storage is initialized in solve_impl()
167  if ( data_.BX.Store != NULL ){
168  /* Cannot use SLU::Destroy_Dense_Matrix routine here, since it attempts to
169  * free the array of non-zero values, but that array has already been
170  * deallocated by the MultiVector object. So we release just the Store
171  * instead.
172  */
173  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
174  }
175 
176  if ( data_.stat.ops != NULL )
177  SLUMT::D::StatFree( &(data_.stat) );
178  }
179 
180  template<class Matrix, class Vector>
181  int
183  {
184  // Use either the column-ordering found in the users perm_c or the requested computed ordering
185  int perm_spec = data_.options.ColPerm;
186  if( perm_spec != SLUMT::MY_PERMC && this->root_ ){
187 #ifdef HAVE_AMESOS2_TIMERS
188  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
189 #endif
190 
191  SLUMT::S::get_perm_c(perm_spec, &(data_.A), data_.perm_c.getRawPtr());
192  }
193  // Ordering will be applied directly before numeric factorization
194 
195  return(0);
196  }
197 
198 
199 
200  template <class Matrix, class Vector>
201  int
203  {
204  // We assume that a call to symbolicFactorization is indicating that
205  // the structure of the matrix has changed. SuperLU doesn't allow
206  // us to do a symbolic factorization directly, but we can leave a
207  // flag that tells it it needs to symbolically factor later.
208  data_.options.refact = SLUMT::NO;
209 
210  if( this->status_.numericFactorizationDone() && this->root_ ){
211  // If we've done a numeric factorization already, then we need to
212  // cleanup the old L and U. Stores and other data will be
213  // allocated during numeric factorization. Only rank 0 has valid
214  // pointers
215  SLUMT::D::Destroy_SuperNode_Matrix( &(data_.L) );
216  SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
217  data_.L.Store = NULL;
218  data_.U.Store = NULL;
219  }
220 
221  return(0);
222  }
223 
224 
225  template <class Matrix, class Vector>
226  int
228  {
229  using Teuchos::as;
230 
231 #ifdef HAVE_AMESOS2_DEBUG
232  const int nprocs = data_.options.nprocs;
233  TEUCHOS_TEST_FOR_EXCEPTION( nprocs <= 0,
234  std::invalid_argument,
235  "The number of threads to spawn should be greater than 0." );
236 #endif
237 
238  int info = 0;
239 
240  if ( this->root_ ) {
241 
242  if( data_.options.fact == SLUMT::EQUILIBRATE ){
243  magnitude_type rowcnd, colcnd, amax;
244  int info;
245 
246  function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
247  data_.C.getRawPtr(), &rowcnd, &colcnd,
248  &amax, &info);
249  TEUCHOS_TEST_FOR_EXCEPTION( info != 0,
250  std::runtime_error,
251  "SuperLU_MT gsequ returned with status " << info );
252 
253  function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
254  data_.C.getRawPtr(), rowcnd, colcnd,
255  amax, &(data_.equed));
256 
257  data_.rowequ = (data_.equed == SLUMT::ROW) || (data_.equed == SLUMT::BOTH);
258  data_.colequ = (data_.equed == SLUMT::COL) || (data_.equed == SLUMT::BOTH);
259 
260  data_.options.fact = SLUMT::DOFACT;
261  }
262 
263  // Cleanup memory allocated during last call to sp_colorder if needed
264  if( data_.AC.Store != NULL ){
265  SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
266  if( data_.options.refact == SLUMT::NO ){ // then we recompute etree; free the old one
267  free( data_.options.etree );
268  free( data_.options.colcnt_h );
269  free( data_.options.part_super_h );
270  }
271  data_.AC.Store = NULL;
272  }
273 
274  // Apply the column ordering, so that AC is the column-permuted A, and compute etree
275  SLUMT::sp_colorder(&(data_.A), data_.perm_c.getRawPtr(),
276  &(data_.options), &(data_.AC));
277 
278 
279  // Allocate and initialize status variable
280  const int n = as<int>(this->globalNumCols_); // n is the number of columns in A
281  if( data_.stat.ops != NULL ){ SLUMT::D::StatFree( &(data_.stat) ); data_.stat.ops = NULL; }
282  SLUMT::D::StatAlloc(n, data_.options.nprocs, data_.options.panel_size,
283  data_.options.relax, &(data_.stat));
284  SLUMT::D::StatInit(n, data_.options.nprocs, &(data_.stat));
285 
286 
287  { // Do factorization
288 #ifdef HAVE_AMESOS2_TIMERS
289  Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
290 #endif
291 
292 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
293  std::cout << "SuperLU_MT:: Before numeric factorization" << std::endl;
294  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
295  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
296  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
297 #endif
298 
299  function_map::gstrf(&(data_.options), &(data_.AC),
300  data_.perm_r.getRawPtr(), &(data_.L), &(data_.U),
301  &(data_.stat), &info);
302  }
303 
304  // Set the number of non-zero values in the L and U factors
305  this->setNnzLU( as<size_t>(((SLUMT::SCformat*)data_.L.Store)->nnz +
306  ((SLUMT::NCformat*)data_.U.Store)->nnz) );
307  }
308 
309  // Check output
310  const global_size_type info_st = as<global_size_type>(info);
311  TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
312  std::runtime_error,
313  "Factorization complete, but matrix is singular. Division by zero eminent");
314  TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
315  std::runtime_error,
316  "Memory allocation failure in SuperLU_MT factorization");
317  // The other option, that info_st < 0 denotes invalid parameters to
318  // the function, but we'll assume for now that that won't happen.
319 
320  data_.options.fact = SLUMT::FACTORED;
321  data_.options.refact = SLUMT::YES;
322 
323  /* All processes should return the same error code */
324  Teuchos::broadcast(*(this->getComm()),0,&info);
325  return(info);
326  }
327 
328 
329  template <class Matrix, class Vector>
330  int
332  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
333  {
334  using Teuchos::as;
335 
336  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
337  const size_t nrhs = X->getGlobalNumVectors();
338 
339  Teuchos::Array<slu_type> bxvals_(ld_rhs * nrhs);
340  size_t ldbx_ = as<size_t>(ld_rhs);
341 
342  { // Get values from B
343 #ifdef HAVE_AMESOS2_TIMERS
344  Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
345  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
346 #endif
347 
348  if ( is_contiguous_ == true ) {
350  slu_type>::do_get(B, bxvals_(),
351  ldbx_,
352  ROOTED, this->rowIndexBase_);
353  }
354  else {
356  slu_type>::do_get(B, bxvals_(),
357  ldbx_,
358  NON_CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
359  }
360  }
361 
362  int info = 0; // returned error code (0 = success)
363  if ( this->root_ ) {
364  // Clean up old B stores if it has already been created
365  if( data_.BX.Store != NULL ){
366  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
367  data_.BX.Store = NULL;
368  }
369 
370  {
371 #ifdef HAVE_AMESOS2_TIMERS
372  Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
373 #endif
374  SLUMT::Dtype_t dtype = type_map::dtype;
375  function_map::create_Dense_Matrix(&(data_.BX), as<int>(ld_rhs), as<int>(nrhs),
376  bxvals_.getRawPtr(), as<int>(ldbx_),
377  SLUMT::SLU_DN, dtype, SLUMT::SLU_GE);
378  }
379 
380  if( data_.options.trans == SLUMT::NOTRANS ){
381  if( data_.rowequ ){ // row equilibration has been done on AC
382  // scale bxvals_ by diag(R)
383  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
384  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
385  }
386  } else if( data_.colequ ){ // column equilibration has been done on AC
387  // scale bxvals_ by diag(C)
388  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
389  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
390  }
391 
392 
393  {
394 #ifdef HAVE_AMESOS2_TIMERS
395  Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
396 #endif
397 
398  function_map::gstrs(data_.options.trans, &(data_.L),
399  &(data_.U), data_.perm_r.getRawPtr(),
400  data_.perm_c.getRawPtr(), &(data_.BX),
401  &(data_.stat), &info);
402  }
403  } // end block for solve time
404 
405  /* All processes should have the same error code */
406  Teuchos::broadcast(*(this->getComm()),0,&info);
407 
408  TEUCHOS_TEST_FOR_EXCEPTION( info < 0,
409  std::runtime_error,
410  "Argument " << -info << " to gstrs had an illegal value" );
411 
412  // "Un-scale" the solution so that it is a solution of the original system
413  if( data_.options.trans == SLUMT::NOTRANS ){
414  if( data_.colequ ){ // column equilibration has been done on AC
415  // scale bxvals_ by diag(C)
416  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
417  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
418  }
419  } else if( data_.rowequ ){ // row equilibration has been done on AC
420  // scale bxvals_ by diag(R)
421  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
422  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
423  }
424 
425  /* Update X's global values */
426  {
427 #ifdef HAVE_AMESOS2_TIMERS
428  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
429 #endif
430 
431  if ( is_contiguous_ == true ) {
433  MultiVecAdapter<Vector>, slu_type>::do_put(X, bxvals_(), ldbx_, ROOTED);
434  }
435  else {
437  MultiVecAdapter<Vector>, slu_type>::do_put(X, bxvals_(), ldbx_, NON_CONTIGUOUS_AND_ROOTED);
438  }
439  }
440 
441  return(info);
442  }
443 
444 
445  template <class Matrix, class Vector>
446  bool
448  {
449  // The SuperLU_MT factorization routines can handle square as well as
450  // rectangular matrices, but SuperLU_MT can only apply the solve routines to
451  // square matrices, so we check the matrix for squareness.
452  return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
453  }
454 
455 
456  template <class Matrix, class Vector>
457  void
458  Superlumt<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
459  {
460  using Teuchos::as;
461  using Teuchos::RCP;
462  using Teuchos::getIntegralValue;
463  using Teuchos::ParameterEntryValidator;
464 
465  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
466 
467 
468  data_.options.nprocs = parameterList->get<int>("nprocs", 1);
469 
470  data_.options.trans = this->control_.useTranspose_ ? SLUMT::TRANS : SLUMT::NOTRANS;
471  // SuperLU_MT "trans" option can override the Amesos2 option
472  if( parameterList->isParameter("trans") ){
473  RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("trans").validator();
474  parameterList->getEntry("trans").setValidator(trans_validator);
475 
476  data_.options.trans = getIntegralValue<SLUMT::trans_t>(*parameterList, "trans");
477  }
478 
479  data_.options.panel_size = parameterList->get<int>("panel_size", SLUMT::D::sp_ienv(1));
480 
481  data_.options.relax = parameterList->get<int>("relax", SLUMT::D::sp_ienv(2));
482 
483  const bool equil = parameterList->get<bool>("Equil", true);
484  data_.options.fact = equil ? SLUMT::EQUILIBRATE : SLUMT::DOFACT;
485 
486  const bool symmetric_mode = parameterList->get<bool>("SymmetricMode", false);
487  data_.options.SymmetricMode = symmetric_mode ? SLUMT::YES : SLUMT::NO;
488 
489  const bool print_stat = parameterList->get<bool>("PrintStat", false);
490  data_.options.PrintStat = print_stat ? SLUMT::YES : SLUMT::NO;
491 
492  data_.options.diag_pivot_thresh = parameterList->get<double>("diag_pivot_thresh", 1.0);
493 
494  if( parameterList->isParameter("ColPerm") ){
495  RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
496  parameterList->getEntry("ColPerm").setValidator(colperm_validator);
497 
498  data_.options.ColPerm = getIntegralValue<SLUMT::colperm_t>(*parameterList, "ColPerm");
499  }
500 
501  // TODO: until we support retrieving col/row permutations from the user
502  data_.options.usepr = SLUMT::NO;
503 
504  if( parameterList->isParameter("IsContiguous") ){
505  is_contiguous_ = parameterList->get<bool>("IsContiguous");
506  }
507  }
508 
509 
510  template <class Matrix, class Vector>
511  Teuchos::RCP<const Teuchos::ParameterList>
513  {
514  using std::string;
515  using Teuchos::tuple;
516  using Teuchos::ParameterList;
517  using Teuchos::EnhancedNumberValidator;
518  using Teuchos::setStringToIntegralParameter;
519  using Teuchos::stringToIntegralParameterEntryValidator;
520 
521  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
522 
523  if( is_null(valid_params) ){
524  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
525 
526  Teuchos::RCP<EnhancedNumberValidator<int> > nprocs_validator
527  = Teuchos::rcp( new EnhancedNumberValidator<int>() );
528  nprocs_validator->setMin(1);
529  pl->set("nprocs", 1, "The number of processors to factorize with", nprocs_validator);
530 
531  setStringToIntegralParameter<SLUMT::trans_t>("trans", "NOTRANS",
532  "Solve for the transpose system or not",
533  tuple<string>("TRANS","NOTRANS","CONJ"),
534  tuple<string>("Solve with transpose",
535  "Do not solve with transpose",
536  "Solve with the conjugate transpose"),
537  tuple<SLUMT::trans_t>(SLUMT::TRANS,
538  SLUMT::NOTRANS,
539  SLUMT::CONJ),
540  pl.getRawPtr());
541 
542  Teuchos::RCP<EnhancedNumberValidator<int> > panel_size_validator
543  = Teuchos::rcp( new EnhancedNumberValidator<int>() );
544  panel_size_validator->setMin(1);
545  pl->set("panel_size", SLUMT::D::sp_ienv(1),
546  "Specifies the number of consecutive columns to be treated as a unit of task.",
547  panel_size_validator);
548 
549  Teuchos::RCP<EnhancedNumberValidator<int> > relax_validator
550  = Teuchos::rcp( new EnhancedNumberValidator<int>() );
551  relax_validator->setMin(1);
552  pl->set("relax", SLUMT::D::sp_ienv(2),
553  "Specifies the number of columns to be grouped as a relaxed supernode",
554  relax_validator);
555 
556  pl->set("Equil", true, "Whether to equilibrate the system before solve");
557 
558  Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
559  = Teuchos::rcp( new EnhancedNumberValidator<double>(0.0, 1.0) );
560  pl->set("diag_pivot_thresh", 1.0,
561  "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
562  diag_pivot_thresh_validator); // partial pivoting
563 
564  // Note: MY_PERMC not yet supported
565  setStringToIntegralParameter<SLUMT::colperm_t>("ColPerm", "COLAMD",
566  "Specifies how to permute the columns of the "
567  "matrix for sparsity preservation",
568  tuple<string>("NATURAL",
569  "MMD_AT_PLUS_A",
570  "MMD_ATA",
571  "COLAMD"),
572  tuple<string>("Natural ordering",
573  "Minimum degree ordering on A^T + A",
574  "Minimum degree ordering on A^T A",
575  "Approximate minimum degree column ordering"),
576  tuple<SLUMT::colperm_t>(SLUMT::NATURAL,
577  SLUMT::MMD_AT_PLUS_A,
578  SLUMT::MMD_ATA,
579  SLUMT::COLAMD),
580  pl.getRawPtr());
581 
582  pl->set("SymmetricMode", false, "Specifies whether to use the symmetric mode");
583 
584  // TODO: until we support getting row/col permutations from user
585  // pl->set("usepr", false);
586 
587  pl->set("PrintStat", false, "Specifies whether to print the solver's statistics");
588 
589  pl->set("IsContiguous", true, "Whether GIDs contiguous");
590 
591  valid_params = pl;
592  }
593 
594  return valid_params;
595  }
596 
597 
598  template <class Matrix, class Vector>
599  bool
601  {
602  using Teuchos::as;
603 
604 #ifdef HAVE_AMESOS2_TIMERS
605  Teuchos::TimeMonitor convTimer( this->timers_.mtxConvTime_ );
606 #endif
607 
608  if( current_phase == SYMBFACT ) return false;
609 
610  // Store is allocated on create_CompCol_Matrix
611  if( data_.A.Store != NULL ){
612  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
613  data_.A.Store = NULL;
614  }
615 
616  if( this->root_ ){
617  nzvals_.resize(this->globalNumNonZeros_);
618  rowind_.resize(this->globalNumNonZeros_);
619  colptr_.resize(this->globalNumCols_ + 1);
620  }
621 
622  int nnz_ret = 0;
623  {
624 #ifdef HAVE_AMESOS2_TIMERS
625  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
626 #endif
627 
628  // Use compressed-column store for SuperLU_MT
629  if ( is_contiguous_ == true ) {
631  MatrixAdapter<Matrix>,slu_type,int,int>::do_get(this->matrixA_.ptr(),
632  nzvals_, rowind_, colptr_,
633  nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_);
634  }
635  else {
637  MatrixAdapter<Matrix>,slu_type,int,int>::do_get(this->matrixA_.ptr(),
638  nzvals_, rowind_, colptr_,
639  nnz_ret, NON_CONTIGUOUS_AND_ROOTED, ARBITRARY, this->rowIndexBase_);
640  }
641  }
642 
643  if( this->root_ ){
644  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
645  std::runtime_error,
646  "rank=0 failed to get all nonzero values in getCcs()");
647 
648  SLUMT::Dtype_t dtype = type_map::dtype;
649  function_map::create_CompCol_Matrix(&(data_.A),
650  as<int>(this->globalNumRows_),
651  as<int>(this->globalNumCols_),
652  nnz_ret,
653  nzvals_.getRawPtr(),
654  rowind_.getRawPtr(),
655  colptr_.getRawPtr(),
656  SLUMT::SLU_NC,
657  dtype, SLUMT::SLU_GE);
658 
659  TEUCHOS_TEST_FOR_EXCEPTION( data_.A.Store == NULL,
660  std::runtime_error,
661  "Failed to allocate matrix store" );
662  }
663 
664  return true;
665  }
666 
667 
668  template<class Matrix, class Vector>
669  const char* Superlumt<Matrix,Vector>::name = "SuperLU_MT";
670 
671 
672 } // end namespace Amesos2
673 
674 #endif // AMESOS2_SUPERLUMT_DEF_HPP
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Superlumt_def.hpp:447
Amesos2 interface to the Multi-threaded version of SuperLU.
Definition: Amesos2_Superlumt_decl.hpp:77
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
SuperLU_MT specific solve.
Definition: Amesos2_Superlumt_def.hpp:331
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:478
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:475
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using SuperLU_MT.
Definition: Amesos2_Superlumt_def.hpp:202
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:266
Utility functions for Amesos2.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Superlumt_def.hpp:512
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a const parameter list of all of the valid parameters that this-&gt;setParameterList(...) will accept.
Definition: Amesos2_SolverCore_def.hpp:314
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:654
Superlumt(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Superlumt_def.hpp:97
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Superlumt_def.hpp:182
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
int numericFactorization_impl()
SuperLU_MT specific numeric factorization.
Definition: Amesos2_Superlumt_def.hpp:227
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_Superlumt_def.hpp:458
~Superlumt()
Destructor.
Definition: Amesos2_Superlumt_def.hpp:136
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:282
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:322
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Superlumt_def.hpp:600