20 #ifndef AMESOS2_SHYLUBASKER_DEF_HPP
21 #define AMESOS2_SHYLUBASKER_DEF_HPP
23 #include <Teuchos_Tuple.hpp>
24 #include <Teuchos_ParameterList.hpp>
25 #include <Teuchos_StandardParameterEntryValidators.hpp>
33 template <
class Matrix,
class Vector>
34 ShyLUBasker<Matrix,Vector>::ShyLUBasker(
35 Teuchos::RCP<const Matrix> A,
36 Teuchos::RCP<Vector> X,
37 Teuchos::RCP<const Vector> B )
38 : SolverCore<Amesos2::ShyLUBasker,Matrix,Vector>(A, X, B)
39 , is_contiguous_(true)
47 #if defined(HAVE_AMESOS2_KOKKOS) && defined(KOKKOS_ENABLE_OPENMP)
52 ShyLUbasker = new ::BaskerNS::BaskerTrilinosInterface<local_ordinal_type, shylubasker_dtype, Exe_Space>();
53 ShyLUbasker->Options.no_pivot = BASKER_FALSE;
54 ShyLUbasker->Options.static_delayed_pivot = 0;
55 ShyLUbasker->Options.symmetric = BASKER_FALSE;
56 ShyLUbasker->Options.realloc = BASKER_TRUE;
57 ShyLUbasker->Options.verbose = BASKER_FALSE;
58 ShyLUbasker->Options.prune = BASKER_TRUE;
59 ShyLUbasker->Options.btf_matching = 2;
60 ShyLUbasker->Options.blk_matching = 1;
61 ShyLUbasker->Options.matrix_scaling = 0;
62 ShyLUbasker->Options.min_block_size = 0;
63 ShyLUbasker->Options.amd_dom = BASKER_TRUE;
64 ShyLUbasker->Options.use_metis = BASKER_TRUE;
65 ShyLUbasker->Options.use_nodeNDP = BASKER_TRUE;
66 ShyLUbasker->Options.run_nd_on_leaves = BASKER_TRUE;
67 ShyLUbasker->Options.run_amd_on_leaves = BASKER_FALSE;
68 ShyLUbasker->Options.transpose = BASKER_FALSE;
69 ShyLUbasker->Options.threaded_solve = BASKER_FALSE;
70 ShyLUbasker->Options.replace_zero_pivot = BASKER_TRUE;
71 ShyLUbasker->Options.replace_tiny_pivot = BASKER_FALSE;
72 ShyLUbasker->Options.verbose_matrix_out = BASKER_FALSE;
74 ShyLUbasker->Options.user_fill = (double)BASKER_FILL_USER;
75 ShyLUbasker->Options.use_sequential_diag_facto = BASKER_FALSE;
76 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
77 num_threads = Kokkos::OpenMP::max_hardware_threads();
79 num_threads = Kokkos::OpenMP::impl_max_hardware_threads();
83 TEUCHOS_TEST_FOR_EXCEPTION(1 != 0,
85 "Amesos2_ShyLUBasker Exception: Do not have supported Kokkos node type (OpenMP) enabled for ShyLUBasker");
90 template <
class Matrix,
class Vector>
91 ShyLUBasker<Matrix,Vector>::~ShyLUBasker( )
94 #if defined(HAVE_AMESOS2_KOKKOS) && defined(KOKKOS_ENABLE_OPENMP)
95 ShyLUbasker->Finalize();
100 template <
class Matrix,
class Vector>
103 return (this->root_ && (this->matrixA_->getComm()->getSize() == 1) && is_contiguous_);
106 template<
class Matrix,
class Vector>
112 #ifdef HAVE_AMESOS2_TIMERS
113 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
120 template <
class Matrix,
class Vector>
128 ShyLUbasker->SetThreads(num_threads);
136 if ( single_proc_optimization() ) {
138 host_ordinal_type_array sp_rowptr;
139 host_ordinal_type_array sp_colind;
141 this->matrixA_->returnRowPtr_kokkos_view(sp_rowptr);
142 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr.data() ==
nullptr,
143 std::runtime_error,
"Amesos2 Runtime Error: sp_rowptr returned null ");
144 this->matrixA_->returnColInd_kokkos_view(sp_colind);
145 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind.data() ==
nullptr,
146 std::runtime_error,
"Amesos2 Runtime Error: sp_colind returned null ");
148 host_value_type_array hsp_values;
149 this->matrixA_->returnValues_kokkos_view(hsp_values);
150 shylubasker_dtype * sp_values = function_map::convert_scalar(hsp_values.data());
152 TEUCHOS_TEST_FOR_EXCEPTION(sp_values ==
nullptr,
153 std::runtime_error,
"Amesos2 Runtime Error: sp_values returned null ");
156 info = ShyLUbasker->Symbolic(this->globalNumRows_,
157 this->globalNumCols_,
158 this->globalNumNonZeros_,
164 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
165 std::runtime_error,
"Error in ShyLUBasker Symbolic");
170 shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
171 info = ShyLUbasker->Symbolic(this->globalNumRows_,
172 this->globalNumCols_,
173 this->globalNumNonZeros_,
178 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
179 std::runtime_error,
"Error in ShyLUBasker Symbolic");
185 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
190 template <
class Matrix,
class Vector>
199 #ifdef HAVE_AMESOS2_TIMERS
200 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
208 if ( single_proc_optimization() ) {
210 host_ordinal_type_array sp_rowptr;
211 host_ordinal_type_array sp_colind;
212 this->matrixA_->returnRowPtr_kokkos_view(sp_rowptr);
213 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr.data() ==
nullptr,
214 std::runtime_error,
"Amesos2 Runtime Error: sp_rowptr returned null ");
215 this->matrixA_->returnColInd_kokkos_view(sp_colind);
216 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind.data() ==
nullptr,
217 std::runtime_error,
"Amesos2 Runtime Error: sp_colind returned null ");
219 host_value_type_array hsp_values;
220 this->matrixA_->returnValues_kokkos_view(hsp_values);
221 shylubasker_dtype * sp_values = function_map::convert_scalar(hsp_values.data());
224 TEUCHOS_TEST_FOR_EXCEPTION(sp_values ==
nullptr,
225 std::runtime_error,
"Amesos2 Runtime Error: sp_values returned null ");
228 info = ShyLUbasker->Factor( this->globalNumRows_,
229 this->globalNumCols_,
230 this->globalNumNonZeros_,
238 shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
239 info = ShyLUbasker->Factor(this->globalNumRows_,
240 this->globalNumCols_,
241 this->globalNumNonZeros_,
249 local_ordinal_type blnnz = local_ordinal_type(0);
250 local_ordinal_type bunnz = local_ordinal_type(0);
251 ShyLUbasker->GetLnnz(blnnz);
252 ShyLUbasker->GetUnnz(bunnz);
256 this->setNnzLU( as<size_t>( blnnz + bunnz ) );
262 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
265 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
266 std::runtime_error,
" ShyLUBasker::numericFactorization failed.");
272 template <
class Matrix,
class Vector>
282 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
283 const size_t nrhs = X->getGlobalNumVectors();
285 const bool ShyluBaskerTransposeRequest = this->control_.useTranspose_;
286 const bool initialize_data =
true;
287 const bool do_not_initialize_data =
false;
288 bool use_gather = use_gather_;
289 use_gather = (use_gather && this->matrixA_->getComm()->getSize() > 1);
290 use_gather = (use_gather && (std::is_same<scalar_type, float>::value || std::is_same<scalar_type, double>::value));
292 #ifdef HAVE_AMESOS2_TIMERS
293 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
295 if ( single_proc_optimization() && nrhs == 1 ) {
298 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
299 host_solve_array_t>::do_get(initialize_data, B, bValues_, as<size_t>(ld_rhs));
301 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
302 host_solve_array_t>::do_get(do_not_initialize_data, X, xValues_, as<size_t>(ld_rhs));
307 int rval = B->gather(bValues_, this->perm_g2l, this->recvCountRows, this->recvDisplRows,
308 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
310 X->gather(xValues_, this->perm_g2l, this->recvCountRows, this->recvDisplRows,
311 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
317 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
318 host_solve_array_t>::do_get(initialize_data, B, bValues_,
320 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
321 this->rowIndexBase_);
324 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
325 host_solve_array_t>::do_get(do_not_initialize_data, X, xValues_,
327 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
328 this->rowIndexBase_);
334 #ifdef HAVE_AMESOS2_TIMERS
335 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
338 shylubasker_dtype * pxValues = function_map::convert_scalar(xValues_.data());
339 shylubasker_dtype * pbValues = function_map::convert_scalar(bValues_.data());
340 if (!ShyluBaskerTransposeRequest)
341 ierr = ShyLUbasker->Solve(nrhs, pbValues, pxValues);
343 ierr = ShyLUbasker->Solve(nrhs, pbValues, pxValues,
true);
346 Teuchos::broadcast(*(this->getComm()), 0, &ierr);
348 TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
350 "Encountered zero diag element at: " << ierr);
351 TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
353 "Could not alloc needed working memory for solve" );
355 #ifdef HAVE_AMESOS2_TIMERS
356 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
359 int rval = X->scatter(xValues_, this->perm_g2l, this->recvCountRows, this->recvDisplRows,
360 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
361 if (rval != 0) use_gather =
false;
364 Util::put_1d_data_helper_kokkos_view<
367 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
374 template <
class Matrix,
class Vector>
379 return( this->globalNumRows_ == this->globalNumCols_ );
383 template <
class Matrix,
class Vector>
388 using Teuchos::getIntegralValue;
389 using Teuchos::ParameterEntryValidator;
391 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
393 if(parameterList->isParameter(
"IsContiguous"))
395 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
398 if(parameterList->isParameter(
"UseCustomGather"))
400 use_gather_ = parameterList->get<
bool>(
"UseCustomGather");
403 if(parameterList->isParameter(
"num_threads"))
405 num_threads = parameterList->get<
int>(
"num_threads");
407 if(parameterList->isParameter(
"pivot"))
409 ShyLUbasker->Options.no_pivot = (!parameterList->get<
bool>(
"pivot"));
411 if(parameterList->isParameter(
"delayed pivot"))
413 ShyLUbasker->Options.static_delayed_pivot = (parameterList->get<
int>(
"delayed pivot"));
415 if(parameterList->isParameter(
"pivot_tol"))
417 ShyLUbasker->Options.pivot_tol = parameterList->get<
double>(
"pivot_tol");
419 if(parameterList->isParameter(
"symmetric"))
421 ShyLUbasker->Options.symmetric = parameterList->get<
bool>(
"symmetric");
423 if(parameterList->isParameter(
"realloc"))
425 ShyLUbasker->Options.realloc = parameterList->get<
bool>(
"realloc");
427 if(parameterList->isParameter(
"verbose"))
429 ShyLUbasker->Options.verbose = parameterList->get<
bool>(
"verbose");
431 if(parameterList->isParameter(
"verbose_matrix"))
433 ShyLUbasker->Options.verbose_matrix_out = parameterList->get<
bool>(
"verbose_matrix");
435 if(parameterList->isParameter(
"btf"))
437 ShyLUbasker->Options.btf = parameterList->get<
bool>(
"btf");
439 if(parameterList->isParameter(
"use_metis"))
441 ShyLUbasker->Options.use_metis = parameterList->get<
bool>(
"use_metis");
443 if(parameterList->isParameter(
"use_nodeNDP"))
445 ShyLUbasker->Options.use_nodeNDP = parameterList->get<
bool>(
"use_nodeNDP");
447 if(parameterList->isParameter(
"run_nd_on_leaves"))
449 ShyLUbasker->Options.run_nd_on_leaves = parameterList->get<
bool>(
"run_nd_on_leaves");
451 if(parameterList->isParameter(
"run_amd_on_leaves"))
453 ShyLUbasker->Options.run_amd_on_leaves = parameterList->get<
bool>(
"run_amd_on_leaves");
455 if(parameterList->isParameter(
"amd_on_blocks"))
457 ShyLUbasker->Options.amd_dom = parameterList->get<
bool>(
"amd_on_blocks");
459 if(parameterList->isParameter(
"transpose"))
462 const auto transpose = parameterList->get<
bool>(
"transpose");
464 this->control_.useTranspose_ =
true;
466 if(parameterList->isParameter(
"threaded_solve"))
468 ShyLUbasker->Options.threaded_solve = parameterList->get<
bool>(
"threaded_solve");
470 if(parameterList->isParameter(
"use_sequential_diag_facto"))
472 ShyLUbasker->Options.use_sequential_diag_facto = parameterList->get<
bool>(
"use_sequential_diag_facto");
474 if(parameterList->isParameter(
"user_fill"))
476 ShyLUbasker->Options.user_fill = parameterList->get<
double>(
"user_fill");
478 if(parameterList->isParameter(
"prune"))
480 ShyLUbasker->Options.prune = parameterList->get<
bool>(
"prune");
482 if(parameterList->isParameter(
"replace_zero_pivot"))
484 ShyLUbasker->Options.replace_zero_pivot = parameterList->get<
bool>(
"replace_zero_pivot");
486 if(parameterList->isParameter(
"replace_tiny_pivot"))
488 ShyLUbasker->Options.replace_tiny_pivot = parameterList->get<
bool>(
"replace_tiny_pivot");
490 if(parameterList->isParameter(
"btf_matching"))
492 ShyLUbasker->Options.btf_matching = parameterList->get<
int>(
"btf_matching");
493 if (ShyLUbasker->Options.btf_matching == 1 || ShyLUbasker->Options.btf_matching == 2) {
494 ShyLUbasker->Options.matching =
true;
496 ShyLUbasker->Options.matching =
false;
499 if(parameterList->isParameter(
"blk_matching"))
501 ShyLUbasker->Options.blk_matching = parameterList->get<
int>(
"blk_matching");
503 if(parameterList->isParameter(
"matrix_scaling"))
505 ShyLUbasker->Options.matrix_scaling = parameterList->get<
int>(
"matrix_scaling");
507 if(parameterList->isParameter(
"min_block_size"))
509 ShyLUbasker->Options.min_block_size = parameterList->get<
int>(
"min_block_size");
513 template <
class Matrix,
class Vector>
514 Teuchos::RCP<const Teuchos::ParameterList>
517 using Teuchos::ParameterList;
519 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
521 if( is_null(valid_params) )
523 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
524 pl->set(
"IsContiguous",
true,
525 "Are GIDs contiguous");
526 pl->set(
"UseCustomGather",
true,
527 "Use Matrix-gather routine");
528 pl->set(
"num_threads", 1,
529 "Number of threads");
530 pl->set(
"pivot",
false,
532 pl->set(
"delayed pivot", 0,
533 "Apply static delayed pivot on a big block");
534 pl->set(
"pivot_tol", .0001,
535 "Tolerance before pivot, currently not used");
536 pl->set(
"symmetric",
false,
537 "Should Symbolic assume symmetric nonzero pattern");
538 pl->set(
"realloc" ,
false,
539 "Should realloc space if not enough");
540 pl->set(
"verbose",
false,
541 "Information about factoring");
542 pl->set(
"verbose_matrix",
false,
543 "Give Permuted Matrices");
546 pl->set(
"prune",
false,
547 "Use prune on BTF blocks (Not Supported)");
548 pl->set(
"btf_matching", 2,
549 "Matching option for BTF: 0 = none, 1 = Basker, 2 = Trilinos (default), (3 = MC64 if enabled)");
550 pl->set(
"blk_matching", 1,
551 "Matching optioon for block: 0 = none, 1 or anything else = Basker (default), (2 = MC64 if enabled)");
552 pl->set(
"matrix_scaling", 0,
553 "Use matrix scaling to biig A BTF block: 0 = no-scaling, 1 = symmetric diagonal scaling, 2 = row-max, and then col-max scaling");
554 pl->set(
"min_block_size", 0,
555 "Size of the minimum diagonal blocks");
556 pl->set(
"replace_zero_pivot",
true,
557 "Replace zero pivots during the numerical factorization");
558 pl->set(
"replace_tiny_pivot",
false,
559 "Replace tiny pivots during the numerical factorization");
560 pl->set(
"use_metis",
true,
562 pl->set(
"use_nodeNDP",
true,
563 "Use nodeND to compute ND partition");
564 pl->set(
"run_nd_on_leaves",
false,
565 "Run ND on the final leaf-nodes for ND factorization");
566 pl->set(
"run_amd_on_leaves",
false,
567 "Run AMD on the final leaf-nodes for ND factorization");
568 pl->set(
"amd_on_blocks",
true,
569 "Run AMD on each diagonal blocks");
570 pl->set(
"transpose",
false,
571 "Solve the transpose A");
572 pl->set(
"threaded_solve",
false,
573 "Use threads for forward/backward solves");
574 pl->set(
"use_sequential_diag_facto",
false,
575 "Use sequential algorithm to factor each diagonal block");
576 pl->set(
"user_fill", (
double)BASKER_FILL_USER,
577 "User-provided padding for the fill ratio");
584 template <
class Matrix,
class Vector>
589 if(current_phase == SOLVE || current_phase == PREORDERING )
return(
false );
591 #ifdef HAVE_AMESOS2_TIMERS
592 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
597 if ( single_proc_optimization() ) {
605 if( this->root_ && current_phase == SYMBFACT )
607 Kokkos::resize(nzvals_view_, this->globalNumNonZeros_);
608 Kokkos::resize(rowind_view_, this->globalNumNonZeros_);
609 Kokkos::resize(colptr_view_, this->globalNumCols_ + 1);
612 local_ordinal_type nnz_ret = -1;
613 bool use_gather = use_gather_;
614 use_gather = (use_gather && this->matrixA_->getComm()->getSize() > 1);
615 use_gather = (use_gather && (std::is_same<scalar_type, float>::value || std::is_same<scalar_type, double>::value));
617 #ifdef HAVE_AMESOS2_TIMERS
618 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
621 bool column_major =
true;
622 if (!is_contiguous_) {
623 auto contig_mat = this->matrixA_->reindex(this->contig_rowmap_, this->contig_colmap_, current_phase);
624 nnz_ret = contig_mat->gather(nzvals_view_, rowind_view_, colptr_view_, this->perm_g2l, this->recvCountRows, this->recvDisplRows, this->recvCounts, this->recvDispls,
625 this->transpose_map, this->nzvals_t, column_major, current_phase);
627 nnz_ret = this->matrixA_->gather(nzvals_view_, rowind_view_, colptr_view_, this->perm_g2l, this->recvCountRows, this->recvDisplRows, this->recvCounts, this->recvDispls,
628 this->transpose_map, this->nzvals_t, column_major, current_phase);
632 if (nnz_ret < 0) use_gather =
false;
637 ::do_get(this->matrixA_.ptr(), nzvals_view_, rowind_view_, colptr_view_, nnz_ret,
638 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
640 this->rowIndexBase_);
645 if (use_gather || this->root_) {
646 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
648 "Amesos2_ShyLUBasker loadA_impl: Did not get the expected number of non-zero vals("
649 +std::to_string(nnz_ret)+
" vs "+std::to_string(this->globalNumNonZeros_)+
")");
656 template<
class Matrix,
class Vector>
662 #endif // AMESOS2_SHYLUBASKER_DEF_HPP
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:618
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
ShyLUBasker specific solve.
Definition: Amesos2_ShyLUBasker_def.hpp:274
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:31
Amesos2 interface to the Baker package.
Definition: Amesos2_ShyLUBasker_decl.hpp:42
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_ShyLUBasker_def.hpp:515
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_ShyLUBasker_def.hpp:108
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_ShyLUBasker_def.hpp:376
Amesos2 ShyLUBasker declarations.
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_ShyLUBasker_def.hpp:586
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:42
int numericFactorization_impl()
ShyLUBasker specific numeric factorization.
Definition: Amesos2_ShyLUBasker_def.hpp:192
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:142
bool single_proc_optimization() const
can we optimize size_type and ordinal_type for straight pass through, also check that is_contiguous_ ...
Definition: Amesos2_ShyLUBasker_def.hpp:102