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> 71 template <
class Matrix,
class Vector>
72 MUMPS<Matrix,Vector>::MUMPS(
73 Teuchos::RCP<const Matrix> A,
74 Teuchos::RCP<Vector> X,
75 Teuchos::RCP<const Vector> B )
76 : SolverCore<
Amesos2::MUMPS,Matrix,Vector>(A, X, B)
82 typedef FunctionMap<MUMPS,scalar_type> function_map;
84 MUMPS_MATRIX_LOAD =
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 if(MUMPS_STRUCT ==
true)
164 template<
class Matrix,
class Vector>
170 #ifdef HAVE_AMESOS2_TIMERS 171 Teuchos::TimeMonitor preOrderTimer(this->
timers_.preOrderTime_);
177 template <
class Matrix,
class Vector>
186 function_map::mumps_c(&(mumps_par));
193 template <
class Matrix,
class Vector>
203 #ifdef HAVE_AMESOS2_TIMERS 204 Teuchos::TimeMonitor numFactTimer(this->
timers_.numFactTime_);
207 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG 208 std::cout <<
"MUMPS:: Before numeric factorization" << std::endl;
209 std::cout <<
"nzvals_ : " <<
nzvals_.toString() << std::endl;
210 std::cout <<
"rowind_ : " <<
rowind_.toString() << std::endl;
211 std::cout <<
"colptr_ : " <<
colptr_.toString() << std::endl;
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_ );
247 slu_type>::do_get(B,
bvals_(),as<size_t>(ld_rhs),
252 mumps_par.nrhs = nrhs;
253 mumps_par.lrhs = mumps_par.n;
258 mumps_par.rhs =
bvals_.getRawPtr();
261 #ifdef HAVE_AMESOS2_TIMERS 262 Teuchos::TimeMonitor solveTimer(this->
timers_.solveTime_);
265 function_map::mumps_c(&(mumps_par));
269 #ifdef HAVE_AMESOS2_TIMERS 270 Teuchos::TimeMonitor redistTimer(this->
timers_.vecRedistTime_);
282 template <
class Matrix,
class Vector>
291 template <
class Matrix,
class Vector>
296 using Teuchos::getIntegralValue;
297 using Teuchos::ParameterEntryValidator;
302 if(parameterList->isParameter(
"ICNTL(1)"))
304 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
307 if(parameterList->isParameter(
"ICNTL(2)"))
309 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
312 if(parameterList->isParameter(
"ICNTL(3)"))
314 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
317 if(parameterList->isParameter(
"ICNTL(4)"))
319 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
322 if(parameterList->isParameter(
"ICNTL(6)"))
324 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
327 if(parameterList->isParameter(
"ICNTL(9)"))
329 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
332 if(parameterList->isParameter(
"ICNTL(11)"))
334 mumps_par.icntl[0] = getIntegralValue<local_ordinal_type>(*parameterList,
340 template <
class Matrix,
class Vector>
341 Teuchos::RCP<const Teuchos::ParameterList>
344 using Teuchos::ParameterList;
346 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
348 if( is_null(valid_params) ){
349 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
351 pl->set(
"ICNTL(1)",
"no",
"See Manual" );
352 pl->set(
"ICNTL(2)",
"no",
"See Manual" );
353 pl->set(
"ICNTL(3)",
"no",
"See Manual" );
354 pl->set(
"ICNTL(4)",
"no",
"See Manual" );
355 pl->set(
"ICNTL(6)",
"no",
"See Manual" );
356 pl->set(
"ICNTL(9)",
"no",
"See Manual" );
357 pl->set(
"ICNTL(11)",
"no",
"See Manual" );
366 template <
class Matrix,
class Vector>
372 #ifdef HAVE_AMESOS2_TIMERS 373 Teuchos::TimeMonitor convTimer(this->
timers_.mtxConvTime_);
376 if(MUMPS_MATRIX_LOAD ==
false)
385 local_ordinal_type nnz_ret = 0;
387 #ifdef HAVE_AMESOS2_TIMERS 388 Teuchos::TimeMonitor mtxRedistTimer( this->
timers_.mtxRedistTime_ );
397 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->
globalNumNonZeros_),
399 "Did not get the expected number of non-zero vals");
408 MUMPS_MATRIX_LOAD =
true;
412 template <
class Matrix,
class Vector>
419 mumps_par.a = (magnitude_type*)malloc(mumps_par.nz *
sizeof(magnitude_type));
420 mumps_par.irn = (local_ordinal_type*)malloc(mumps_par.nz *
sizeof(local_ordinal_type));
421 mumps_par.jcn = (local_ordinal_type*)malloc(mumps_par.nz *
sizeof(local_ordinal_type));
423 if((mumps_par.a == NULL) || (mumps_par.irn == NULL)
424 || (mumps_par.jcn == NULL))
430 local_ordinal_type tri_count = 0;
431 local_ordinal_type i,j;
437 mumps_par.jcn[tri_count] = (local_ordinal_type)i+1;
438 mumps_par.irn[tri_count] = (local_ordinal_type)
rowind_[j]+1;
439 mumps_par.a[tri_count] =
nzvals_[j];
447 template<
class Matrix,
class Vector>
451 if(mumps_par.info[0] < 0)
453 TEUCHOS_TEST_FOR_EXCEPTION(
false,
460 template<
class Matrix,
class Vector>
465 #endif // AMESOS2_MUMPS_DEF_HPP int numericFactorization_impl()
Basker specific numeric factorization.
Definition: Amesos2_MUMPS_def.hpp:195
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
bool root_
If true, then this is the root processor.
Definition: Amesos2_SolverCore_decl.hpp:507
Teuchos::Array< local_ordinal_type > rowind_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_MUMPS_decl.hpp:207
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
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:166
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
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
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_MUMPS_def.hpp:284
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:78
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
Teuchos::Array< slu_type > xvals_
Persisting 1D store for X.
Definition: Amesos2_MUMPS_decl.hpp:212
Teuchos::Array< slu_type > bvals_
Persisting 1D store for B.
Definition: Amesos2_MUMPS_decl.hpp:214
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_MUMPS_def.hpp:342
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_MUMPS_def.hpp:293
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
Passes functions to TPL functions based on type.
Definition: Amesos2_FunctionMap.hpp:76
Teuchos::Array< slu_type > nzvals_
Stores the values of the nonzero entries for MUMPS.
Definition: Amesos2_MUMPS_decl.hpp:205
Definition: Amesos2_TypeDecl.hpp:127
Teuchos::Array< local_ordinal_type > colptr_
Stores the row indices of the nonzero entries.
Definition: Amesos2_MUMPS_decl.hpp:209
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
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