52 #ifndef AMESOS2_KLU2_DEF_HPP
53 #define AMESOS2_KLU2_DEF_HPP
55 #include <Teuchos_Tuple.hpp>
56 #include <Teuchos_ParameterList.hpp>
57 #include <Teuchos_StandardParameterEntryValidators.hpp>
65 template <
class Matrix,
class Vector>
67 Teuchos::RCP<const Matrix> A,
68 Teuchos::RCP<Vector> X,
69 Teuchos::RCP<const Vector> B )
75 , is_contiguous_(true)
77 ::KLU2::klu_defaults<slu_type, local_ordinal_type> (&(data_.common_)) ;
78 data_.symbolic_ = NULL;
79 data_.numeric_ = NULL;
86 template <
class Matrix,
class Vector>
94 if (data_.symbolic_ != NULL)
95 ::KLU2::klu_free_symbolic<slu_type, local_ordinal_type>
96 (&(data_.symbolic_), &(data_.common_)) ;
97 if (data_.numeric_ != NULL)
98 ::KLU2::klu_free_numeric<slu_type, local_ordinal_type>
99 (&(data_.numeric_), &(data_.common_)) ;
112 template <
class Matrix,
class Vector>
115 return (this->root_ && (this->matrixA_->getComm()->getSize() == 1) && is_contiguous_);
118 template<
class Matrix,
class Vector>
124 #ifdef HAVE_AMESOS2_TIMERS
125 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
132 template <
class Matrix,
class Vector>
136 if (data_.symbolic_ != NULL) {
137 ::KLU2::klu_free_symbolic<slu_type, local_ordinal_type>
138 (&(data_.symbolic_), &(data_.common_)) ;
141 if ( single_proc_optimization() ) {
143 auto sp_rowptr = this->matrixA_->returnRowPtr();
144 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr ==
nullptr,
145 std::runtime_error,
"Amesos2 Runtime Error: sp_rowptr returned null ");
146 auto sp_colind = this->matrixA_->returnColInd();
147 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind ==
nullptr,
148 std::runtime_error,
"Amesos2 Runtime Error: sp_colind returned null ");
149 #ifndef HAVE_TEUCHOS_COMPLEX
150 auto sp_values = this->matrixA_->returnValues();
154 using complex_type =
typename Util::getStdCplxType< magnitude_type, typename matrix_adapter_type::spmtx_vals_t >::type;
155 complex_type * sp_values =
nullptr;
156 sp_values =
reinterpret_cast< complex_type *
> ( this->matrixA_->returnValues() );
158 TEUCHOS_TEST_FOR_EXCEPTION(sp_values ==
nullptr,
159 std::runtime_error,
"Amesos2 Runtime Error: sp_values returned null ");
164 Teuchos::Array< local_ordinal_type > sp_rowptr_with_common_type( this->globalNumRows_ + 1 );
165 for ( global_size_type i = 0; i < this->globalNumRows_+1; ++i ) {
167 sp_rowptr_with_common_type[i] = sp_rowptr[i];
170 data_.symbolic_ = ::KLU2::klu_analyze<slu_type, local_ordinal_type>
171 ((local_ordinal_type)this->globalNumCols_, sp_rowptr_with_common_type.getRawPtr(),
172 sp_colind, &(data_.common_)) ;
176 data_.symbolic_ = ::KLU2::klu_analyze<slu_type, local_ordinal_type>
177 ((local_ordinal_type)this->globalNumCols_, colptr_.getRawPtr(),
178 rowind_.getRawPtr(), &(data_.common_)) ;
186 template <
class Matrix,
class Vector>
200 #ifdef HAVE_AMESOS2_TIMERS
201 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
204 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
205 std::cout <<
"KLU2:: Before numeric factorization" << std::endl;
206 std::cout <<
"nzvals_ : " << nzvals_.toString() << std::endl;
207 std::cout <<
"rowind_ : " << rowind_.toString() << std::endl;
208 std::cout <<
"colptr_ : " << colptr_.toString() << std::endl;
211 if (data_.numeric_ != NULL) {
212 ::KLU2::klu_free_numeric<slu_type, local_ordinal_type>
213 (&(data_.numeric_), &(data_.common_)) ;
216 if ( single_proc_optimization() ) {
218 auto sp_rowptr = this->matrixA_->returnRowPtr();
219 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr ==
nullptr,
220 std::runtime_error,
"Amesos2 Runtime Error: sp_rowptr returned null ");
221 auto sp_colind = this->matrixA_->returnColInd();
222 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind ==
nullptr,
223 std::runtime_error,
"Amesos2 Runtime Error: sp_colind returned null ");
224 #ifndef HAVE_TEUCHOS_COMPLEX
225 auto sp_values = this->matrixA_->returnValues();
229 using complex_type =
typename Util::getStdCplxType< magnitude_type, typename matrix_adapter_type::spmtx_vals_t >::type;
230 complex_type * sp_values =
nullptr;
231 sp_values =
reinterpret_cast< complex_type *
> ( this->matrixA_->returnValues() );
233 TEUCHOS_TEST_FOR_EXCEPTION(sp_values ==
nullptr,
234 std::runtime_error,
"Amesos2 Runtime Error: sp_values returned null ");
239 Teuchos::Array< local_ordinal_type > sp_rowptr_with_common_type( this->globalNumRows_ + 1 );
240 for ( global_size_type i = 0; i < this->globalNumRows_+1; ++i ) {
242 sp_rowptr_with_common_type[i] = sp_rowptr[i];
245 data_.numeric_ = ::KLU2::klu_factor<slu_type, local_ordinal_type>
246 (sp_rowptr_with_common_type.getRawPtr(), sp_colind, sp_values,
247 data_.symbolic_, &(data_.common_)) ;
250 data_.numeric_ = ::KLU2::klu_factor<slu_type, local_ordinal_type>
251 (colptr_.getRawPtr(), rowind_.getRawPtr(), nzvals_.getRawPtr(),
252 data_.symbolic_, &(data_.common_)) ;
260 if(data_.numeric_ ==
nullptr) {
267 this->setNnzLU( as<size_t>((data_.numeric_)->lnz) + as<size_t>((data_.numeric_)->unz) );
274 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
276 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::runtime_error,
277 "KLU2 numeric factorization failed");
295 template <
class Matrix,
class Vector>
304 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
305 const size_t nrhs = X->getGlobalNumVectors();
307 if ( single_proc_optimization() && nrhs == 1 ) {
308 #ifdef HAVE_AMESOS2_TIMERS
309 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
312 #ifndef HAVE_TEUCHOS_COMPLEX
318 using complex_type =
typename Util::getStdCplxType< magnitude_type, typename matrix_adapter_type::spmtx_vals_t >::type;
322 TEUCHOS_TEST_FOR_EXCEPTION(b_vector ==
nullptr,
323 std::runtime_error,
"Amesos2 Runtime Error: b_vector returned null ");
325 TEUCHOS_TEST_FOR_EXCEPTION(x_vector ==
nullptr,
326 std::runtime_error,
"Amesos2 Runtime Error: x_vector returned null ");
333 ::KLU2::klu_tsolve2<slu_type, local_ordinal_type>
334 (data_.symbolic_, data_.numeric_,
335 (local_ordinal_type)this->globalNumCols_,
336 (local_ordinal_type)nrhs,
337 b_vector, x_vector, &(data_.common_)) ;
340 ::KLU2::klu_solve2<slu_type, local_ordinal_type>
341 (data_.symbolic_, data_.numeric_,
342 (local_ordinal_type)this->globalNumCols_,
343 (local_ordinal_type)nrhs,
344 b_vector, x_vector, &(data_.common_)) ;
354 const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
355 Teuchos::Array<slu_type> bValues(val_store_size);
358 #ifdef HAVE_AMESOS2_TIMERS
359 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
360 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
362 if ( is_contiguous_ ==
true ) {
364 slu_type>::do_get(B, bValues(),
366 ROOTED, this->rowIndexBase_);
370 slu_type>::do_get(B, bValues(),
372 CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
379 #ifdef HAVE_AMESOS2_TIMERS
380 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
387 if ( single_proc_optimization() ) {
388 ::KLU2::klu_tsolve<slu_type, local_ordinal_type>
389 (data_.symbolic_, data_.numeric_,
390 (local_ordinal_type)this->globalNumCols_,
391 (local_ordinal_type)nrhs,
392 bValues.getRawPtr(), &(data_.common_)) ;
395 ::KLU2::klu_solve<slu_type, local_ordinal_type>
396 (data_.symbolic_, data_.numeric_,
397 (local_ordinal_type)this->globalNumCols_,
398 (local_ordinal_type)nrhs,
399 bValues.getRawPtr(), &(data_.common_)) ;
407 if ( single_proc_optimization() ) {
408 ::KLU2::klu_solve<slu_type, local_ordinal_type>
409 (data_.symbolic_, data_.numeric_,
410 (local_ordinal_type)this->globalNumCols_,
411 (local_ordinal_type)nrhs,
412 bValues.getRawPtr(), &(data_.common_)) ;
415 ::KLU2::klu_tsolve<slu_type, local_ordinal_type>
416 (data_.symbolic_, data_.numeric_,
417 (local_ordinal_type)this->globalNumCols_,
418 (local_ordinal_type)nrhs,
419 bValues.getRawPtr(), &(data_.common_)) ;
443 #ifdef HAVE_AMESOS2_TIMERS
444 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
447 if ( is_contiguous_ ==
true ) {
451 ROOTED, this->rowIndexBase_);
457 CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
466 template <
class Matrix,
class Vector>
473 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
477 template <
class Matrix,
class Vector>
482 using Teuchos::getIntegralValue;
483 using Teuchos::ParameterEntryValidator;
485 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
487 transFlag_ = this->control_.useTranspose_ ? 1: 0;
489 if( parameterList->isParameter(
"Trans") ){
490 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"Trans").validator();
491 parameterList->getEntry(
"Trans").setValidator(trans_validator);
493 transFlag_ = getIntegralValue<int>(*parameterList,
"Trans");
496 if( parameterList->isParameter(
"IsContiguous") ){
497 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
502 template <
class Matrix,
class Vector>
503 Teuchos::RCP<const Teuchos::ParameterList>
507 using Teuchos::tuple;
508 using Teuchos::ParameterList;
509 using Teuchos::setStringToIntegralParameter;
511 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
513 if( is_null(valid_params) )
515 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
517 pl->set(
"Equil",
true,
"Whether to equilibrate the system before solve, does nothing now");
518 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
520 setStringToIntegralParameter<int>(
"Trans",
"NOTRANS",
521 "Solve for the transpose system or not",
522 tuple<string>(
"NOTRANS",
"TRANS",
"CONJ"),
523 tuple<string>(
"Solve with transpose",
524 "Do not solve with transpose",
525 "Solve with the conjugate transpose"),
535 template <
class Matrix,
class Vector>
541 if(current_phase == SOLVE)
return(
false);
543 if ( single_proc_optimization() ) {
549 #ifdef HAVE_AMESOS2_TIMERS
550 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
555 nzvals_.resize(this->globalNumNonZeros_);
556 rowind_.resize(this->globalNumNonZeros_);
557 colptr_.resize(this->globalNumCols_ + 1);
560 local_ordinal_type nnz_ret = 0;
562 #ifdef HAVE_AMESOS2_TIMERS
563 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
566 if ( is_contiguous_ ==
true ) {
569 ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
570 nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_);
575 ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
576 nnz_ret, CONTIGUOUS_AND_ROOTED, ARBITRARY, this->rowIndexBase_);
582 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
584 "Did not get the expected number of non-zero vals");
593 template<
class Matrix,
class Vector>
599 #endif // AMESOS2_KLU2_DEF_HPP
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
KLU2(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_KLU2_def.hpp:66
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
KLU2 specific solve.
Definition: Amesos2_KLU2_def.hpp:297
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
Amesos2 KLU2 declarations.
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_KLU2_def.hpp:537
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:266
~KLU2()
Destructor.
Definition: Amesos2_KLU2_def.hpp:87
Helper struct for getting pointers to the MV data - only used when number of vectors = 1 and single M...
Definition: Amesos2_MultiVecAdapter_decl.hpp:218
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_KLU2_def.hpp:479
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using KLU2.
Definition: Amesos2_KLU2_def.hpp:134
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:654
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_KLU2_def.hpp:120
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_KLU2_def.hpp:468
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
int numericFactorization_impl()
KLU2 specific numeric factorization.
Definition: Amesos2_KLU2_def.hpp:188
Amesos2 interface to the KLU2 package.
Definition: Amesos2_KLU2_decl.hpp:72
bool single_proc_optimization() const
can we optimize size_type and ordinal_type for straight pass through, also check that is_contiguous_ ...
Definition: Amesos2_KLU2_def.hpp:114
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_KLU2_def.hpp:504
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:322
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176