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