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