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 
65 #include "Amesos2_MUMPS_decl.hpp"
66 
67 namespace Amesos2
68 {
69 
70 
71  template <class Matrix, class Vector>
72  MUMPS<Matrix,Vector>::MUMPS(
73  Teuchos::RCP<const Matrix> A,
74  Teuchos::RCP<Vector> X,
75  Teuchos::RCP<const Vector> B )
76  : SolverCore<Amesos2::MUMPS,Matrix,Vector>(A, X, B)
77  , nzvals_() // initialize to empty arrays
78  , rowind_()
79  , colptr_()
80  {
81 
82  typedef FunctionMap<MUMPS,scalar_type> function_map;
83 
84  MUMPS_MATRIX_LOAD = false;
85  MUMPS_STRUCT = 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  if(MUMPS_STRUCT == true)
157  {
158  free(mumps_par.a);
159  free(mumps_par.jcn);
160  free(mumps_par.irn);
161  }
162  }
163 
164  template<class Matrix, class Vector>
165  int
167  {
168  /* TODO: Define what it means for MUMPS
169  */
170  #ifdef HAVE_AMESOS2_TIMERS
171  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
172  #endif
173 
174  return(0);
175  }//end preOrdering_impl()
176 
177  template <class Matrix, class Vector>
178  int
180  {
181 
182  typedef FunctionMap<MUMPS,scalar_type> function_map;
183 
184  mumps_par.par = 1;
185  mumps_par.job = 1; // sym factor
186  function_map::mumps_c(&(mumps_par));
187  MUMPS_ERROR();
188 
189  return(0);
190  }//end symblicFactortion_impl()
191 
192 
193  template <class Matrix, class Vector>
194  int
196  {
197  using Teuchos::as;
198  typedef FunctionMap<MUMPS,scalar_type> function_map;
199 
200  if ( this->root_ )
201  {
202  { // Do factorization
203  #ifdef HAVE_AMESOS2_TIMERS
204  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
205  #endif
206 
207  #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
208  std::cout << "MUMPS:: Before numeric factorization" << std::endl;
209  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
210  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
211  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
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
247  slu_type>::do_get(B, bvals_(),as<size_t>(ld_rhs),
248  ROOTED);
249 
250 
251  int ierr = 0; // returned error code
252  mumps_par.nrhs = nrhs;
253  mumps_par.lrhs = mumps_par.n;
254  mumps_par.job = 3;
255 
256  if ( this->root_ )
257  {
258  mumps_par.rhs = bvals_.getRawPtr();
259  }
260 
261  #ifdef HAVE_AMESOS2_TIMERS
262  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
263  #endif
264 
265  function_map::mumps_c(&(mumps_par));
266  MUMPS_ERROR();
267 
268 
269  #ifdef HAVE_AMESOS2_TIMERS
270  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
271  #endif
272 
274  MultiVecAdapter<Vector>,slu_type>::do_put(X, bvals_(),
275  as<size_t>(ld_rhs),
276  ROOTED);
277 
278  return(ierr);
279  }//end solve()
280 
281 
282  template <class Matrix, class Vector>
283  bool
285  {
286  // The Basker can only handle square for right now
287  return( this->globalNumRows_ == this->globalNumCols_ );
288  }
289 
290 
291  template <class Matrix, class Vector>
292  void
293  MUMPS<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
294  {
295  using Teuchos::RCP;
296  using Teuchos::getIntegralValue;
297  using Teuchos::ParameterEntryValidator;
298 
299  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
300  /*To Do --- add support for parameters */
301 
302  if(parameterList->isParameter("ICNTL(1)"))
303  {
304  mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
305  "ICNTL(1)");
306  }
307  if(parameterList->isParameter("ICNTL(2)"))
308  {
309  mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
310  "ICNTL(2)");
311  }
312  if(parameterList->isParameter("ICNTL(3)"))
313  {
314  mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
315  "ICNTL(3)");
316  }
317  if(parameterList->isParameter("ICNTL(4)"))
318  {
319  mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
320  "ICNTL(4)");
321  }
322  if(parameterList->isParameter("ICNTL(6)"))
323  {
324  mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
325  "ICNTL(6)");
326  }
327  if(parameterList->isParameter("ICNTL(9)"))
328  {
329  mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
330  "ICNTL(9)");
331  }
332  if(parameterList->isParameter("ICNTL(11)"))
333  {
334  mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
335  "ICNTL(11)");
336  }
337  }//end set parameters()
338 
339 
340  template <class Matrix, class Vector>
341  Teuchos::RCP<const Teuchos::ParameterList>
343  {
344  using Teuchos::ParameterList;
345 
346  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
347 
348  if( is_null(valid_params) ){
349  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
350 
351  pl->set("ICNTL(1)", "no", "See Manual" );
352  pl->set("ICNTL(2)", "no", "See Manual" );
353  pl->set("ICNTL(3)", "no", "See Manual" );
354  pl->set("ICNTL(4)", "no", "See Manual" );
355  pl->set("ICNTL(6)", "no", "See Manual" );
356  pl->set("ICNTL(9)", "no", "See Manual" );
357  pl->set("ICNTL(11)", "no", "See Manual" );
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 
376  if(MUMPS_MATRIX_LOAD == false)
377  {
378  // Only the root image needs storage allocated
379  if( this->root_ ){
380  nzvals_.resize(this->globalNumNonZeros_);
381  rowind_.resize(this->globalNumNonZeros_);
382  colptr_.resize(this->globalNumCols_ + 1);
383  }
384 
385  local_ordinal_type nnz_ret = 0;
386 
387  #ifdef HAVE_AMESOS2_TIMERS
388  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
389  #endif
390 
392  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
393  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
394  nnz_ret, ROOTED, ARBITRARY);
395 
396  if( this->root_ ){
397  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
398  std::runtime_error,
399  "Did not get the expected number of non-zero vals");
400  }
401 
402 
403  if( this->root_ ){
404  ConvertToTriplet();
405  }
406  }
407 
408  MUMPS_MATRIX_LOAD = true;
409  return (true);
410  }//end loadA_impl()
411 
412  template <class Matrix, class Vector>
413  int
415  {
416  MUMPS_STRUCT = true;
417  mumps_par.n = this->globalNumCols_;
418  mumps_par.nz = this->globalNumNonZeros_;
419  mumps_par.a = (magnitude_type*)malloc(mumps_par.nz * sizeof(magnitude_type));
420  mumps_par.irn = (local_ordinal_type*)malloc(mumps_par.nz *sizeof(local_ordinal_type));
421  mumps_par.jcn = (local_ordinal_type*)malloc(mumps_par.nz * sizeof(local_ordinal_type));
422 
423  if((mumps_par.a == NULL) || (mumps_par.irn == NULL)
424  || (mumps_par.jcn == NULL))
425  {
426  return -1;
427  }
428  /* Going from full CSC to full Triplet */
429  /* Will have to add support for symmetric case*/
430  local_ordinal_type tri_count = 0;
431  local_ordinal_type i,j;
432 
433  for(i = 0; i < (local_ordinal_type)this->globalNumCols_; i++)
434  {
435  for( j = colptr_[i]; j < colptr_[i+1]; j++)
436  {
437  mumps_par.jcn[tri_count] = (local_ordinal_type)i+1; //Fortran index
438  mumps_par.irn[tri_count] = (local_ordinal_type)rowind_[j]+1; //Fortran index
439  mumps_par.a[tri_count] = nzvals_[j];
440  tri_count++;
441  }
442  }
443 
444  return 0;
445  }//end Convert to Trip()
446 
447  template<class Matrix, class Vector>
448  void
450  {
451  if(mumps_par.info[0] < 0)
452  {
453  TEUCHOS_TEST_FOR_EXCEPTION(false,
454  std::runtime_error,
455  "MUMPS error");
456  }
457  }//end MUMPS_ERROR()
458 
459 
460  template<class Matrix, class Vector>
461  const char* MUMPS<Matrix,Vector>::name = "MUMPS";
462 
463 } // end namespace Amesos2
464 
465 #endif // AMESOS2_MUMPS_DEF_HPP
int numericFactorization_impl()
Basker specific numeric factorization.
Definition: Amesos2_MUMPS_def.hpp:195
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
bool root_
If true, then this is the root processor.
Definition: Amesos2_SolverCore_decl.hpp:507
Teuchos::Array< local_ordinal_type > rowind_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_MUMPS_decl.hpp:207
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
Amesos2 MUMPS declarations.
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
Definition: Amesos2_TypeDecl.hpp:142
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_MUMPS_def.hpp:166
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
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
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_MUMPS_def.hpp:284
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:78
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
Teuchos::Array< slu_type > xvals_
Persisting 1D store for X.
Definition: Amesos2_MUMPS_decl.hpp:212
Teuchos::Array< slu_type > bvals_
Persisting 1D store for B.
Definition: Amesos2_MUMPS_decl.hpp:214
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_MUMPS_def.hpp:342
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_MUMPS_def.hpp:293
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
Passes functions to TPL functions based on type.
Definition: Amesos2_FunctionMap.hpp:76
Teuchos::Array< slu_type > nzvals_
Stores the values of the nonzero entries for MUMPS.
Definition: Amesos2_MUMPS_decl.hpp:205
Definition: Amesos2_TypeDecl.hpp:127
Teuchos::Array< local_ordinal_type > colptr_
Stores the row indices of the nonzero entries.
Definition: Amesos2_MUMPS_decl.hpp:209
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