Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_MUMPS_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Amesos2: Templated Direct Sparse Solver Package
4 //
5 // Copyright 2011 NTESS and the Amesos2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
18 #ifndef AMESOS2_MUMPS_DEF_HPP
19 #define AMESOS2_MUMPS_DEF_HPP
20 
21 #include <Teuchos_Tuple.hpp>
22 #include <Teuchos_ParameterList.hpp>
23 #include <Teuchos_StandardParameterEntryValidators.hpp>
24 
25 #ifdef HAVE_MPI
26 #include <Teuchos_DefaultMpiComm.hpp>
27 #endif
28 
29 #include <limits>
30 
32 #include "Amesos2_MUMPS_decl.hpp"
33 
34 namespace Amesos2
35 {
36 
37 
38  template <class Matrix, class Vector>
39  MUMPS<Matrix,Vector>::MUMPS(
40  Teuchos::RCP<const Matrix> A,
41  Teuchos::RCP<Vector> X,
42  Teuchos::RCP<const Vector> B )
43  : SolverCore<Amesos2::MUMPS,Matrix,Vector>(A, X, B)
44  , is_contiguous_(true)
45  {
46 
47  typedef FunctionMap<MUMPS,scalar_type> function_map;
48 
49  MUMPS_MATRIX_LOAD = false;
50  MUMPS_STRUCT = false;
51  MUMPS_MATRIX_LOAD_PREORDERING = false;
52 
53  #ifdef HAVE_MPI
54  using Teuchos::Comm;
55  using Teuchos::MpiComm;
56  using Teuchos::RCP;
57  using Teuchos::rcp;
58  using Teuchos::rcp_dynamic_cast;
59 
60  //Comm
61  mumps_par.comm_fortran = -987654;
62  RCP<const Comm<int> > matComm = this->matrixA_->getComm();
63  //Add Exception Checking
64  TEUCHOS_TEST_FOR_EXCEPTION(
65  matComm.is_null(), std::logic_error, "Amesos2::Comm");
66  RCP<const MpiComm<int> > matMpiComm =
67  rcp_dynamic_cast<const MpiComm<int> >(matComm);
68  //Add Exception Checking
69  TEUCHOS_TEST_FOR_EXCEPTION(
70  matMpiComm->getRawMpiComm().is_null(),
71  std::logic_error, "Amesos2::MPI");
72  MPI_Comm rawMpiComm = (* (matMpiComm->getRawMpiComm()) )();
73  mumps_par.comm_fortran = (int) MPI_Comm_c2f(rawMpiComm);
74  #endif
75 
76  mumps_par.job = -1;
77  mumps_par.par = 1;
78  mumps_par.sym = 0;
79 
80  function_map::mumps_c(&(mumps_par)); //init
81  MUMPS_ERROR();
82 
83  mumps_par.n = this->globalNumCols_;
84 
85  /*Default parameters*/
86  mumps_par.icntl[0] = -1; // Turn off error messages
87  mumps_par.icntl[1] = -1; // Turn off diagnositic printing
88  mumps_par.icntl[2] = -1; // Turn off global information messages
89  mumps_par.icntl[3] = 1; // No messages
90  mumps_par.icntl[4] = 0; // Matrix given in assembled (Triplet) form
91  mumps_par.icntl[5] = 7; // Choose column permuation automatically
92  mumps_par.icntl[6] = 7; // Choose ordering methods automatically
93  mumps_par.icntl[7] = 7; // Choose scaling automatically
94  mumps_par.icntl[8] = 1; // Compuate Ax = b;
95  mumps_par.icntl[9] = 0; // iterative refinement
96  mumps_par.icntl[10] = 0; //Do not collect statistics
97  mumps_par.icntl[11] = 0; // Automatic choice of ordering strategy
98  mumps_par.icntl[12] = 0; //Use ScaLAPACK for root node
99  mumps_par.icntl[13] = 20; // Increase memory allocation 20% at a time
100  mumps_par.icntl[17] = 0;
101  mumps_par.icntl[18] = 0; // do not provide back the Schur complement
102  mumps_par.icntl[19] = 0; // RHS is given in dense form
103  mumps_par.icntl[20] = 0; //RHS is in dense form
104  mumps_par.icntl[21] = 0; // Do all compuations in-core
105  mumps_par.icntl[22] = 0; // No max MB for work
106  mumps_par.icntl[23] = 0; // Do not perform null pivot detection
107  mumps_par.icntl[24] = 0; // No null space basis compuation
108  mumps_par.icntl[25] = 0; // Do not condense/reduce Schur RHS
109  mumps_par.icntl[27] = 1; // sequential analysis
110  mumps_par.icntl[28] = 0; //
111  mumps_par.icntl[29] = 0; //
112  mumps_par.icntl[30] = 0; //
113  mumps_par.icntl[31] = 0;
114  mumps_par.icntl[32] = 0;
115  }
116 
117  template <class Matrix, class Vector>
118  MUMPS<Matrix,Vector>::~MUMPS( )
119  {
120  /* Clean up the struc*/
121  typedef FunctionMap<MUMPS,scalar_type> function_map;
122 
123  if(MUMPS_STRUCT == true)
124  {
125  free(mumps_par.a);
126  free(mumps_par.jcn);
127  free(mumps_par.irn);
128  }
129  mumps_par.job = -2;
130  if (this->rank_ < this->nprocs_) {
131  function_map::mumps_c(&(mumps_par));
132  }
133 
134  }
135 
136  template<class Matrix, class Vector>
137  int
139  {
140  /* TODO: Define what it means for MUMPS
141  */
142  #ifdef HAVE_AMESOS2_TIMERS
143  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
144  #endif
145 
146  return(0);
147  }//end preOrdering_impl()
148 
149  template <class Matrix, class Vector>
150  int
152  {
153  #ifdef HAVE_AMESOS2_TIMERS
154  Teuchos::TimeMonitor symFactTime( this->timers_.symFactTime_ );
155  #endif
156 
157  typedef FunctionMap<MUMPS,scalar_type> function_map;
158  if ( this->globalNumRows_ > 0 ) {
159  mumps_par.par = 1;
160  mumps_par.job = 1; // sym factor
161  function_map::mumps_c(&(mumps_par));
162  MUMPS_ERROR();
163  }
164  return(0);
165  }//end symblicFactortion_impl()
166 
167 
168  template <class Matrix, class Vector>
169  int
171  {
172  #ifdef HAVE_AMESOS2_TIMERS
173  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
174  #endif
175 
176  typedef FunctionMap<MUMPS,scalar_type> function_map;
177  if ( this->globalNumRows_ > 0 ) {
178  mumps_par.job = 2;
179  function_map::mumps_c(&(mumps_par));
180  MUMPS_ERROR();
181  }
182  return(0);
183  }//end numericFactorization_impl()
184 
185 
186  template <class Matrix, class Vector>
187  int
189  const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
190  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
191  {
192  typedef FunctionMap<MUMPS,scalar_type> function_map;
193 
194  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
195  const size_t nrhs = X->getGlobalNumVectors();
196  const size_t val_store_size = Teuchos::as<size_t>(ld_rhs * nrhs);
197  {
198  #ifdef HAVE_AMESOS2_TIMERS
199  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
200  #endif
201 
202  xvals_.resize(val_store_size);
203  bvals_.resize(val_store_size);
204 
206  mumps_type>::do_get(B, bvals_(), Teuchos::as<size_t>(ld_rhs),
207  (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
208  this->rowIndexBase_);
209  }
210  int ierr = 0; // returned error code
211  if ( this->globalNumRows_ > 0 ) {
212  mumps_par.nrhs = nrhs;
213  mumps_par.lrhs = mumps_par.n;
214  mumps_par.job = 3;
215 
216  if ( this->root_ )
217  {
218  mumps_par.rhs = bvals_.getRawPtr();
219  }
220 
221  #ifdef HAVE_AMESOS2_TIMERS
222  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
223  #endif
224 
225  function_map::mumps_c(&(mumps_par));
226  MUMPS_ERROR();
227  }
228  {
229  #ifdef HAVE_AMESOS2_TIMERS
230  Teuchos::TimeMonitor redistTimer2(this->timers_.vecRedistTime_);
231  #endif
232 
234  MultiVecAdapter<Vector>,mumps_type>::do_put(X, bvals_(), Teuchos::as<size_t>(ld_rhs),
235  (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
236  }
237  // ch: see function loadA_impl()
238  MUMPS_MATRIX_LOAD_PREORDERING = false;
239  return(ierr);
240  }//end solve()
241 
242 
243  template <class Matrix, class Vector>
244  bool
246  {
247  // The MUMPS can only handle square for right now
248  return( this->globalNumRows_ == this->globalNumCols_ );
249  }
250 
251 
252  template <class Matrix, class Vector>
253  void
254  MUMPS<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
255  {
256  using Teuchos::RCP;
257  using Teuchos::getIntegralValue;
258  using Teuchos::ParameterEntryValidator;
259 
260  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
261  /*To Do --- add support for parameters */
262  if(parameterList->isParameter("ICNTL(1)")){
263  mumps_par.icntl[0] = parameterList->get<int>("ICNTL(1)", -1);
264  }
265  if(parameterList->isParameter("ICNTL(2)")){
266  mumps_par.icntl[1] = parameterList->get<int>("ICNTL(2)", -1);
267  }
268  if(parameterList->isParameter("ICNTL(3)")){
269  mumps_par.icntl[2] = parameterList->get<int>("ICNTL(3)", -1);
270  }
271  if(parameterList->isParameter("ICNTL(4)")){
272  mumps_par.icntl[3] = parameterList->get<int>("ICNTL(4)", 1);
273  }
274  if(parameterList->isParameter("ICNTL(6)")){
275  mumps_par.icntl[5] = parameterList->get<int>("ICNTL(6)", 0);
276  }
277  if(parameterList->isParameter("ICNTL(7)")){
278  mumps_par.icntl[6] = parameterList->get<int>("ICNTL(7)", 7);
279  }
280  if(parameterList->isParameter("ICNTL(9)")){
281  mumps_par.icntl[8] = parameterList->get<int>("ICNTL(9)", 1);
282  }
283  if(parameterList->isParameter("ICNTL(11)")){
284  mumps_par.icntl[10] = parameterList->get<int>("ICNTL(11)", 0);
285  }
286  if(parameterList->isParameter("ICNTL(14)")){
287  mumps_par.icntl[13] = parameterList->get<int>("ICNTL(14)", 20);
288  }
289  if( parameterList->isParameter("IsContiguous") ){
290  is_contiguous_ = parameterList->get<bool>("IsContiguous");
291  }
292  }//end set parameters()
293 
294 
295  template <class Matrix, class Vector>
296  Teuchos::RCP<const Teuchos::ParameterList>
298  {
299  using Teuchos::ParameterList;
300 
301  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
302 
303  if( is_null(valid_params) ){
304  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
305 
306  pl->set("ICNTL(1)", -1, "See Manual" );
307  pl->set("ICNTL(2)", -1, "See Manual" );
308  pl->set("ICNTL(3)", -1, "See Manual" );
309  pl->set("ICNTL(4)", 1, "See Manual" );
310  pl->set("ICNTL(6)", 0, "See Manual" );
311  pl->set("ICNTL(9)", 1, "See Manual" );
312  pl->set("ICNTL(11)", 0, "See Manual" );
313  pl->set("ICNTL(14)", 20, "See Manual" );
314  pl->set("IsContiguous", true, "Whether GIDs contiguous");
315 
316  valid_params = pl;
317  }
318 
319  return valid_params;
320  }//end getValidParmaeters_impl()
321 
322 
323  template <class Matrix, class Vector>
324  bool
326  {
327  #ifdef HAVE_AMESOS2_TIMERS
328  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
329  #endif
330  if(MUMPS_MATRIX_LOAD == false || current_phase==NUMFACT)
331  {
332  // Only the root image needs storage allocated
333  if( !MUMPS_MATRIX_LOAD && this->root_ ){
334  Kokkos::resize(host_nzvals_view_, this->globalNumNonZeros_);
335  Kokkos::resize(host_rows_view_, this->globalNumNonZeros_);
336  Kokkos::resize(host_col_ptr_view_, this->globalNumRows_ + 1);
337  }
338 
339  local_ordinal_type nnz_ret = 0;
340 
341  #ifdef HAVE_AMESOS2_TIMERS
342  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
343  #endif
344 
346  MatrixAdapter<Matrix>,host_value_type_view,host_ordinal_type_view,host_ordinal_type_view>
347  ::do_get(this->matrixA_.ptr(), host_nzvals_view_, host_rows_view_, host_col_ptr_view_, nnz_ret,
348  (is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
349  ARBITRARY,
350  this->rowIndexBase_);
351 
352  if( this->root_ ){
353  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != Teuchos::as<local_ordinal_type>(this->globalNumNonZeros_),
354  std::runtime_error,
355  "Did not get the expected number of non-zero vals");
356  }
357 
358 
359  if( this->root_ ){
360  ConvertToTriplet();
361  }
362  /* ch: In general, the matrix is loaded during the preordering phase.
363  However, if the matrix pattern has not changed during consecutive calls of numeric factorizations
364  we can reuse the previous symbolic factorization. In this case, the matrix is not loaded in the preordering phase,
365  because it is not called. Therefore, we need to load the matrix in the numeric factorization phase.
366  */
367  if (current_phase==PREORDERING){
368  MUMPS_MATRIX_LOAD_PREORDERING = true;
369  }
370  }
371 
372 
373 
374  MUMPS_MATRIX_LOAD = true;
375  return (true);
376  }//end loadA_impl()
377 
378  template <class Matrix, class Vector>
379  int
381  {
382  if ( !MUMPS_STRUCT ) {
383  MUMPS_STRUCT = true;
384  mumps_par.n = this->globalNumCols_;
385  mumps_par.nz = this->globalNumNonZeros_;
386  mumps_par.a = (mumps_type*)malloc(mumps_par.nz * sizeof(mumps_type));
387  mumps_par.irn = (MUMPS_INT*)malloc(mumps_par.nz *sizeof(MUMPS_INT));
388  mumps_par.jcn = (MUMPS_INT*)malloc(mumps_par.nz * sizeof(MUMPS_INT));
389  }
390  if((mumps_par.a == NULL) || (mumps_par.irn == NULL)
391  || (mumps_par.jcn == NULL))
392  {
393  return -1;
394  }
395  /* Going from full CSC to full Triplet */
396  /* Will have to add support for symmetric case*/
397  local_ordinal_type tri_count = 0;
398  local_ordinal_type i,j;
399  local_ordinal_type max_local_ordinal = 0;
400 
401  for(i = 0; i < (local_ordinal_type)this->globalNumCols_; i++)
402  {
403  for( j = host_col_ptr_view_(i); j < host_col_ptr_view_(i+1)-1; j++)
404  {
405  mumps_par.jcn[tri_count] = (MUMPS_INT)i+1; //Fortran index
406  mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1; //Fortran index
407  mumps_par.a[tri_count] = host_nzvals_view_(j);
408 
409  tri_count++;
410  }
411 
412  j = host_col_ptr_view_(i+1)-1;
413  mumps_par.jcn[tri_count] = (MUMPS_INT)i+1; //Fortran index
414  mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1; //Fortran index
415  mumps_par.a[tri_count] = host_nzvals_view_(j);
416 
417  tri_count++;
418 
419  if(host_rows_view_(j) > max_local_ordinal)
420  {
421  max_local_ordinal = host_rows_view_(j);
422  }
423  }
424  TEUCHOS_TEST_FOR_EXCEPTION(std::numeric_limits<MUMPS_INT>::max() <= max_local_ordinal,
425  std::runtime_error,
426  "Matrix index larger than MUMPS_INT");
427 
428  return 0;
429  }//end Convert to Trip()
430 
431  template<class Matrix, class Vector>
432  void
433  MUMPS<Matrix,Vector>::MUMPS_ERROR()const
434  {
435  using Teuchos::Comm;
436  using Teuchos::RCP;
437  bool Wrong = ((mumps_par.info[0] != 0) || (mumps_par.infog[0] != 0)) && (this->rank_ < this->nprocs_);
438  if(Wrong){
439  if (this->rank_==0) {
440  std::cerr << "Amesos2_Mumps : ERROR" << std::endl;
441  if ( this->status_.getNumSolve() > 0) {
442  std::cerr << " Last Phase : SOLVE" << std::endl;
443  } else if( this->status_.numericFactorizationDone() ){
444  std::cerr << " Last Phase : NUMFACT" << std::endl;
445  } else if( this->status_.symbolicFactorizationDone() ){
446  std::cerr << " Last Phase : SYMBFACT" << std::endl;
447  } else if( this->status_.preOrderingDone() ){
448  std::cerr << " Last Phase : PREORDERING" << std::endl;
449  } else {
450  std::cerr << " Last Phase : CLEAN" << std::endl;
451  }
452  std::cerr << "Amesos2_Mumps : INFOG(1) = " << mumps_par.infog[0] << std::endl;
453  std::cerr << "Amesos2_Mumps : INFOG(2) = " << mumps_par.infog[1] << std::endl;
454  }
455  if (mumps_par.info[0] != 0 && Wrong) {
456  std::cerr << "Amesos2_Mumps : On process " << this->matrixA_->getComm()->getRank()
457  << ", INFO(1) = " << mumps_par.info[0] << std::endl;
458  std::cerr << "Amesos2_Mumps : On process " << this->matrixA_->getComm()->getRank()
459  << ", INFO(2) = " << mumps_par.info[1] << std::endl;
460  }
461 
462 
463  }
464  // Throw on all ranks
465  int WrongInt = Wrong;
466  RCP<const Comm<int> > matComm = this->matrixA_->getComm();
467  Teuchos::broadcast<int,int>(*matComm,0,1,&WrongInt);
468  TEUCHOS_TEST_FOR_EXCEPTION(WrongInt>0,
469  std::runtime_error,
470  "MUMPS error");
471 
472  }//end MUMPS_ERROR()
473 
474 
475  template<class Matrix, class Vector>
476  const char* MUMPS<Matrix,Vector>::name = "MUMPS";
477 
478 } // end namespace Amesos2
479 
480 #endif // AMESOS2_MUMPS_DEF_HPP
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:632
int numericFactorization_impl()
MUMPS specific numeric factorization.
Definition: Amesos2_MUMPS_def.hpp:170
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:31
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:445
Amesos2 MUMPS declarations.
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:233
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_MUMPS_def.hpp:138
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
MUMPS specific solve.
Definition: Amesos2_MUMPS_def.hpp:188
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_MUMPS_def.hpp:245
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_MUMPS_def.hpp:325
Amesos2 interface to the MUMPS package.
Definition: Amesos2_MUMPS_decl.hpp:50
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:42
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_MUMPS_def.hpp:297
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_MUMPS_def.hpp:254
Passes functions to TPL functions based on type.
Definition: Amesos2_FunctionMap.hpp:42
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:339
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:142
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:421