52 #ifndef AMESOS2_MUMPS_DEF_HPP
53 #define AMESOS2_MUMPS_DEF_HPP
55 #include <Teuchos_Tuple.hpp>
56 #include <Teuchos_ParameterList.hpp>
57 #include <Teuchos_StandardParameterEntryValidators.hpp>
60 #include <Teuchos_DefaultMpiComm.hpp>
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 , is_contiguous_(true)
81 typedef FunctionMap<MUMPS,scalar_type> function_map;
83 MUMPS_MATRIX_LOAD =
false;
85 MUMPS_MATRIX_LOAD_PREORDERING =
false;
89 using Teuchos::MpiComm;
92 using Teuchos::rcp_dynamic_cast;
95 mumps_par.comm_fortran = -987654;
96 RCP<const Comm<int> > matComm = this->
matrixA_->getComm();
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);
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);
114 function_map::mumps_c(&(mumps_par));
120 mumps_par.icntl[0] = -1;
121 mumps_par.icntl[1] = -1;
122 mumps_par.icntl[2] = -1;
123 mumps_par.icntl[3] = 1;
124 mumps_par.icntl[4] = 0;
125 mumps_par.icntl[5] = 7;
126 mumps_par.icntl[6] = 7;
127 mumps_par.icntl[7] = 7;
128 mumps_par.icntl[8] = 1;
129 mumps_par.icntl[9] = 0;
130 mumps_par.icntl[10] = 0;
131 mumps_par.icntl[11] = 0;
132 mumps_par.icntl[12] = 0;
133 mumps_par.icntl[13] = 20;
134 mumps_par.icntl[17] = 0;
135 mumps_par.icntl[18] = 0;
136 mumps_par.icntl[19] = 0;
137 mumps_par.icntl[20] = 0;
138 mumps_par.icntl[21] = 0;
139 mumps_par.icntl[22] = 0;
140 mumps_par.icntl[23] = 0;
141 mumps_par.icntl[24] = 0;
142 mumps_par.icntl[25] = 0;
143 mumps_par.icntl[27] = 1;
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;
152 template <
class Matrix,
class Vector>
153 MUMPS<Matrix,Vector>::~MUMPS( )
156 typedef FunctionMap<MUMPS,scalar_type> function_map;
158 if(MUMPS_STRUCT ==
true)
165 if (this->rank_ < this->nprocs_) {
166 function_map::mumps_c(&(mumps_par));
171 template<
class Matrix,
class Vector>
177 #ifdef HAVE_AMESOS2_TIMERS
178 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
184 template <
class Matrix,
class Vector>
193 function_map::mumps_c(&(mumps_par));
200 template <
class Matrix,
class Vector>
210 #ifdef HAVE_AMESOS2_TIMERS
211 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
216 function_map::mumps_c(&(mumps_par));
223 template <
class Matrix,
class Vector>
234 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
235 const size_t nrhs = X->getGlobalNumVectors();
237 const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
239 xvals_.resize(val_store_size);
240 bvals_.resize(val_store_size);
242 #ifdef HAVE_AMESOS2_TIMERS
243 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
244 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
248 mumps_type>::do_get(B, bvals_(),
250 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
251 this->rowIndexBase_);
254 mumps_par.nrhs = nrhs;
255 mumps_par.lrhs = mumps_par.n;
260 mumps_par.rhs = bvals_.getRawPtr();
263 #ifdef HAVE_AMESOS2_TIMERS
264 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
267 function_map::mumps_c(&(mumps_par));
271 #ifdef HAVE_AMESOS2_TIMERS
272 Teuchos::TimeMonitor redistTimer2(this->timers_.vecRedistTime_);
278 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
281 MUMPS_MATRIX_LOAD_PREORDERING =
false;
286 template <
class Matrix,
class Vector>
291 return( this->globalNumRows_ == this->globalNumCols_ );
295 template <
class Matrix,
class Vector>
300 using Teuchos::getIntegralValue;
301 using Teuchos::ParameterEntryValidator;
303 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
305 if(parameterList->isParameter(
"ICNTL(1)")){
306 mumps_par.icntl[0] = parameterList->get<
int>(
"ICNTL(1)", -1);
308 if(parameterList->isParameter(
"ICNTL(2)")){
309 mumps_par.icntl[1] = parameterList->get<
int>(
"ICNTL(2)", -1);
311 if(parameterList->isParameter(
"ICNTL(3)")){
312 mumps_par.icntl[2] = parameterList->get<
int>(
"ICNTL(3)", -1);
314 if(parameterList->isParameter(
"ICNTL(4)")){
315 mumps_par.icntl[3] = parameterList->get<
int>(
"ICNTL(4)", 1);
317 if(parameterList->isParameter(
"ICNTL(6)")){
318 mumps_par.icntl[5] = parameterList->get<
int>(
"ICNTL(6)", 0);
320 if(parameterList->isParameter(
"ICNTL(7)")){
321 mumps_par.icntl[6] = parameterList->get<
int>(
"ICNTL(7)", 7);
323 if(parameterList->isParameter(
"ICNTL(9)")){
324 mumps_par.icntl[8] = parameterList->get<
int>(
"ICNTL(9)", 1);
326 if(parameterList->isParameter(
"ICNTL(11)")){
327 mumps_par.icntl[10] = parameterList->get<
int>(
"ICNTL(11)", 0);
329 if(parameterList->isParameter(
"ICNTL(14)")){
330 mumps_par.icntl[13] = parameterList->get<
int>(
"ICNTL(14)", 20);
332 if( parameterList->isParameter(
"IsContiguous") ){
333 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
338 template <
class Matrix,
class Vector>
339 Teuchos::RCP<const Teuchos::ParameterList>
342 using Teuchos::ParameterList;
344 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
346 if( is_null(valid_params) ){
347 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
349 pl->set(
"ICNTL(1)", -1,
"See Manual" );
350 pl->set(
"ICNTL(2)", -1,
"See Manual" );
351 pl->set(
"ICNTL(3)", -1,
"See Manual" );
352 pl->set(
"ICNTL(4)", 1,
"See Manual" );
353 pl->set(
"ICNTL(6)", 0,
"See Manual" );
354 pl->set(
"ICNTL(9)", 1,
"See Manual" );
355 pl->set(
"ICNTL(11)", 0,
"See Manual" );
356 pl->set(
"ICNTL(14)", 20,
"See Manual" );
357 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
366 template <
class Matrix,
class Vector>
372 #ifdef HAVE_AMESOS2_TIMERS
373 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
375 if(MUMPS_MATRIX_LOAD ==
false || (current_phase==NUMFACT && !MUMPS_MATRIX_LOAD_PREORDERING))
378 if( !MUMPS_MATRIX_LOAD && this->root_ ){
379 Kokkos::resize(host_nzvals_view_, this->globalNumNonZeros_);
380 Kokkos::resize(host_rows_view_, this->globalNumNonZeros_);
381 Kokkos::resize(host_col_ptr_view_, this->globalNumRows_ + 1);
384 local_ordinal_type nnz_ret = 0;
386 #ifdef HAVE_AMESOS2_TIMERS
387 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
392 ::do_get(this->matrixA_.ptr(), host_nzvals_view_, host_rows_view_, host_col_ptr_view_, nnz_ret,
393 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
395 this->rowIndexBase_);
398 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
400 "Did not get the expected number of non-zero vals");
412 if (current_phase==PREORDERING){
413 MUMPS_MATRIX_LOAD_PREORDERING =
true;
419 MUMPS_MATRIX_LOAD =
true;
423 template <
class Matrix,
class Vector>
427 if ( !MUMPS_STRUCT ) {
429 mumps_par.n = this->globalNumCols_;
430 mumps_par.nz = this->globalNumNonZeros_;
431 mumps_par.a = (mumps_type*)malloc(mumps_par.nz *
sizeof(mumps_type));
432 mumps_par.irn = (MUMPS_INT*)malloc(mumps_par.nz *
sizeof(MUMPS_INT));
433 mumps_par.jcn = (MUMPS_INT*)malloc(mumps_par.nz *
sizeof(MUMPS_INT));
435 if((mumps_par.a == NULL) || (mumps_par.irn == NULL)
436 || (mumps_par.jcn == NULL))
442 local_ordinal_type tri_count = 0;
443 local_ordinal_type i,j;
444 local_ordinal_type max_local_ordinal = 0;
446 for(i = 0; i < (local_ordinal_type)this->globalNumCols_; i++)
448 for( j = host_col_ptr_view_(i); j < host_col_ptr_view_(i+1)-1; j++)
450 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1;
451 mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1;
452 mumps_par.a[tri_count] = host_nzvals_view_(j);
457 j = host_col_ptr_view_(i+1)-1;
458 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1;
459 mumps_par.irn[tri_count] = (MUMPS_INT)host_rows_view_(j)+1;
460 mumps_par.a[tri_count] = host_nzvals_view_(j);
464 if(host_rows_view_(j) > max_local_ordinal)
466 max_local_ordinal = host_rows_view_(j);
469 TEUCHOS_TEST_FOR_EXCEPTION(std::numeric_limits<MUMPS_INT>::max() <= max_local_ordinal,
471 "Matrix index larger than MUMPS_INT");
476 template<
class Matrix,
class Vector>
478 MUMPS<Matrix,Vector>::MUMPS_ERROR()
const
482 bool Wrong = ((mumps_par.info[0] != 0) || (mumps_par.infog[0] != 0)) && (this->rank_ < this->nprocs_);
484 if (this->rank_==0) {
485 std::cerr <<
"Amesos_Mumps : ERROR" << std::endl;
486 std::cerr <<
"Amesos_Mumps : INFOG(1) = " << mumps_par.infog[0] << std::endl;
487 std::cerr <<
"Amesos_Mumps : INFOG(2) = " << mumps_par.infog[1] << std::endl;
489 if (mumps_par.info[0] != 0 && Wrong) {
490 std::cerr <<
"Amesos_Mumps : On process " << this->matrixA_->getComm()->getRank()
491 <<
", INFO(1) = " << mumps_par.info[0] << std::endl;
492 std::cerr <<
"Amesos_Mumps : On process " << this->matrixA_->getComm()->getRank()
493 <<
", INFO(2) = " << mumps_par.info[1] << std::endl;
499 int WrongInt = Wrong;
500 RCP<const Comm<int> > matComm = this->matrixA_->getComm();
501 Teuchos::broadcast<int,int>(*matComm,0,1,&WrongInt);
502 TEUCHOS_TEST_FOR_EXCEPTION(WrongInt>0,
509 template<
class Matrix,
class Vector>
510 const char* MUMPS<Matrix,Vector>::name =
"MUMPS";
514 #endif // AMESOS2_MUMPS_DEF_HPP
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:648
int numericFactorization_impl()
MUMPS specific numeric factorization.
Definition: Amesos2_MUMPS_def.hpp:202
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:267
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_MUMPS_def.hpp:173
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
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_MUMPS_def.hpp:288
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: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:340
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_MUMPS_def.hpp:297
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:373
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:455