18 #ifndef AMESOS2_MUMPS_DEF_HPP
19 #define AMESOS2_MUMPS_DEF_HPP
21 #include <Teuchos_Tuple.hpp>
22 #include <Teuchos_ParameterList.hpp>
23 #include <Teuchos_StandardParameterEntryValidators.hpp>
26 #include <Teuchos_DefaultMpiComm.hpp>
38 template <
class Matrix,
class Vector>
39 MUMPS<Matrix,Vector>::MUMPS(
40 Teuchos::RCP<const Matrix> A,
41 Teuchos::RCP<Vector> X,
42 Teuchos::RCP<const Vector> B )
43 : SolverCore<Amesos2::MUMPS,Matrix,Vector>(A, X, B)
44 , is_contiguous_(true)
47 typedef FunctionMap<MUMPS,scalar_type> function_map;
49 MUMPS_MATRIX_LOAD =
false;
51 MUMPS_MATRIX_LOAD_PREORDERING =
false;
55 using Teuchos::MpiComm;
58 using Teuchos::rcp_dynamic_cast;
61 mumps_par.comm_fortran = -987654;
62 RCP<const Comm<int> > matComm = this->
matrixA_->getComm();
64 TEUCHOS_TEST_FOR_EXCEPTION(
65 matComm.is_null(), std::logic_error,
"Amesos2::Comm");
66 RCP<const MpiComm<int> > matMpiComm =
67 rcp_dynamic_cast<
const MpiComm<int> >(matComm);
69 TEUCHOS_TEST_FOR_EXCEPTION(
70 matMpiComm->getRawMpiComm().is_null(),
71 std::logic_error,
"Amesos2::MPI");
72 MPI_Comm rawMpiComm = (* (matMpiComm->getRawMpiComm()) )();
73 mumps_par.comm_fortran = (int) MPI_Comm_c2f(rawMpiComm);
80 function_map::mumps_c(&(mumps_par));
86 mumps_par.icntl[0] = -1;
87 mumps_par.icntl[1] = -1;
88 mumps_par.icntl[2] = -1;
89 mumps_par.icntl[3] = 1;
90 mumps_par.icntl[4] = 0;
91 mumps_par.icntl[5] = 7;
92 mumps_par.icntl[6] = 7;
93 mumps_par.icntl[7] = 7;
94 mumps_par.icntl[8] = 1;
95 mumps_par.icntl[9] = 0;
96 mumps_par.icntl[10] = 0;
97 mumps_par.icntl[11] = 0;
98 mumps_par.icntl[12] = 0;
99 mumps_par.icntl[13] = 20;
100 mumps_par.icntl[17] = 0;
101 mumps_par.icntl[18] = 0;
102 mumps_par.icntl[19] = 0;
103 mumps_par.icntl[20] = 0;
104 mumps_par.icntl[21] = 0;
105 mumps_par.icntl[22] = 0;
106 mumps_par.icntl[23] = 0;
107 mumps_par.icntl[24] = 0;
108 mumps_par.icntl[25] = 0;
109 mumps_par.icntl[27] = 1;
110 mumps_par.icntl[28] = 0;
111 mumps_par.icntl[29] = 0;
112 mumps_par.icntl[30] = 0;
113 mumps_par.icntl[31] = 0;
114 mumps_par.icntl[32] = 0;
117 template <
class Matrix,
class Vector>
118 MUMPS<Matrix,Vector>::~MUMPS( )
121 typedef FunctionMap<MUMPS,scalar_type> function_map;
123 if(MUMPS_STRUCT ==
true)
130 if (this->rank_ < this->nprocs_) {
131 function_map::mumps_c(&(mumps_par));
136 template<
class Matrix,
class Vector>
142 #ifdef HAVE_AMESOS2_TIMERS
143 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
149 template <
class Matrix,
class Vector>
153 #ifdef HAVE_AMESOS2_TIMERS
154 Teuchos::TimeMonitor symFactTime( this->timers_.symFactTime_ );
158 if ( this->globalNumRows_ > 0 ) {
161 function_map::mumps_c(&(mumps_par));
168 template <
class Matrix,
class Vector>
172 #ifdef HAVE_AMESOS2_TIMERS
173 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
177 if ( this->globalNumRows_ > 0 ) {
179 function_map::mumps_c(&(mumps_par));
186 template <
class Matrix,
class Vector>
194 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
195 const size_t nrhs = X->getGlobalNumVectors();
196 const size_t val_store_size = Teuchos::as<size_t>(ld_rhs * nrhs);
198 #ifdef HAVE_AMESOS2_TIMERS
199 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
202 xvals_.resize(val_store_size);
203 bvals_.resize(val_store_size);
206 mumps_type>::do_get(B, bvals_(), Teuchos::as<size_t>(ld_rhs),
207 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
208 this->rowIndexBase_);
211 if ( this->globalNumRows_ > 0 ) {
212 mumps_par.nrhs = nrhs;
213 mumps_par.lrhs = mumps_par.n;
218 mumps_par.rhs = bvals_.getRawPtr();
221 #ifdef HAVE_AMESOS2_TIMERS
222 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
225 function_map::mumps_c(&(mumps_par));
229 #ifdef HAVE_AMESOS2_TIMERS
230 Teuchos::TimeMonitor redistTimer2(this->timers_.vecRedistTime_);
235 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
238 MUMPS_MATRIX_LOAD_PREORDERING =
false;
243 template <
class Matrix,
class Vector>
248 return( this->globalNumRows_ == this->globalNumCols_ );
252 template <
class Matrix,
class Vector>
257 using Teuchos::getIntegralValue;
258 using Teuchos::ParameterEntryValidator;
260 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
262 if(parameterList->isParameter(
"ICNTL(1)")){
263 mumps_par.icntl[0] = parameterList->get<
int>(
"ICNTL(1)", -1);
265 if(parameterList->isParameter(
"ICNTL(2)")){
266 mumps_par.icntl[1] = parameterList->get<
int>(
"ICNTL(2)", -1);
268 if(parameterList->isParameter(
"ICNTL(3)")){
269 mumps_par.icntl[2] = parameterList->get<
int>(
"ICNTL(3)", -1);
271 if(parameterList->isParameter(
"ICNTL(4)")){
272 mumps_par.icntl[3] = parameterList->get<
int>(
"ICNTL(4)", 1);
274 if(parameterList->isParameter(
"ICNTL(6)")){
275 mumps_par.icntl[5] = parameterList->get<
int>(
"ICNTL(6)", 0);
277 if(parameterList->isParameter(
"ICNTL(7)")){
278 mumps_par.icntl[6] = parameterList->get<
int>(
"ICNTL(7)", 7);
280 if(parameterList->isParameter(
"ICNTL(9)")){
281 mumps_par.icntl[8] = parameterList->get<
int>(
"ICNTL(9)", 1);
283 if(parameterList->isParameter(
"ICNTL(11)")){
284 mumps_par.icntl[10] = parameterList->get<
int>(
"ICNTL(11)", 0);
286 if(parameterList->isParameter(
"ICNTL(14)")){
287 mumps_par.icntl[13] = parameterList->get<
int>(
"ICNTL(14)", 20);
289 if( parameterList->isParameter(
"IsContiguous") ){
290 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
295 template <
class Matrix,
class Vector>
296 Teuchos::RCP<const Teuchos::ParameterList>
299 using Teuchos::ParameterList;
301 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
303 if( is_null(valid_params) ){
304 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
306 pl->set(
"ICNTL(1)", -1,
"See Manual" );
307 pl->set(
"ICNTL(2)", -1,
"See Manual" );
308 pl->set(
"ICNTL(3)", -1,
"See Manual" );
309 pl->set(
"ICNTL(4)", 1,
"See Manual" );
310 pl->set(
"ICNTL(6)", 0,
"See Manual" );
311 pl->set(
"ICNTL(9)", 1,
"See Manual" );
312 pl->set(
"ICNTL(11)", 0,
"See Manual" );
313 pl->set(
"ICNTL(14)", 20,
"See Manual" );
314 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
323 template <
class Matrix,
class Vector>
327 #ifdef HAVE_AMESOS2_TIMERS
328 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
330 if(MUMPS_MATRIX_LOAD ==
false || current_phase==NUMFACT)
333 if( !MUMPS_MATRIX_LOAD && this->root_ ){
334 Kokkos::resize(host_nzvals_view_, this->globalNumNonZeros_);
335 Kokkos::resize(host_rows_view_, this->globalNumNonZeros_);
336 Kokkos::resize(host_col_ptr_view_, this->globalNumRows_ + 1);
339 local_ordinal_type nnz_ret = 0;
341 #ifdef HAVE_AMESOS2_TIMERS
342 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
347 ::do_get(this->matrixA_.ptr(), host_nzvals_view_, host_rows_view_, host_col_ptr_view_, nnz_ret,
348 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
350 this->rowIndexBase_);
353 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != Teuchos::as<local_ordinal_type>(this->globalNumNonZeros_),
355 "Did not get the expected number of non-zero vals");
367 if (current_phase==PREORDERING){
368 MUMPS_MATRIX_LOAD_PREORDERING =
true;
374 MUMPS_MATRIX_LOAD =
true;
378 template <
class Matrix,
class Vector>
382 if ( !MUMPS_STRUCT ) {
384 mumps_par.n = this->globalNumCols_;
385 mumps_par.nz = this->globalNumNonZeros_;
386 mumps_par.a = (mumps_type*)malloc(mumps_par.nz *
sizeof(mumps_type));
387 mumps_par.irn = (MUMPS_INT*)malloc(mumps_par.nz *
sizeof(MUMPS_INT));
388 mumps_par.jcn = (MUMPS_INT*)malloc(mumps_par.nz *
sizeof(MUMPS_INT));
390 if((mumps_par.a == NULL) || (mumps_par.irn == NULL)
391 || (mumps_par.jcn == NULL))
397 local_ordinal_type tri_count = 0;
398 local_ordinal_type i,j;
399 local_ordinal_type max_local_ordinal = 0;
401 for(i = 0; i < (local_ordinal_type)this->globalNumCols_; i++)
403 for( j = host_col_ptr_view_(i); j < host_col_ptr_view_(i+1)-1; j++)
405 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1;
406 mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1;
407 mumps_par.a[tri_count] = host_nzvals_view_(j);
412 j = host_col_ptr_view_(i+1)-1;
413 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1;
414 mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1;
415 mumps_par.a[tri_count] = host_nzvals_view_(j);
419 if(host_rows_view_(j) > max_local_ordinal)
421 max_local_ordinal = host_rows_view_(j);
424 TEUCHOS_TEST_FOR_EXCEPTION(std::numeric_limits<MUMPS_INT>::max() <= max_local_ordinal,
426 "Matrix index larger than MUMPS_INT");
431 template<
class Matrix,
class Vector>
433 MUMPS<Matrix,Vector>::MUMPS_ERROR()
const
437 bool Wrong = ((mumps_par.info[0] != 0) || (mumps_par.infog[0] != 0)) && (this->rank_ < this->nprocs_);
439 if (this->rank_==0) {
440 std::cerr <<
"Amesos2_Mumps : ERROR" << std::endl;
441 if ( this->status_.getNumSolve() > 0) {
442 std::cerr <<
" Last Phase : SOLVE" << std::endl;
443 }
else if( this->status_.numericFactorizationDone() ){
444 std::cerr <<
" Last Phase : NUMFACT" << std::endl;
445 }
else if( this->status_.symbolicFactorizationDone() ){
446 std::cerr <<
" Last Phase : SYMBFACT" << std::endl;
447 }
else if( this->status_.preOrderingDone() ){
448 std::cerr <<
" Last Phase : PREORDERING" << std::endl;
450 std::cerr <<
" Last Phase : CLEAN" << std::endl;
452 std::cerr <<
"Amesos2_Mumps : INFOG(1) = " << mumps_par.infog[0] << std::endl;
453 std::cerr <<
"Amesos2_Mumps : INFOG(2) = " << mumps_par.infog[1] << std::endl;
455 if (mumps_par.info[0] != 0 && Wrong) {
456 std::cerr <<
"Amesos2_Mumps : On process " << this->matrixA_->getComm()->getRank()
457 <<
", INFO(1) = " << mumps_par.info[0] << std::endl;
458 std::cerr <<
"Amesos2_Mumps : On process " << this->matrixA_->getComm()->getRank()
459 <<
", INFO(2) = " << mumps_par.info[1] << std::endl;
465 int WrongInt = Wrong;
466 RCP<const Comm<int> > matComm = this->matrixA_->getComm();
467 Teuchos::broadcast<int,int>(*matComm,0,1,&WrongInt);
468 TEUCHOS_TEST_FOR_EXCEPTION(WrongInt>0,
475 template<
class Matrix,
class Vector>
476 const char* MUMPS<Matrix,Vector>::name =
"MUMPS";
480 #endif // AMESOS2_MUMPS_DEF_HPP
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:632
int numericFactorization_impl()
MUMPS specific numeric factorization.
Definition: Amesos2_MUMPS_def.hpp:170
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:31
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:445
Amesos2 MUMPS declarations.
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:233
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_MUMPS_def.hpp:138
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:188
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_MUMPS_def.hpp:245
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_MUMPS_def.hpp:325
Amesos2 interface to the MUMPS package.
Definition: Amesos2_MUMPS_decl.hpp:50
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:42
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_MUMPS_def.hpp:297
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_MUMPS_def.hpp:254
Passes functions to TPL functions based on type.
Definition: Amesos2_FunctionMap.hpp:42
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:339
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:142
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:421