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)
83 typedef FunctionMap<MUMPS,scalar_type> function_map;
85 MUMPS_MATRIX_LOAD =
false;
90 using Teuchos::MpiComm;
93 using Teuchos::rcp_dynamic_cast;
96 mumps_par.comm_fortran = -987654;
97 RCP<const Comm<int> > matComm = this->
matrixA_->getComm();
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);
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);
115 function_map::mumps_c(&(mumps_par));
121 mumps_par.icntl[0] = -1;
122 mumps_par.icntl[1] = -1;
123 mumps_par.icntl[2] = -1;
124 mumps_par.icntl[3] = 1;
125 mumps_par.icntl[4] = 0;
126 mumps_par.icntl[5] = 7;
127 mumps_par.icntl[6] = 7;
128 mumps_par.icntl[7] = 7;
129 mumps_par.icntl[8] = 1;
130 mumps_par.icntl[9] = 0;
131 mumps_par.icntl[10] = 0;
132 mumps_par.icntl[11] = 0;
133 mumps_par.icntl[12] = 0;
134 mumps_par.icntl[13] = 20;
135 mumps_par.icntl[17] = 0;
136 mumps_par.icntl[18] = 0;
137 mumps_par.icntl[19] = 0;
138 mumps_par.icntl[20] = 0;
139 mumps_par.icntl[21] = 0;
140 mumps_par.icntl[22] = 0;
141 mumps_par.icntl[23] = 0;
142 mumps_par.icntl[24] = 0;
143 mumps_par.icntl[25] = 0;
144 mumps_par.icntl[27] = 1;
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;
153 template <
class Matrix,
class Vector>
154 MUMPS<Matrix,Vector>::~MUMPS( )
157 if(MUMPS_STRUCT ==
true)
165 template<
class Matrix,
class Vector>
171 #ifdef HAVE_AMESOS2_TIMERS 172 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
178 template <
class Matrix,
class Vector>
187 function_map::mumps_c(&(mumps_par));
194 template <
class Matrix,
class Vector>
204 #ifdef HAVE_AMESOS2_TIMERS 205 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
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;
217 function_map::mumps_c(&(mumps_par));
224 template <
class Matrix,
class Vector>
235 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
236 const size_t nrhs = X->getGlobalNumVectors();
238 const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
240 xvals_.resize(val_store_size);
241 bvals_.resize(val_store_size);
243 #ifdef HAVE_AMESOS2_TIMERS 244 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
245 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
248 slu_type>::do_get(B, bvals_(),as<size_t>(ld_rhs),
253 mumps_par.nrhs = nrhs;
254 mumps_par.lrhs = mumps_par.n;
259 mumps_par.rhs = bvals_.getRawPtr();
262 #ifdef HAVE_AMESOS2_TIMERS 263 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
266 function_map::mumps_c(&(mumps_par));
270 #ifdef HAVE_AMESOS2_TIMERS 271 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
283 template <
class Matrix,
class Vector>
288 return( this->globalNumRows_ == this->globalNumCols_ );
292 template <
class Matrix,
class Vector>
297 using Teuchos::getIntegralValue;
298 using Teuchos::ParameterEntryValidator;
300 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
303 if(parameterList->isParameter(
"ICNTL(1)"))
305 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
308 if(parameterList->isParameter(
"ICNTL(2)"))
310 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
313 if(parameterList->isParameter(
"ICNTL(3)"))
315 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
318 if(parameterList->isParameter(
"ICNTL(4)"))
320 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
323 if(parameterList->isParameter(
"ICNTL(6)"))
325 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
328 if(parameterList->isParameter(
"ICNTL(9)"))
330 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
333 if(parameterList->isParameter(
"ICNTL(11)"))
335 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
341 template <
class Matrix,
class Vector>
342 Teuchos::RCP<const Teuchos::ParameterList>
345 using Teuchos::ParameterList;
347 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
349 if( is_null(valid_params) ){
350 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
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" );
367 template <
class Matrix,
class Vector>
373 #ifdef HAVE_AMESOS2_TIMERS 374 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
377 if(MUMPS_MATRIX_LOAD ==
false)
381 nzvals_.resize(this->globalNumNonZeros_);
382 rowind_.resize(this->globalNumNonZeros_);
383 colptr_.resize(this->globalNumCols_ + 1);
386 local_ordinal_type nnz_ret = 0;
388 #ifdef HAVE_AMESOS2_TIMERS 389 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
394 ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
398 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
400 "Did not get the expected number of non-zero vals");
409 MUMPS_MATRIX_LOAD =
true;
413 template <
class Matrix,
class Vector>
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));
424 if((mumps_par.a == NULL) || (mumps_par.irn == NULL)
425 || (mumps_par.jcn == NULL))
431 local_ordinal_type tri_count = 0;
432 local_ordinal_type i,j;
433 local_ordinal_type max_local_ordinal = 0;
435 for(i = 0; i < (local_ordinal_type)this->globalNumCols_; i++)
437 for( j = colptr_[i]; j < colptr_[i+1]-1; j++)
439 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1;
440 mumps_par.irn[tri_count] = (MUMPS_INT)rowind_[j]+1;
441 mumps_par.a[tri_count] = nzvals_[j];
447 mumps_par.jcn[tri_count] = (MUMPS_INT)i+1;
448 mumps_par.irn[tri_count] = (MUMPS_INT)rowind_[j]+1;
449 mumps_par.a[tri_count] = nzvals_[j];
453 if(rowind_[j] > max_local_ordinal)
455 max_local_ordinal = rowind_[j];
458 TEUCHOS_TEST_FOR_EXCEPTION(std::numeric_limits<MUMPS_INT>::max() <= max_local_ordinal,
460 "Matrix index larger than MUMPS_INT");
465 template<
class Matrix,
class Vector>
467 MUMPS<Matrix,Vector>::MUMPS_ERROR()
const 469 if(mumps_par.info[0] < 0)
471 TEUCHOS_TEST_FOR_EXCEPTION(
false,
478 template<
class Matrix,
class Vector>
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 > ¶meterList)
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