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  {
105  Teuchos::RCP<Teuchos::ParameterList> default_params
106  = Teuchos::parameterList( *(this->getValidParameters()) );
107  this->setParameters(default_params);
108 
109  data_.options.lwork = 0; // Use system memory for factorization
110 
111  data_.perm_r.resize(this->globalNumRows_);
112  data_.perm_c.resize(this->globalNumCols_);
113  data_.options.perm_r = data_.perm_r.getRawPtr();
114  data_.options.perm_c = data_.perm_c.getRawPtr();
115 
116  data_.R.resize(this->globalNumRows_);
117  data_.C.resize(this->globalNumCols_);
118 
119  data_.options.refact = SLUMT::NO; // initially we are not refactoring
120  data_.equed = SLUMT::NOEQUIL; // No equilibration has yet been performed
121  data_.rowequ = false;
122  data_.colequ = false;
123 
124  data_.A.Store = NULL;
125  data_.AC.Store = NULL;
126  data_.BX.Store = NULL;
127  data_.L.Store = NULL;
128  data_.U.Store = NULL;
129 
130  data_.stat.ops = NULL;
131  }
132 
133 
134  template <class Matrix, class Vector>
136  {
137  /* Free SuperLU_MT data_types
138  * - Matrices
139  * - Vectors
140  * - Stat object
141  */
142 
143  // Storage is initialized in numericFactorization_impl()
144  if ( data_.A.Store != NULL ){
145  // Our Teuchos::Array's will destroy rowind, colptr, and nzval for us
146  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
147  }
148 
149  // Cleanup memory allocated during last call to sp_colorder if needed
150  if( data_.AC.Store != NULL ){
151  SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
152  }
153 
154  if ( data_.L.Store != NULL ){ // will only ever be true for this->root_
155  SLUMT::D::Destroy_SuperNode_SCP( &(data_.L) );
156  SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
157 
158  // memory alloc'd in sp_colorder
159  free( data_.options.etree );
160  free( data_.options.colcnt_h );
161  free( data_.options.part_super_h );
162  }
163 
164 
165  // Storage is initialized in solve_impl()
166  if ( data_.BX.Store != NULL ){
167  /* Cannot use SLU::Destroy_Dense_Matrix routine here, since it attempts to
168  * free the array of non-zero values, but that array has already been
169  * deallocated by the MultiVector object. So we release just the Store
170  * instead.
171  */
172  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
173  }
174 
175  if ( data_.stat.ops != NULL )
176  SLUMT::D::StatFree( &(data_.stat) );
177  }
178 
179  template<class Matrix, class Vector>
180  int
182  {
183  // Use either the column-ordering found in the users perm_c or the requested computed ordering
184  int perm_spec = data_.options.ColPerm;
185  if( perm_spec != SLUMT::MY_PERMC && this->root_ ){
186 #ifdef HAVE_AMESOS2_TIMERS
187  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
188 #endif
189 
190  SLUMT::S::get_perm_c(perm_spec, &(data_.A), data_.perm_c.getRawPtr());
191  }
192  // Ordering will be applied directly before numeric factorization
193 
194  return(0);
195  }
196 
197 
198 
199  template <class Matrix, class Vector>
200  int
202  {
203  // We assume that a call to symbolicFactorization is indicating that
204  // the structure of the matrix has changed. SuperLU doesn't allow
205  // us to do a symbolic factorization directly, but we can leave a
206  // flag that tells it it needs to symbolically factor later.
207  data_.options.refact = SLUMT::NO;
208 
209  if( this->status_.numericFactorizationDone() && this->root_ ){
210  // If we've done a numeric factorization already, then we need to
211  // cleanup the old L and U. Stores and other data will be
212  // allocated during numeric factorization. Only rank 0 has valid
213  // pointers
214  SLUMT::D::Destroy_SuperNode_Matrix( &(data_.L) );
215  SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
216  data_.L.Store = NULL;
217  data_.U.Store = NULL;
218  }
219 
220  return(0);
221  }
222 
223 
224  template <class Matrix, class Vector>
225  int
227  {
228  using Teuchos::as;
229 
230 #ifdef HAVE_AMESOS2_DEBUG
231  const int nprocs = data_.options.nprocs;
232  TEUCHOS_TEST_FOR_EXCEPTION( nprocs <= 0,
233  std::invalid_argument,
234  "The number of threads to spawn should be greater than 0." );
235 #endif
236 
237  int info = 0;
238 
239  if ( this->root_ ) {
240 
241  if( data_.options.fact == SLUMT::EQUILIBRATE ){
242  magnitude_type rowcnd, colcnd, amax;
243  int info;
244 
245  function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
246  data_.C.getRawPtr(), &rowcnd, &colcnd,
247  &amax, &info);
248  TEUCHOS_TEST_FOR_EXCEPTION( info != 0,
249  std::runtime_error,
250  "SuperLU_MT gsequ returned with status " << info );
251 
252  function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
253  data_.C.getRawPtr(), rowcnd, colcnd,
254  amax, &(data_.equed));
255 
256  data_.rowequ = (data_.equed == SLUMT::ROW) || (data_.equed == SLUMT::BOTH);
257  data_.colequ = (data_.equed == SLUMT::COL) || (data_.equed == SLUMT::BOTH);
258 
259  data_.options.fact = SLUMT::DOFACT;
260  }
261 
262  // Cleanup memory allocated during last call to sp_colorder if needed
263  if( data_.AC.Store != NULL ){
264  SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) ); // free's colbeg, colend, and Store
265  if( data_.options.refact == SLUMT::NO ){ // then we recompute etree; free the old one
266  free( data_.options.etree );
267  free( data_.options.colcnt_h );
268  free( data_.options.part_super_h );
269  }
270  data_.AC.Store = NULL;
271  }
272 
273  // Apply the column ordering, so that AC is the column-permuted A, and compute etree
274  SLUMT::sp_colorder(&(data_.A), data_.perm_c.getRawPtr(),
275  &(data_.options), &(data_.AC));
276 
277 
278  // Allocate and initialize status variable
279  const int n = as<int>(this->globalNumCols_); // n is the number of columns in A
280  if( data_.stat.ops != NULL ){ SLUMT::D::StatFree( &(data_.stat) ); data_.stat.ops = NULL; }
281  SLUMT::D::StatAlloc(n, data_.options.nprocs, data_.options.panel_size,
282  data_.options.relax, &(data_.stat));
283  SLUMT::D::StatInit(n, data_.options.nprocs, &(data_.stat));
284 
285 
286  { // Do factorization
287 #ifdef HAVE_AMESOS2_TIMERS
288  Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
289 #endif
290 
291 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
292  std::cout << "SuperLU_MT:: Before numeric factorization" << std::endl;
293  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
294  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
295  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
296 #endif
297 
298  function_map::gstrf(&(data_.options), &(data_.AC),
299  data_.perm_r.getRawPtr(), &(data_.L), &(data_.U),
300  &(data_.stat), &info);
301  }
302 
303  // Set the number of non-zero values in the L and U factors
304  this->setNnzLU( as<size_t>(((SLUMT::SCformat*)data_.L.Store)->nnz +
305  ((SLUMT::NCformat*)data_.U.Store)->nnz) );
306  }
307 
308  // Check output
309  const global_size_type info_st = as<global_size_type>(info);
310  TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
311  std::runtime_error,
312  "Factorization complete, but matrix is singular. Division by zero eminent");
313  TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
314  std::runtime_error,
315  "Memory allocation failure in SuperLU_MT factorization");
316  // The other option, that info_st < 0 denotes invalid parameters to
317  // the function, but we'll assume for now that that won't happen.
318 
319  data_.options.fact = SLUMT::FACTORED;
320  data_.options.refact = SLUMT::YES;
321 
322  /* All processes should return the same error code */
323  Teuchos::broadcast(*(this->getComm()),0,&info);
324  return(info);
325  }
326 
327 
328  template <class Matrix, class Vector>
329  int
331  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
332  {
333  using Teuchos::as;
334 
335  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
336  const size_t nrhs = X->getGlobalNumVectors();
337 
338  Teuchos::Array<slu_type> bxvals_(ld_rhs * nrhs);
339  size_t ldbx_ = as<size_t>(ld_rhs);
340 
341  { // Get values from B
342 #ifdef HAVE_AMESOS2_TIMERS
343  Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
344  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
345 #endif
346 
348  slu_type>::do_get(B, bxvals_(),
349  ldbx_,
350  ROOTED);
351  }
352 
353  int info = 0; // returned error code (0 = success)
354  if ( this->root_ ) {
355  // Clean up old B stores if it has already been created
356  if( data_.BX.Store != NULL ){
357  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
358  data_.BX.Store = NULL;
359  }
360 
361  {
362 #ifdef HAVE_AMESOS2_TIMERS
363  Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
364 #endif
365  SLUMT::Dtype_t dtype = type_map::dtype;
366  function_map::create_Dense_Matrix(&(data_.BX), as<int>(ld_rhs), as<int>(nrhs),
367  bxvals_.getRawPtr(), as<int>(ldbx_),
368  SLUMT::SLU_DN, dtype, SLUMT::SLU_GE);
369  }
370 
371  if( data_.options.trans == SLUMT::NOTRANS ){
372  if( data_.rowequ ){ // row equilibration has been done on AC
373  // scale bxvals_ by diag(R)
374  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
375  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
376  }
377  } else if( data_.colequ ){ // column equilibration has been done on AC
378  // scale bxvals_ by diag(C)
379  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
380  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
381  }
382 
383 
384  {
385 #ifdef HAVE_AMESOS2_TIMERS
386  Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
387 #endif
388 
389  function_map::gstrs(data_.options.trans, &(data_.L),
390  &(data_.U), data_.perm_r.getRawPtr(),
391  data_.perm_c.getRawPtr(), &(data_.BX),
392  &(data_.stat), &info);
393  }
394  } // end block for solve time
395 
396  /* All processes should have the same error code */
397  Teuchos::broadcast(*(this->getComm()),0,&info);
398 
399  TEUCHOS_TEST_FOR_EXCEPTION( info < 0,
400  std::runtime_error,
401  "Argument " << -info << " to gstrs had an illegal value" );
402 
403  // "Un-scale" the solution so that it is a solution of the original system
404  if( data_.options.trans == SLUMT::NOTRANS ){
405  if( data_.colequ ){ // column equilibration has been done on AC
406  // scale bxvals_ by diag(C)
407  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
408  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
409  }
410  } else if( data_.rowequ ){ // row equilibration has been done on AC
411  // scale bxvals_ by diag(R)
412  Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
413  SLUMT::slu_mt_mult<slu_type,magnitude_type>());
414  }
415 
416  /* Update X's global values */
417  {
418 #ifdef HAVE_AMESOS2_TIMERS
419  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
420 #endif
421 
423  MultiVecAdapter<Vector>, slu_type>::do_put(X, bxvals_(), ldbx_, ROOTED);
424  }
425 
426  return(info);
427  }
428 
429 
430  template <class Matrix, class Vector>
431  bool
433  {
434  // The SuperLU_MT factorization routines can handle square as well as
435  // rectangular matrices, but SuperLU_MT can only apply the solve routines to
436  // square matrices, so we check the matrix for squareness.
437  return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
438  }
439 
440 
441  template <class Matrix, class Vector>
442  void
443  Superlumt<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
444  {
445  using Teuchos::as;
446  using Teuchos::RCP;
447  using Teuchos::getIntegralValue;
448  using Teuchos::ParameterEntryValidator;
449 
450  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
451 
452 
453  data_.options.nprocs = parameterList->get<int>("nprocs", 1);
454 
455  data_.options.trans = this->control_.useTranspose_ ? SLUMT::TRANS : SLUMT::NOTRANS;
456  // SuperLU_MT "trans" option can override the Amesos2 option
457  if( parameterList->isParameter("trans") ){
458  RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("trans").validator();
459  parameterList->getEntry("trans").setValidator(trans_validator);
460 
461  data_.options.trans = getIntegralValue<SLUMT::trans_t>(*parameterList, "trans");
462  }
463 
464  data_.options.panel_size = parameterList->get<int>("panel_size", SLUMT::D::sp_ienv(1));
465 
466  data_.options.relax = parameterList->get<int>("relax", SLUMT::D::sp_ienv(2));
467 
468  const bool equil = parameterList->get<bool>("Equil", true);
469  data_.options.fact = equil ? SLUMT::EQUILIBRATE : SLUMT::DOFACT;
470 
471  const bool symmetric_mode = parameterList->get<bool>("SymmetricMode", false);
472  data_.options.SymmetricMode = symmetric_mode ? SLUMT::YES : SLUMT::NO;
473 
474  const bool print_stat = parameterList->get<bool>("PrintStat", false);
475  data_.options.PrintStat = print_stat ? SLUMT::YES : SLUMT::NO;
476 
477  data_.options.diag_pivot_thresh = parameterList->get<double>("diag_pivot_thresh", 1.0);
478 
479  if( parameterList->isParameter("ColPerm") ){
480  RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
481  parameterList->getEntry("ColPerm").setValidator(colperm_validator);
482 
483  data_.options.ColPerm = getIntegralValue<SLUMT::colperm_t>(*parameterList, "ColPerm");
484  }
485 
486  // TODO: until we support retrieving col/row permutations from the user
487  data_.options.usepr = SLUMT::NO;
488  }
489 
490 
491  template <class Matrix, class Vector>
492  Teuchos::RCP<const Teuchos::ParameterList>
494  {
495  using std::string;
496  using Teuchos::tuple;
497  using Teuchos::ParameterList;
498  using Teuchos::EnhancedNumberValidator;
499  using Teuchos::setStringToIntegralParameter;
500  using Teuchos::stringToIntegralParameterEntryValidator;
501 
502  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
503 
504  if( is_null(valid_params) ){
505  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
506 
507  Teuchos::RCP<EnhancedNumberValidator<int> > nprocs_validator
508  = Teuchos::rcp( new EnhancedNumberValidator<int>() );
509  nprocs_validator->setMin(1);
510  pl->set("nprocs", 1, "The number of processors to factorize with", nprocs_validator);
511 
512  setStringToIntegralParameter<SLUMT::trans_t>("trans", "NOTRANS",
513  "Solve for the transpose system or not",
514  tuple<string>("TRANS","NOTRANS","CONJ"),
515  tuple<string>("Solve with transpose",
516  "Do not solve with transpose",
517  "Solve with the conjugate transpose"),
518  tuple<SLUMT::trans_t>(SLUMT::TRANS,
519  SLUMT::NOTRANS,
520  SLUMT::CONJ),
521  pl.getRawPtr());
522 
523  Teuchos::RCP<EnhancedNumberValidator<int> > panel_size_validator
524  = Teuchos::rcp( new EnhancedNumberValidator<int>() );
525  panel_size_validator->setMin(1);
526  pl->set("panel_size", SLUMT::D::sp_ienv(1),
527  "Specifies the number of consecutive columns to be treated as a unit of task.",
528  panel_size_validator);
529 
530  Teuchos::RCP<EnhancedNumberValidator<int> > relax_validator
531  = Teuchos::rcp( new EnhancedNumberValidator<int>() );
532  relax_validator->setMin(1);
533  pl->set("relax", SLUMT::D::sp_ienv(2),
534  "Specifies the number of columns to be grouped as a relaxed supernode",
535  relax_validator);
536 
537  pl->set("Equil", true, "Whether to equilibrate the system before solve");
538 
539  Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
540  = Teuchos::rcp( new EnhancedNumberValidator<double>(0.0, 1.0) );
541  pl->set("diag_pivot_thresh", 1.0,
542  "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
543  diag_pivot_thresh_validator); // partial pivoting
544 
545  // Note: MY_PERMC not yet supported
546  setStringToIntegralParameter<SLUMT::colperm_t>("ColPerm", "COLAMD",
547  "Specifies how to permute the columns of the "
548  "matrix for sparsity preservation",
549  tuple<string>("NATURAL",
550  "MMD_AT_PLUS_A",
551  "MMD_ATA",
552  "COLAMD"),
553  tuple<string>("Natural ordering",
554  "Minimum degree ordering on A^T + A",
555  "Minimum degree ordering on A^T A",
556  "Approximate minimum degree column ordering"),
557  tuple<SLUMT::colperm_t>(SLUMT::NATURAL,
558  SLUMT::MMD_AT_PLUS_A,
559  SLUMT::MMD_ATA,
560  SLUMT::COLAMD),
561  pl.getRawPtr());
562 
563  pl->set("SymmetricMode", false, "Specifies whether to use the symmetric mode");
564 
565  // TODO: until we support getting row/col permutations from user
566  // pl->set("usepr", false);
567 
568  pl->set("PrintStat", false, "Specifies whether to print the solver's statistics");
569 
570  valid_params = pl;
571  }
572 
573  return valid_params;
574  }
575 
576 
577  template <class Matrix, class Vector>
578  bool
580  {
581  using Teuchos::as;
582 
583 #ifdef HAVE_AMESOS2_TIMERS
584  Teuchos::TimeMonitor convTimer( this->timers_.mtxConvTime_ );
585 #endif
586 
587  if( current_phase == SYMBFACT ) return false;
588 
589  // Store is allocated on create_CompCol_Matrix
590  if( data_.A.Store != NULL ){
591  SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
592  data_.A.Store = NULL;
593  }
594 
595  if( this->root_ ){
596  nzvals_.resize(this->globalNumNonZeros_);
597  rowind_.resize(this->globalNumNonZeros_);
598  colptr_.resize(this->globalNumCols_ + 1);
599  }
600 
601  int nnz_ret = 0;
602  {
603 #ifdef HAVE_AMESOS2_TIMERS
604  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
605 #endif
606 
607  // Use compressed-column store for SuperLU_MT
609  MatrixAdapter<Matrix>,slu_type,int,int>::do_get(this->matrixA_.ptr(),
611  nnz_ret, ROOTED, ARBITRARY);
612  }
613 
614  if( this->root_ ){
615  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
616  std::runtime_error,
617  "rank=0 failed to get all nonzero values in getCcs()");
618 
619  SLUMT::Dtype_t dtype = type_map::dtype;
620  function_map::create_CompCol_Matrix(&(data_.A),
621  as<int>(this->globalNumRows_),
622  as<int>(this->globalNumCols_),
623  nnz_ret,
624  nzvals_.getRawPtr(),
625  rowind_.getRawPtr(),
626  colptr_.getRawPtr(),
627  SLUMT::SLU_NC,
628  dtype, SLUMT::SLU_GE);
629 
630  TEUCHOS_TEST_FOR_EXCEPTION( data_.A.Store == NULL,
631  std::runtime_error,
632  "Failed to allocate matrix store" );
633  }
634 
635  return true;
636  }
637 
638 
639  template<class Matrix, class Vector>
640  const char* Superlumt<Matrix,Vector>::name = "SuperLU_MT";
641 
642 
643 } // end namespace Amesos2
644 
645 #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:432
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:330
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
Teuchos::Array< int > colptr_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_Superlumt_decl.hpp:274
bool root_
If true, then this is the root processor.
Definition: Amesos2_SolverCore_decl.hpp:507
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using SuperLU_MT.
Definition: Amesos2_Superlumt_def.hpp:201
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
Definition: Amesos2_TypeDecl.hpp:142
Utility functions for Amesos2.
Control control_
Parameters for solving.
Definition: Amesos2_SolverCore_decl.hpp:495
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Superlumt_def.hpp:493
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a const parameter list of all of the valid parameters that this->setParameterList(...) will accept.
Definition: Amesos2_SolverCore_def.hpp:307
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
Definition: Amesos2_Superlumt_def.hpp:63
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Superlumt_def.hpp:181
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns a pointer to the Teuchos::Comm communicator with this operator.
Definition: Amesos2_SolverCore_decl.hpp:363
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:226
bool numericFactorizationDone() const
If true , then numeric factorization has been performed.
Definition: Amesos2_Status.hpp:118
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_Superlumt_def.hpp:443
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
Teuchos::Array< typename TypeMap< Amesos2::Superlumt, scalar_type >::type > nzvals_
Stores the values of the nonzero entries for SuperLU.
Definition: Amesos2_Superlumt_decl.hpp:270
~Superlumt()
Destructor.
Definition: Amesos2_Superlumt_def.hpp:135
void setNnzLU(size_t nnz)
Set the number of non-zero values in the and factors.
Definition: Amesos2_SolverCore_decl.hpp:452
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:275
Definition: Amesos2_TypeDecl.hpp:127
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:296
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:175
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Superlumt_def.hpp:579
Status status_
Holds status information about a solver.
Definition: Amesos2_SolverCore_decl.hpp:492
Teuchos::Array< int > rowind_
Stores the row indices of the nonzero entries.
Definition: Amesos2_Superlumt_decl.hpp:272