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