Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Superludist_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_SUPERLUDIST_DEF_HPP
53 #define AMESOS2_SUPERLUDIST_DEF_HPP
54 
55 #include <Teuchos_Tuple.hpp>
56 #include <Teuchos_StandardParameterEntryValidators.hpp>
57 #include <Teuchos_DefaultMpiComm.hpp>
58 
61 #include "Amesos2_Util.hpp"
62 
63 
64 namespace Amesos2 {
65 
66 
67  template <class Matrix, class Vector>
68  Superludist<Matrix,Vector>::Superludist(Teuchos::RCP<const Matrix> A,
69  Teuchos::RCP<Vector> X,
70  Teuchos::RCP<const Vector> B)
71  : SolverCore<Amesos2::Superludist,Matrix,Vector>(A, X, B)
72  , nzvals_() // initialization to empty arrays
73  , colind_()
74  , rowptr_()
75  , bvals_()
76  , xvals_()
77  , in_grid_(false)
78  {
79  using Teuchos::Comm;
80  // It's OK to depend on MpiComm explicitly here, because
81  // SuperLU_DIST requires MPI anyway.
82  using Teuchos::MpiComm;
83  using Teuchos::outArg;
84  using Teuchos::ParameterList;
85  using Teuchos::parameterList;
86  using Teuchos::RCP;
87  using Teuchos::rcp;
88  using Teuchos::rcp_dynamic_cast;
89  using Teuchos::REDUCE_SUM;
90  using Teuchos::reduceAll;
91  typedef global_ordinal_type GO;
92  typedef Tpetra::Map<local_ordinal_type, GO, node_type> map_type;
93 
95  // Set up the SuperLU_DIST processor grid //
97 
98  RCP<const Comm<int> > comm = this->getComm ();
99  const int myRank = comm->getRank ();
100  const int numProcs = comm->getSize ();
101 
102  SLUD::int_t nprow, npcol;
103  get_default_grid_size (numProcs, nprow, npcol);
104 
105  {
106  // FIXME (mfh 16 Dec 2014) getComm() just returns
107  // matrixA_->getComm(), so it's not clear why we need to ask for
108  // the matrix's communicator separately here.
109  RCP<const Comm<int> > matComm = this->matrixA_->getComm ();
110  TEUCHOS_TEST_FOR_EXCEPTION(
111  matComm.is_null (), std::logic_error, "Amesos2::Superlustdist "
112  "constructor: The matrix's communicator is null!");
113  RCP<const MpiComm<int> > matMpiComm =
114  rcp_dynamic_cast<const MpiComm<int> > (matComm);
115  // FIXME (mfh 16 Dec 2014) If the matrix's communicator is a
116  // SerialComm, we probably could just use MPI_COMM_SELF here.
117  // I'm not sure if SuperLU_DIST is smart enough to handle that
118  // case, though.
119  TEUCHOS_TEST_FOR_EXCEPTION(
120  matMpiComm.is_null (), std::logic_error, "Amesos2::Superlustdist "
121  "constructor: The matrix's communicator is not an MpiComm!");
122  TEUCHOS_TEST_FOR_EXCEPTION(
123  matMpiComm->getRawMpiComm ().is_null (), std::logic_error, "Amesos2::"
124  "Superlustdist constructor: The matrix's communicator claims to be a "
125  "Teuchos::MpiComm<int>, but its getRawPtrComm() method returns "
126  "Teuchos::null! This means that the underlying MPI_Comm doesn't even "
127  "exist, which likely implies that the Teuchos::MpiComm was constructed "
128  "incorrectly. It means something different than if the MPI_Comm were "
129  "MPI_COMM_NULL.");
130  MPI_Comm rawMpiComm = (* (matMpiComm->getRawMpiComm ())) ();
131  data_.mat_comm = rawMpiComm;
132  // This looks a bit like ScaLAPACK's grid initialization (which
133  // technically takes place in the BLACS, not in ScaLAPACK
134  // proper). See http://netlib.org/scalapack/slug/node34.html.
135  // The main difference is that SuperLU_DIST depends explicitly
136  // on MPI, while the BLACS hides its communication protocol.
137  SLUD::superlu_gridinit(data_.mat_comm, nprow, npcol, &(data_.grid));
138  }
139 
141  // Set some default parameters. //
142  // //
143  // Must do this after grid has been created in //
144  // case user specifies the nprow and npcol parameters //
146  RCP<ParameterList> default_params =
147  parameterList (* (this->getValidParameters ()));
148  this->setParameters (default_params);
149 
150  // Set some internal options
151  data_.options.Fact = SLUD::DOFACT;
152  data_.equed = SLUD::NOEQUIL; // No equilibration has yet been performed
153  data_.options.SolveInitialized = SLUD::NO;
154  data_.options.RefineInitialized = SLUD::NO;
155  data_.rowequ = false;
156  data_.colequ = false;
157  data_.perm_r.resize(this->globalNumRows_);
158  data_.perm_c.resize(this->globalNumCols_);
159 
161  // Set up a communicator for the parallel column ordering and //
162  // parallel symbolic factorization. //
164  data_.symb_comm = MPI_COMM_NULL;
165 
166  // domains is the next power of 2 less than nprow*npcol. This
167  // value will be used for creating an MPI communicator for the
168  // pre-ordering and symbolic factorization methods.
169  data_.domains = (int) ( pow(2.0, floor(log10((double)nprow*npcol)/log10(2.0))) );
170 
171  const int color = (myRank < data_.domains) ? 0 : MPI_UNDEFINED;
172  MPI_Comm_split (data_.mat_comm, color, myRank, &(data_.symb_comm));
173 
175  // Set up a row Map that only includes processes that are in the
176  // SuperLU process grid. This will be used for redistributing A.
178 
179  // mfh 16 Dec 2014: We could use createWeightedContigMapWithNode
180  // with myProcParticipates as the weight, but that costs an extra
181  // all-reduce.
182 
183  // Set to 1 if I am in the grid, and I get some of the matrix rows.
184  int myProcParticipates = 0;
185  if (myRank < nprow * npcol) {
186  in_grid_ = true;
187  myProcParticipates = 1;
188  }
189 
190  // Compute how many processes in the communicator belong to the
191  // process grid.
192  int numParticipatingProcs = 0;
193  reduceAll<int, int> (*comm, REDUCE_SUM, myProcParticipates,
194  outArg (numParticipatingProcs));
195  TEUCHOS_TEST_FOR_EXCEPTION(
196  this->globalNumRows_ != 0 && numParticipatingProcs == 0,
197  std::logic_error, "Amesos2::Superludist constructor: The matrix has "
198  << this->globalNumRows_ << " > 0 global row(s), but no processes in the "
199  "communicator want to participate in its factorization! nprow = "
200  << nprow << " and npcol = " << npcol << ".");
201 
202  // Divide up the rows among the participating processes.
203  size_t myNumRows = 0;
204  {
205  const GO GNR = static_cast<GO> (this->globalNumRows_);
206  const GO quotient = (numParticipatingProcs == 0) ? static_cast<GO> (0) :
207  GNR / static_cast<GO> (numParticipatingProcs);
208  const GO remainder =
209  GNR - quotient * static_cast<GO> (numParticipatingProcs);
210  const GO lclNumRows = (static_cast<GO> (myRank) < remainder) ?
211  (quotient + static_cast<GO> (1)) : quotient;
212  myNumRows = static_cast<size_t> (lclNumRows);
213  }
214 
215  // TODO: might only need to initialize if parallel symbolic factorization is requested.
216  const GO indexBase = this->matrixA_->getRowMap ()->getIndexBase ();
218  rcp (new map_type (this->globalNumRows_, myNumRows, indexBase, comm));
219 
221  // Do some other initialization //
223 
224  data_.A.Store = NULL;
225  function_map::LUstructInit(this->globalNumRows_, this->globalNumCols_, &(data_.lu));
226  SLUD::PStatInit(&(data_.stat));
227  // We do not use ScalePermstructInit because we will use our own
228  // arrays for storing perm_r and perm_c
229  data_.scale_perm.perm_r = data_.perm_r.getRawPtr();
230  data_.scale_perm.perm_c = data_.perm_c.getRawPtr();
231  }
232 
233 
234  template <class Matrix, class Vector>
236  {
237  /* Free SuperLU_DIST data_types
238  * - Matrices
239  * - Vectors
240  * - Stat object
241  * - ScalePerm, LUstruct, grid, and solve objects
242  *
243  * Note: the function definitions are the same regardless whether
244  * complex or real, so we arbitrarily use the D namespace
245  */
246  if ( this->status_.getNumPreOrder() > 0 ){
247  free( data_.sizes );
248  free( data_.fstVtxSep );
249  }
250 
251  // Cleanup old matrix store memory if it's non-NULL. Our
252  // Teuchos::Array's will destroy rowind, colptr, and nzval for us
253  if( data_.A.Store != NULL ){
254  SLUD::Destroy_SuperMatrix_Store_dist( &(data_.A) );
255  }
256 
257  // LU data is initialized in numericFactorization_impl()
258  if ( this->status_.getNumNumericFact() > 0 ){
259  function_map::Destroy_LU(this->globalNumRows_, &(data_.grid), &(data_.lu));
260  }
261  function_map::LUstructFree(&(data_.lu));
262 
263  // If a symbolic factorization is ever performed without a
264  // follow-up numericfactorization, there are some arrays in the
265  // Pslu_freeable struct which will never be free'd by
266  // SuperLU_DIST.
267  if ( this->status_.symbolicFactorizationDone() &&
269  if ( data_.pslu_freeable.xlsub != NULL ){
270  free( data_.pslu_freeable.xlsub );
271  free( data_.pslu_freeable.lsub );
272  }
273  if ( data_.pslu_freeable.xusub != NULL ){
274  free( data_.pslu_freeable.xusub );
275  free( data_.pslu_freeable.usub );
276  }
277  if ( data_.pslu_freeable.supno_loc != NULL ){
278  free( data_.pslu_freeable.supno_loc );
279  free( data_.pslu_freeable.xsup_beg_loc );
280  free( data_.pslu_freeable.xsup_end_loc );
281  }
282  free( data_.pslu_freeable.globToLoc );
283  }
284 
285  SLUD::PStatFree( &(data_.stat) ) ;
286 
287  // Teuchos::Arrays will free R, C, perm_r, and perm_c
288  // SLUD::D::ScalePermstructFree(&(data_.scale_perm));
289 
290  if ( data_.options.SolveInitialized == SLUD::YES )
291  function_map::SolveFinalize(&(data_.options), &(data_.solve_struct));
292 
293  SLUD::superlu_gridexit(&(data_.grid)); // TODO: are there any
294  // cases where grid
295  // wouldn't be initialized?
296 
297  if ( data_.symb_comm != MPI_COMM_NULL ) MPI_Comm_free(&(data_.symb_comm));
298  }
299 
300  template<class Matrix, class Vector>
301  int
303  {
304  // We will always use the NATURAL row ordering to avoid the
305  // sequential bottleneck present when doing any other row
306  // ordering scheme from SuperLU_DIST
307  //
308  // Set perm_r to be the natural ordering
309  SLUD::int_t slu_rows_ub = Teuchos::as<SLUD::int_t>(this->globalNumRows_);
310  for( SLUD::int_t i = 0; i < slu_rows_ub; ++i ) data_.perm_r[i] = i;
311 
312  // loadA_impl(); // Refresh matrix values
313 
314  if( in_grid_ ){
315  // If this function has been called at least once, then the
316  // sizes, and fstVtxSep arrays were allocated in
317  // get_perm_c_parmetis. Delete them before calling that
318  // function again. These arrays will also be dealloc'd in the
319  // deconstructor.
320  if( this->status_.getNumPreOrder() > 0 ){
321  free( data_.sizes );
322  free( data_.fstVtxSep );
323  }
324 #ifdef HAVE_AMESOS2_TIMERS
325  Teuchos::TimeMonitor preOrderTime( this->timers_.preOrderTime_ );
326 #endif
327 
328  float info = 0.0;
329  info = SLUD::get_perm_c_parmetis( &(data_.A),
330  data_.perm_r.getRawPtr(), data_.perm_c.getRawPtr(),
331  data_.grid.nprow * data_.grid.npcol, data_.domains,
332  &(data_.sizes), &(data_.fstVtxSep),
333  &(data_.grid), &(data_.symb_comm) );
334 
335  TEUCHOS_TEST_FOR_EXCEPTION( info > 0.0,
336  std::runtime_error,
337  "SuperLU_DIST pre-ordering ran out of memory after allocating "
338  << info << " bytes of memory" );
339  }
340 
341  // Ordering will be applied directly before numeric factorization,
342  // after we have a chance to get updated coefficients from the
343  // matrix
344 
345  return EXIT_SUCCESS;
346  }
347 
348 
349 
350  template <class Matrix, class Vector>
351  int
353  {
354  // loadA_impl(); // Refresh matrix values
355 
356  if( in_grid_ ){
357 
358 #ifdef HAVE_AMESOS2_TIMERS
359  Teuchos::TimeMonitor symFactTime( this->timers_.symFactTime_ );
360 #endif
361 
362  float info = 0.0;
363  info = SLUD::symbfact_dist((data_.grid.nprow) * (data_.grid.npcol),
364  data_.domains, &(data_.A), data_.perm_c.getRawPtr(),
365  data_.perm_r.getRawPtr(), data_.sizes,
366  data_.fstVtxSep, &(data_.pslu_freeable),
367  &(data_.grid.comm), &(data_.symb_comm),
368  &(data_.mem_usage));
369 
370  TEUCHOS_TEST_FOR_EXCEPTION( info > 0.0,
371  std::runtime_error,
372  "SuperLU_DIST symbolic factorization ran out of memory after"
373  " allocating " << info << " bytes of memory" );
374  }
375  same_symbolic_ = false;
376  same_solve_struct_ = false;
377 
378  return EXIT_SUCCESS;
379  }
380 
381 
382  template <class Matrix, class Vector>
383  int
385  using Teuchos::as;
386 
387  // loadA_impl(); // Refresh the matrix values
388 
389  // if( data_.options.Equil == SLUD::YES ){
390  // // Apply the scalings computed in preOrdering
391  // function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
392  // data_.C.getRawPtr(), data_.rowcnd, data_.colcnd,
393  // data_.amax, &(data_.equed));
394 
395  // data_.rowequ = (data_.equed == SLUD::ROW) || (data_.equed == SLUD::BOTH);
396  // data_.colequ = (data_.equed == SLUD::COL) || (data_.equed == SLUD::BOTH);
397  // }
398 
399  if( in_grid_ ){
400  // Apply the column ordering, so that AC is the column-permuted A, and compute etree
401  size_t nnz_loc = ((SLUD::NRformat_loc*)data_.A.Store)->nnz_loc;
402  for( size_t i = 0; i < nnz_loc; ++i ) colind_[i] = data_.perm_c[colind_[i]];
403 
404  // Distribute data from the symbolic factorization
405  if( same_symbolic_ ){
406  // Note: with the SamePattern_SameRowPerm options, it does not
407  // matter that the glu_freeable member has never been
408  // initialized, because it is never accessed. It is a
409  // placeholder arg. The real work is done in data_.lu
410  function_map::pdistribute(SLUD::SamePattern_SameRowPerm,
411  as<SLUD::int_t>(this->globalNumRows_), // aka "n"
412  &(data_.A), &(data_.scale_perm),
413  &(data_.glu_freeable), &(data_.lu),
414  &(data_.grid));
415  } else {
416  function_map::dist_psymbtonum(SLUD::DOFACT,
417  as<SLUD::int_t>(this->globalNumRows_), // aka "n"
418  &(data_.A), &(data_.scale_perm),
419  &(data_.pslu_freeable), &(data_.lu),
420  &(data_.grid));
421  }
422 
423  // Retrieve the normI of A (required by gstrf).
424  double anorm = function_map::plangs((char *)"I", &(data_.A), &(data_.grid));
425 
426  int info = 0;
427  {
428 #ifdef HAVE_AMESOS2_TIMERS
429  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
430 #endif
431 
432  function_map::gstrf(&(data_.options), this->globalNumRows_,
433  this->globalNumCols_, anorm, &(data_.lu),
434  &(data_.grid), &(data_.stat), &info);
435  }
436 
437  // Check output
438  TEUCHOS_TEST_FOR_EXCEPTION( info > 0,
439  std::runtime_error,
440  "L and U factors have been computed but U("
441  << info << "," << info << ") is exactly zero "
442  "(i.e. U is singular)");
443  }
444 
445  // The other option, that info_st < 0, denotes invalid parameters
446  // to the function, but we'll assume for now that that won't
447  // happen.
448 
449  data_.options.Fact = SLUD::FACTORED;
450  same_symbolic_ = true;
451 
452  return EXIT_SUCCESS;
453  }
454 
455 
456  template <class Matrix, class Vector>
457  int
459  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
460  {
461  using Teuchos::as;
462 
463  // local_len_rhs is how many of the multivector rows belong to
464  // this processor in the SuperLU_DIST processor grid.
465  const size_t local_len_rhs = superlu_rowmap_->getNodeNumElements();
466  const global_size_type nrhs = X->getGlobalNumVectors();
467  const global_ordinal_type first_global_row_b = superlu_rowmap_->getMinGlobalIndex();
468 
469  // make sure our multivector storage is sized appropriately
470  bvals_.resize(nrhs * local_len_rhs);
471  xvals_.resize(nrhs * local_len_rhs);
472 
473  // We assume the global length of the two vectors have already been
474  // checked for compatibility
475 
476  { // get the values from B
477 #ifdef HAVE_AMESOS2_TIMERS
478  Teuchos::TimeMonitor convTimer(this->timers_.vecConvTime_);
479 #endif
480 
481  {
482  // The input dense matrix for B should be distributed in the
483  // same manner as the superlu_dist matrix. That is, if a
484  // processor has m_loc rows of A, then it should also have
485  // m_loc rows of B (and the same rows). We accomplish this by
486  // distributing the multivector rows with the same Map that
487  // the matrix A's rows are distributed.
488 #ifdef HAVE_AMESOS2_TIMERS
489  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
490 #endif
491 
492  // get grid-distributed mv data. The multivector data will be
493  // distributed across the processes in the SuperLU_DIST grid.
494  typedef Util::get_1d_copy_helper<MultiVecAdapter<Vector>,slu_type> copy_helper;
495  copy_helper::do_get(B,
496  bvals_(),
497  local_len_rhs,
498  Teuchos::ptrInArg(*superlu_rowmap_));
499  }
500  } // end block for conversion time
501 
502  if( in_grid_ ){
503  // if( data_.options.trans == SLUD::NOTRANS ){
504  // if( data_.rowequ ){ // row equilibration has been done on AC
505  // // scale bxvals_ by diag(R)
506  // Util::scale(bxvals_(), as<size_t>(len_rhs), ldbx_, data_.R(),
507  // SLUD::slu_mt_mult<slu_type,magnitude_type>());
508  // }
509  // } else if( data_.colequ ){ // column equilibration has been done on AC
510  // // scale bxvals_ by diag(C)
511  // Util::scale(bxvals_(), as<size_t>(len_rhs), ldbx_, data_.C(),
512  // SLUD::slu_mt_mult<slu_type,magnitude_type>());
513  // }
514 
515  // Initialize the SOLVEstruct_t.
516  //
517  // We are able to reuse the solve struct if we have not changed
518  // the sparsity pattern of L and U since the last solve
519  if( !same_solve_struct_ ){
520  if( data_.options.SolveInitialized == SLUD::YES ){
521  function_map::SolveFinalize(&(data_.options), &(data_.solve_struct));
522  }
523  function_map::SolveInit(&(data_.options), &(data_.A), data_.perm_r.getRawPtr(),
524  data_.perm_c.getRawPtr(), as<SLUD::int_t>(nrhs), &(data_.lu),
525  &(data_.grid), &(data_.solve_struct));
526  // Flag that we can reuse this solve_struct unless another
527  // symbolicFactorization is called between here and the next
528  // solve.
529  same_solve_struct_ = true;
530  }
531 
532  int ierr = 0; // returned error code
533  {
534 #ifdef HAVE_AMESOS2_TIMERS
535  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
536 #endif
537 
538  function_map::gstrs(as<SLUD::int_t>(this->globalNumRows_), &(data_.lu),
539  &(data_.scale_perm), &(data_.grid), bvals_.getRawPtr(),
540  as<SLUD::int_t>(local_len_rhs), as<SLUD::int_t>(first_global_row_b),
541  as<SLUD::int_t>(local_len_rhs), as<int>(nrhs),
542  &(data_.solve_struct), &(data_.stat), &ierr);
543  } // end block for solve time
544 
545  TEUCHOS_TEST_FOR_EXCEPTION( ierr < 0,
546  std::runtime_error,
547  "Argument " << -ierr << " to gstrs had an illegal value" );
548 
549  // "Un-scale" the solution so that it is a solution of the original system
550  // if( data_.options.trans == SLUD::NOTRANS ){
551  // if( data_.colequ ){ // column equilibration has been done on AC
552  // // scale bxvals_ by diag(C)
553  // Util::scale(bxvals_(), as<size_t>(len_rhs), ldbx_, data_.C(),
554  // SLUD::slu_mt_mult<slu_type,magnitude_type>());
555  // }
556  // } else if( data_.rowequ ){ // row equilibration has been done on AC
557  // // scale bxvals_ by diag(R)
558  // Util::scale(bxvals_(), as<size_t>(len_rhs), ldbx_, data_.R(),
559  // SLUD::slu_mt_mult<slu_type,magnitude_type>());
560  // }
561  { // permute B to a solution of the original system
562 #ifdef HAVE_AMESOS2_TIMERS
563  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
564 #endif
565  SLUD::int_t ld = as<SLUD::int_t>(local_len_rhs);
566  function_map::permute_Dense_Matrix(as<SLUD::int_t>(first_global_row_b),
567  as<SLUD::int_t>(local_len_rhs),
568  data_.solve_struct.row_to_proc,
569  data_.solve_struct.inv_perm_c,
570  bvals_.getRawPtr(), ld,
571  xvals_.getRawPtr(), ld,
572  as<int>(nrhs),
573  &(data_.grid));
574  }
575  }
576 
577  /* Update X's global values */
578  {
579 #ifdef HAVE_AMESOS2_TIMERS
580  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
581 #endif
582 
583  typedef Util::put_1d_data_helper<MultiVecAdapter<Vector>,slu_type> put_helper;
584  put_helper::do_put(X,
585  xvals_(),
586  local_len_rhs,
587  Teuchos::ptrInArg(*superlu_rowmap_));
588  }
589 
590  return EXIT_SUCCESS;
591  }
592 
593 
594  template <class Matrix, class Vector>
595  bool
597  {
598  // SuperLU_DIST requires square matrices
599  return( this->globalNumRows_ == this->globalNumCols_ );
600  }
601 
602 
603  template <class Matrix, class Vector>
604  void
605  Superludist<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
606  {
607  using Teuchos::as;
608  using Teuchos::RCP;
609  using Teuchos::getIntegralValue;
610  using Teuchos::ParameterEntryValidator;
611 
612  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
613 
614  if( parameterList->isParameter("npcol") || parameterList->isParameter("nprow") ){
615  TEUCHOS_TEST_FOR_EXCEPTION( !(parameterList->isParameter("nprow") &&
616  parameterList->isParameter("npcol")),
617  std::invalid_argument,
618  "nprow and npcol must be set together" );
619 
620  SLUD::int_t nprow = parameterList->template get<SLUD::int_t>("nprow");
621  SLUD::int_t npcol = parameterList->template get<SLUD::int_t>("npcol");
622 
623  TEUCHOS_TEST_FOR_EXCEPTION( nprow * npcol > this->getComm()->getSize(),
624  std::invalid_argument,
625  "nprow and npcol combination invalid" );
626 
627  if( (npcol != data_.grid.npcol) || (nprow != data_.grid.nprow) ){
628  // De-allocate the default grid that was initialized in the constructor
629  SLUD::superlu_gridexit(&(data_.grid));
630  // Create a new grid
631  SLUD::superlu_gridinit(data_.mat_comm, nprow, npcol, &(data_.grid));
632  } // else our grid has not changed size since the last initialization
633  }
634 
635  TEUCHOS_TEST_FOR_EXCEPTION( this->control_.useTranspose_,
636  std::invalid_argument,
637  "SuperLU_DIST does not support solving the tranpose system" );
638 
639  data_.options.Trans = SLUD::NOTRANS; // should always be set this way;
640 
641  // TODO: Uncomment when supported
642  // bool equil = parameterList->get<bool>("Equil", true);
643  // data_.options.Equil = equil ? SLUD::YES : SLUD::NO;
644  data_.options.Equil = SLUD::NO;
645 
646  if( parameterList->isParameter("ColPerm") ){
647  RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
648  parameterList->getEntry("ColPerm").setValidator(colperm_validator);
649 
650  data_.options.ColPerm = getIntegralValue<SLUD::colperm_t>(*parameterList, "ColPerm");
651  }
652 
653  // Always use the "NOROWPERM" option to avoid a serial bottleneck
654  // with the weighted bipartite matching algorithm used for the
655  // "LargeDiag" RowPerm. Note the inconsistency with the SuperLU
656  // User guide (which states that the value should be "NATURAL").
657  data_.options.RowPerm = SLUD::NOROWPERM;
658 
659  // TODO: Uncomment when supported
660  // if( parameterList->isParameter("IterRefine") ){
661  // RCP<const ParameterEntryValidator> iter_refine_validator = valid_params->getEntry("IterRefine").validator();
662  // parameterList->getEntry("IterRefine").setValidator(iter_refine_validator);
663 
664  // data_.options.IterRefine = getIntegralValue<SLUD::IterRefine_t>(*parameterList, "IterRefine");
665  // }
666  data_.options.IterRefine = SLUD::NOREFINE;
667 
668  bool replace_tiny = parameterList->get<bool>("ReplaceTinyPivot", true);
669  data_.options.ReplaceTinyPivot = replace_tiny ? SLUD::YES : SLUD::NO;
670  }
671 
672 
673  template <class Matrix, class Vector>
674  Teuchos::RCP<const Teuchos::ParameterList>
676  {
677  using std::string;
678  using Teuchos::tuple;
679  using Teuchos::ParameterList;
680  using Teuchos::EnhancedNumberValidator;
681  using Teuchos::setStringToIntegralParameter;
682  using Teuchos::stringToIntegralParameterEntryValidator;
683 
684  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
685 
686  if( is_null(valid_params) ){
687  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
688 
689  Teuchos::RCP<EnhancedNumberValidator<SLUD::int_t> > col_row_validator
690  = Teuchos::rcp( new EnhancedNumberValidator<SLUD::int_t>() );
691  col_row_validator->setMin(1);
692 
693  pl->set("npcol", data_.grid.npcol,
694  "Number of columns in the processor grid. "
695  "Must be set with nprow", col_row_validator);
696  pl->set("nprow", data_.grid.nprow,
697  "Number of rows in the SuperLU_DIST processor grid. "
698  "Must be set together with npcol", col_row_validator);
699 
700  // validator will catch any value besides NOTRANS
701  setStringToIntegralParameter<SLUD::trans_t>("Trans", "NOTRANS",
702  "Solve for the transpose system or not",
703  tuple<string>("NOTRANS"),
704  tuple<string>("Do not solve with transpose"),
705  tuple<SLUD::trans_t>(SLUD::NOTRANS),
706  pl.getRawPtr());
707 
708  // TODO: uncomment when supported
709  // pl->set("Equil", false, "Whether to equilibrate the system before solve");
710 
711  // TODO: uncomment when supported
712  // setStringToIntegralParameter<SLUD::IterRefine_t>("IterRefine", "NOREFINE",
713  // "Type of iterative refinement to use",
714  // tuple<string>("NOREFINE", "DOUBLE"),
715  // tuple<string>("Do not use iterative refinement",
716  // "Do double iterative refinement"),
717  // tuple<SLUD::IterRefine_t>(SLUD::NOREFINE,
718  // SLUD::DOUBLE),
719  // pl.getRawPtr());
720 
721  pl->set("ReplaceTinyPivot", true,
722  "Specifies whether to replace tiny diagonals during LU factorization");
723 
724  setStringToIntegralParameter<SLUD::colperm_t>("ColPerm", "PARMETIS",
725  "Specifies how to permute the columns of the "
726  "matrix for sparsity preservation",
727  tuple<string>("NATURAL", "PARMETIS"),
728  tuple<string>("Natural ordering",
729  "ParMETIS ordering on A^T + A"),
730  tuple<SLUD::colperm_t>(SLUD::NATURAL,
731  SLUD::PARMETIS),
732  pl.getRawPtr());
733 
734  valid_params = pl;
735  }
736 
737  return valid_params;
738  }
739 
740 
741  template <class Matrix, class Vector>
742  void
744  SLUD::int_t& nprow,
745  SLUD::int_t& npcol) const {
746  TEUCHOS_TEST_FOR_EXCEPTION( nprocs < 1,
747  std::invalid_argument,
748  "Number of MPI processes must be at least 1" );
749  SLUD::int_t c, r = 1;
750  while( r*r <= nprocs ) r++;
751  nprow = npcol = --r; // fall back to square grid
752  c = nprocs / r;
753  while( (r--)*c != nprocs ){
754  c = nprocs / r; // note integer division
755  }
756  ++r;
757  // prefer the square grid over a single row (which will only happen
758  // in the case of a prime nprocs
759  if( r > 1 || nprocs < 9){ // nprocs < 9 is a heuristic for the small cases
760  nprow = r;
761  npcol = c;
762  }
763  }
764 
765 
766  template <class Matrix, class Vector>
767  bool
769  // Extract the necessary information from mat and call SLU function
770  using Teuchos::Array;
771  using Teuchos::ArrayView;
772  using Teuchos::ptrInArg;
773  using Teuchos::as;
774 
775  using SLUD::int_t;
776 
777 #ifdef HAVE_AMESOS2_TIMERS
778  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
779 #endif
780 
781  // Cleanup old store memory if it's non-NULL
782  if( data_.A.Store != NULL ){
783  SLUD::Destroy_SuperMatrix_Store_dist( &(data_.A) );
784  data_.A.Store = NULL;
785  }
786 
787  Teuchos::RCP<const MatrixAdapter<Matrix> > redist_mat
788  = this->matrixA_->get(ptrInArg(*superlu_rowmap_));
789 
790  int_t l_nnz, l_rows, g_rows, g_cols, fst_global_row;
791  l_nnz = as<int_t>(redist_mat->getLocalNNZ());
792  l_rows = as<int_t>(redist_mat->getLocalNumRows());
793  g_rows = as<int_t>(redist_mat->getGlobalNumRows());
794  g_cols = g_rows; // we deal with square matrices
795  fst_global_row = as<int_t>(superlu_rowmap_->getMinGlobalIndex());
796 
797  nzvals_.resize(l_nnz);
798  colind_.resize(l_nnz);
799  rowptr_.resize(l_rows + 1);
800 
801  int_t nnz_ret = 0;
802  {
803 #ifdef HAVE_AMESOS2_TIMERS
804  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
805 #endif
806 
809  slu_type, int_t, int_t >::do_get(redist_mat.ptr(),
810  nzvals_(), colind_(), rowptr_(),
811  nnz_ret,
812  ptrInArg(*superlu_rowmap_),
813  ARBITRARY);
814  }
815 
816  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != l_nnz,
817  std::runtime_error,
818  "Did not get the expected number of non-zero vals");
819 
820  // Get the SLU data type for this type of matrix
821  SLUD::Dtype_t dtype = type_map::dtype;
822 
823  if( in_grid_ ){
824  function_map::create_CompRowLoc_Matrix(&(data_.A),
825  g_rows, g_cols,
826  l_nnz, l_rows, fst_global_row,
827  nzvals_.getRawPtr(),
828  colind_.getRawPtr(),
829  rowptr_.getRawPtr(),
830  SLUD::SLU_NR_loc,
831  dtype, SLUD::SLU_GE);
832  }
833 
834  return true;
835 }
836 
837 
838  template<class Matrix, class Vector>
839  const char* Superludist<Matrix,Vector>::name = "SuperLU_DIST";
840 
841 
842 } // end namespace Amesos2
843 
844 #endif // AMESOS2_SUPERLUDIST_DEF_HPP
int getNumNumericFact() const
Returns the number of numeric factorizations performed by the owning solver.
Definition: Amesos2_Status.hpp:102
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
Amesos2 interface to the distributed memory version of SuperLU.
Definition: Amesos2_Superludist_decl.hpp:90
int getNumPreOrder() const
Returns the number of pre-orderings performed by the owning solver.
Definition: Amesos2_Status.hpp:94
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
Teuchos::Array< int > colind_
Stores the row indices of the nonzero entries.
Definition: Amesos2_Superludist_decl.hpp:311
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:591
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
Superludist(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Superludist_def.hpp:68
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_Superludist_def.hpp:605
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
Teuchos::Array< slu_type > nzvals_
Stores the values of the nonzero entries for SuperLU_DIST.
Definition: Amesos2_Superludist_decl.hpp:309
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_Superludist_def.hpp:675
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
Provides definition of SuperLU_DIST types as well as conversions and type traits. ...
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
bool symbolicFactorizationDone() const
If true , then symbolic factorization has been performed.
Definition: Amesos2_Status.hpp:114
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Superludist_def.hpp:302
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
Teuchos::Array< slu_type > bvals_
1D store for B values
Definition: Amesos2_Superludist_decl.hpp:315
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Superludist_def.hpp:596
bool numericFactorizationDone() const
If true , then numeric factorization has been performed.
Definition: Amesos2_Status.hpp:118
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > superlu_rowmap_
Maps rows of the matrix to processors in the SuperLU_DIST processor grid.
Definition: Amesos2_Superludist_decl.hpp:327
~Superludist()
Destructor.
Definition: Amesos2_Superludist_def.hpp:235
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:275
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using SuperLU_DIST.
Definition: Amesos2_Superludist_def.hpp:352
Teuchos::Array< int > rowptr_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_Superludist_decl.hpp:313
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
void get_default_grid_size(int nprocs, SLUD::int_t &nprow, SLUD::int_t &npcol) const
Definition: Amesos2_Superludist_def.hpp:743
bool in_grid_
true if this processor is in SuperLU_DISTS&#39;s 2D process grid
Definition: Amesos2_Superludist_decl.hpp:320
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:296
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
SuperLU_DIST specific solve.
Definition: Amesos2_Superludist_def.hpp:458
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:175
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal solver structures.
Definition: Amesos2_Superludist_def.hpp:768
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455
Teuchos::Array< slu_type > xvals_
1D store for X values
Definition: Amesos2_Superludist_decl.hpp:317
int numericFactorization_impl()
SuperLU_DIST specific numeric factorization.
Definition: Amesos2_Superludist_def.hpp:384
Status status_
Holds status information about a solver.
Definition: Amesos2_SolverCore_decl.hpp:492