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