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.replace_zero_pivot = BASKER_TRUE;
70 ShyLUbasker->Options.replace_tiny_pivot = BASKER_FALSE;
71 ShyLUbasker->Options.verbose_matrix_out = BASKER_FALSE;
73 ShyLUbasker->Options.user_fill = (double)BASKER_FILL_USER;
74 ShyLUbasker->Options.use_sequential_diag_facto = BASKER_FALSE;
75 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
76 num_threads = Kokkos::OpenMP::max_hardware_threads();
78 num_threads = Kokkos::OpenMP::impl_max_hardware_threads();
82 TEUCHOS_TEST_FOR_EXCEPTION(1 != 0,
84 "Amesos2_ShyLUBasker Exception: Do not have supported Kokkos node type (OpenMP) enabled for ShyLUBasker");
89 template <
class Matrix,
class Vector>
90 ShyLUBasker<Matrix,Vector>::~ShyLUBasker( )
93 #if defined(HAVE_AMESOS2_KOKKOS) && defined(KOKKOS_ENABLE_OPENMP)
94 ShyLUbasker->Finalize();
99 template <
class Matrix,
class Vector>
102 return (this->root_ && (this->matrixA_->getComm()->getSize() == 1) && is_contiguous_);
105 template<
class Matrix,
class Vector>
111 #ifdef HAVE_AMESOS2_TIMERS
112 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
119 template <
class Matrix,
class Vector>
127 ShyLUbasker->SetThreads(num_threads);
135 if ( single_proc_optimization() ) {
137 host_ordinal_type_array sp_rowptr;
138 host_ordinal_type_array sp_colind;
140 this->matrixA_->returnRowPtr_kokkos_view(sp_rowptr);
141 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr.data() ==
nullptr,
142 std::runtime_error,
"Amesos2 Runtime Error: sp_rowptr returned null ");
143 this->matrixA_->returnColInd_kokkos_view(sp_colind);
144 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind.data() ==
nullptr,
145 std::runtime_error,
"Amesos2 Runtime Error: sp_colind returned null ");
147 host_value_type_array hsp_values;
148 this->matrixA_->returnValues_kokkos_view(hsp_values);
149 shylubasker_dtype * sp_values = function_map::convert_scalar(hsp_values.data());
151 TEUCHOS_TEST_FOR_EXCEPTION(sp_values ==
nullptr,
152 std::runtime_error,
"Amesos2 Runtime Error: sp_values returned null ");
155 info = ShyLUbasker->Symbolic(this->globalNumRows_,
156 this->globalNumCols_,
157 this->globalNumNonZeros_,
163 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
164 std::runtime_error,
"Error in ShyLUBasker Symbolic");
169 shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
170 info = ShyLUbasker->Symbolic(this->globalNumRows_,
171 this->globalNumCols_,
172 this->globalNumNonZeros_,
177 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
178 std::runtime_error,
"Error in ShyLUBasker Symbolic");
184 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
189 template <
class Matrix,
class Vector>
198 #ifdef HAVE_AMESOS2_TIMERS
199 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
207 if ( single_proc_optimization() ) {
209 host_ordinal_type_array sp_rowptr;
210 host_ordinal_type_array sp_colind;
211 this->matrixA_->returnRowPtr_kokkos_view(sp_rowptr);
212 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr.data() ==
nullptr,
213 std::runtime_error,
"Amesos2 Runtime Error: sp_rowptr returned null ");
214 this->matrixA_->returnColInd_kokkos_view(sp_colind);
215 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind.data() ==
nullptr,
216 std::runtime_error,
"Amesos2 Runtime Error: sp_colind returned null ");
218 host_value_type_array hsp_values;
219 this->matrixA_->returnValues_kokkos_view(hsp_values);
220 shylubasker_dtype * sp_values = function_map::convert_scalar(hsp_values.data());
223 TEUCHOS_TEST_FOR_EXCEPTION(sp_values ==
nullptr,
224 std::runtime_error,
"Amesos2 Runtime Error: sp_values returned null ");
227 info = ShyLUbasker->Factor( this->globalNumRows_,
228 this->globalNumCols_,
229 this->globalNumNonZeros_,
237 shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
238 info = ShyLUbasker->Factor(this->globalNumRows_,
239 this->globalNumCols_,
240 this->globalNumNonZeros_,
248 local_ordinal_type blnnz = local_ordinal_type(0);
249 local_ordinal_type bunnz = local_ordinal_type(0);
250 ShyLUbasker->GetLnnz(blnnz);
251 ShyLUbasker->GetUnnz(bunnz);
255 this->setNnzLU( as<size_t>( blnnz + bunnz ) );
261 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
264 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
265 std::runtime_error,
" ShyLUBasker::numericFactorization failed.");
271 template <
class Matrix,
class Vector>
281 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
282 const size_t nrhs = X->getGlobalNumVectors();
284 const bool ShyluBaskerTransposeRequest = this->control_.useTranspose_;
285 const bool initialize_data =
true;
286 const bool do_not_initialize_data =
false;
287 bool use_gather = use_gather_;
288 use_gather = (use_gather && this->matrixA_->getComm()->getSize() > 1);
289 use_gather = (use_gather && (std::is_same<scalar_type, float>::value || std::is_same<scalar_type, double>::value));
291 #ifdef HAVE_AMESOS2_TIMERS
292 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
294 if ( single_proc_optimization() && nrhs == 1 ) {
297 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
298 host_solve_array_t>::do_get(initialize_data, B, bValues_, as<size_t>(ld_rhs));
300 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
301 host_solve_array_t>::do_get(do_not_initialize_data, X, xValues_, as<size_t>(ld_rhs));
306 int rval = B->gather(bValues_, this->perm_g2l, this->recvCountRows, this->recvDisplRows,
307 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
309 X->gather(xValues_, this->perm_g2l, this->recvCountRows, this->recvDisplRows,
310 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
316 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
317 host_solve_array_t>::do_get(initialize_data, B, bValues_,
319 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
320 this->rowIndexBase_);
323 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
324 host_solve_array_t>::do_get(do_not_initialize_data, X, xValues_,
326 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
327 this->rowIndexBase_);
333 #ifdef HAVE_AMESOS2_TIMERS
334 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
337 shylubasker_dtype * pxValues = function_map::convert_scalar(xValues_.data());
338 shylubasker_dtype * pbValues = function_map::convert_scalar(bValues_.data());
339 if (!ShyluBaskerTransposeRequest)
340 ierr = ShyLUbasker->Solve(nrhs, pbValues, pxValues);
342 ierr = ShyLUbasker->Solve(nrhs, pbValues, pxValues,
true);
345 Teuchos::broadcast(*(this->getComm()), 0, &ierr);
347 TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
349 "Encountered zero diag element at: " << ierr);
350 TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
352 "Could not alloc needed working memory for solve" );
354 #ifdef HAVE_AMESOS2_TIMERS
355 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
358 int rval = X->scatter(xValues_, this->perm_g2l, this->recvCountRows, this->recvDisplRows,
359 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
360 if (rval != 0) use_gather =
false;
363 Util::put_1d_data_helper_kokkos_view<
366 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
373 template <
class Matrix,
class Vector>
378 return( this->globalNumRows_ == this->globalNumCols_ );
382 template <
class Matrix,
class Vector>
387 using Teuchos::getIntegralValue;
388 using Teuchos::ParameterEntryValidator;
390 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
392 if(parameterList->isParameter(
"IsContiguous"))
394 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
397 if(parameterList->isParameter(
"UseCustomGather"))
399 use_gather_ = parameterList->get<
bool>(
"UseCustomGather");
402 if(parameterList->isParameter(
"num_threads"))
404 num_threads = parameterList->get<
int>(
"num_threads");
406 if(parameterList->isParameter(
"pivot"))
408 ShyLUbasker->Options.no_pivot = (!parameterList->get<
bool>(
"pivot"));
410 if(parameterList->isParameter(
"delayed pivot"))
412 ShyLUbasker->Options.static_delayed_pivot = (parameterList->get<
int>(
"delayed pivot"));
414 if(parameterList->isParameter(
"pivot_tol"))
416 ShyLUbasker->Options.pivot_tol = parameterList->get<
double>(
"pivot_tol");
418 if(parameterList->isParameter(
"symmetric"))
420 ShyLUbasker->Options.symmetric = parameterList->get<
bool>(
"symmetric");
422 if(parameterList->isParameter(
"realloc"))
424 ShyLUbasker->Options.realloc = parameterList->get<
bool>(
"realloc");
426 if(parameterList->isParameter(
"verbose"))
428 ShyLUbasker->Options.verbose = parameterList->get<
bool>(
"verbose");
430 if(parameterList->isParameter(
"verbose_matrix"))
432 ShyLUbasker->Options.verbose_matrix_out = parameterList->get<
bool>(
"verbose_matrix");
434 if(parameterList->isParameter(
"btf"))
436 ShyLUbasker->Options.btf = parameterList->get<
bool>(
"btf");
438 if(parameterList->isParameter(
"use_metis"))
440 ShyLUbasker->Options.use_metis = parameterList->get<
bool>(
"use_metis");
442 if(parameterList->isParameter(
"use_nodeNDP"))
444 ShyLUbasker->Options.use_nodeNDP = parameterList->get<
bool>(
"use_nodeNDP");
446 if(parameterList->isParameter(
"run_nd_on_leaves"))
448 ShyLUbasker->Options.run_nd_on_leaves = parameterList->get<
bool>(
"run_nd_on_leaves");
450 if(parameterList->isParameter(
"run_amd_on_leaves"))
452 ShyLUbasker->Options.run_amd_on_leaves = parameterList->get<
bool>(
"run_amd_on_leaves");
454 if(parameterList->isParameter(
"amd_on_blocks"))
456 ShyLUbasker->Options.amd_dom = parameterList->get<
bool>(
"amd_on_blocks");
458 if(parameterList->isParameter(
"transpose"))
461 const auto transpose = parameterList->get<
bool>(
"transpose");
463 this->control_.useTranspose_ =
true;
465 if(parameterList->isParameter(
"use_sequential_diag_facto"))
467 ShyLUbasker->Options.use_sequential_diag_facto = parameterList->get<
bool>(
"use_sequential_diag_facto");
469 if(parameterList->isParameter(
"user_fill"))
471 ShyLUbasker->Options.user_fill = parameterList->get<
double>(
"user_fill");
473 if(parameterList->isParameter(
"prune"))
475 ShyLUbasker->Options.prune = parameterList->get<
bool>(
"prune");
477 if(parameterList->isParameter(
"replace_zero_pivot"))
479 ShyLUbasker->Options.replace_zero_pivot = parameterList->get<
bool>(
"replace_zero_pivot");
481 if(parameterList->isParameter(
"replace_tiny_pivot"))
483 ShyLUbasker->Options.replace_tiny_pivot = parameterList->get<
bool>(
"replace_tiny_pivot");
485 if(parameterList->isParameter(
"btf_matching"))
487 ShyLUbasker->Options.btf_matching = parameterList->get<
int>(
"btf_matching");
488 if (ShyLUbasker->Options.btf_matching == 1 || ShyLUbasker->Options.btf_matching == 2) {
489 ShyLUbasker->Options.matching =
true;
491 ShyLUbasker->Options.matching =
false;
494 if(parameterList->isParameter(
"blk_matching"))
496 ShyLUbasker->Options.blk_matching = parameterList->get<
int>(
"blk_matching");
498 if(parameterList->isParameter(
"matrix_scaling"))
500 ShyLUbasker->Options.matrix_scaling = parameterList->get<
int>(
"matrix_scaling");
502 if(parameterList->isParameter(
"min_block_size"))
504 ShyLUbasker->Options.min_block_size = parameterList->get<
int>(
"min_block_size");
508 template <
class Matrix,
class Vector>
509 Teuchos::RCP<const Teuchos::ParameterList>
512 using Teuchos::ParameterList;
514 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
516 if( is_null(valid_params) )
518 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
519 pl->set(
"IsContiguous",
true,
520 "Are GIDs contiguous");
521 pl->set(
"UseCustomGather",
true,
522 "Use Matrix-gather routine");
523 pl->set(
"num_threads", 1,
524 "Number of threads");
525 pl->set(
"pivot",
false,
527 pl->set(
"delayed pivot", 0,
528 "Apply static delayed pivot on a big block");
529 pl->set(
"pivot_tol", .0001,
530 "Tolerance before pivot, currently not used");
531 pl->set(
"symmetric",
false,
532 "Should Symbolic assume symmetric nonzero pattern");
533 pl->set(
"realloc" ,
false,
534 "Should realloc space if not enough");
535 pl->set(
"verbose",
false,
536 "Information about factoring");
537 pl->set(
"verbose_matrix",
false,
538 "Give Permuted Matrices");
541 pl->set(
"prune",
false,
542 "Use prune on BTF blocks (Not Supported)");
543 pl->set(
"btf_matching", 2,
544 "Matching option for BTF: 0 = none, 1 = Basker, 2 = Trilinos (default), (3 = MC64 if enabled)");
545 pl->set(
"blk_matching", 1,
546 "Matching optioon for block: 0 = none, 1 or anything else = Basker (default), (2 = MC64 if enabled)");
547 pl->set(
"matrix_scaling", 0,
548 "Use matrix scaling to biig A BTF block: 0 = no-scaling, 1 = symmetric diagonal scaling, 2 = row-max, and then col-max scaling");
549 pl->set(
"min_block_size", 0,
550 "Size of the minimum diagonal blocks");
551 pl->set(
"replace_zero_pivot",
true,
552 "Replace zero pivots during the numerical factorization");
553 pl->set(
"replace_tiny_pivot",
false,
554 "Replace tiny pivots during the numerical factorization");
555 pl->set(
"use_metis",
true,
557 pl->set(
"use_nodeNDP",
true,
558 "Use nodeND to compute ND partition");
559 pl->set(
"run_nd_on_leaves",
false,
560 "Run ND on the final leaf-nodes for ND factorization");
561 pl->set(
"run_amd_on_leaves",
false,
562 "Run AMD on the final leaf-nodes for ND factorization");
563 pl->set(
"amd_on_blocks",
true,
564 "Run AMD on each diagonal blocks");
565 pl->set(
"transpose",
false,
566 "Solve the transpose A");
567 pl->set(
"use_sequential_diag_facto",
false,
568 "Use sequential algorithm to factor each diagonal block");
569 pl->set(
"user_fill", (
double)BASKER_FILL_USER,
570 "User-provided padding for the fill ratio");
577 template <
class Matrix,
class Vector>
582 if(current_phase == SOLVE || current_phase == PREORDERING )
return(
false );
584 #ifdef HAVE_AMESOS2_TIMERS
585 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
590 if ( single_proc_optimization() ) {
598 if( this->root_ && current_phase == SYMBFACT )
600 Kokkos::resize(nzvals_view_, this->globalNumNonZeros_);
601 Kokkos::resize(rowind_view_, this->globalNumNonZeros_);
602 Kokkos::resize(colptr_view_, this->globalNumCols_ + 1);
605 local_ordinal_type nnz_ret = -1;
606 bool use_gather = use_gather_;
607 use_gather = (use_gather && this->matrixA_->getComm()->getSize() > 1);
608 use_gather = (use_gather && (std::is_same<scalar_type, float>::value || std::is_same<scalar_type, double>::value));
610 #ifdef HAVE_AMESOS2_TIMERS
611 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
614 bool column_major =
true;
615 if (!is_contiguous_) {
616 auto contig_mat = this->matrixA_->reindex(this->contig_rowmap_, this->contig_colmap_, current_phase);
617 nnz_ret = contig_mat->gather(nzvals_view_, rowind_view_, colptr_view_, this->perm_g2l, this->recvCountRows, this->recvDisplRows, this->recvCounts, this->recvDispls,
618 this->transpose_map, this->nzvals_t, column_major, current_phase);
620 nnz_ret = this->matrixA_->gather(nzvals_view_, rowind_view_, colptr_view_, this->perm_g2l, this->recvCountRows, this->recvDisplRows, this->recvCounts, this->recvDispls,
621 this->transpose_map, this->nzvals_t, column_major, current_phase);
625 if (nnz_ret < 0) use_gather =
false;
630 ::do_get(this->matrixA_.ptr(), nzvals_view_, rowind_view_, colptr_view_, nnz_ret,
631 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
633 this->rowIndexBase_);
638 if (use_gather || this->root_) {
639 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
641 "Amesos2_ShyLUBasker loadA_impl: Did not get the expected number of non-zero vals("
642 +std::to_string(nnz_ret)+
" vs "+std::to_string(this->globalNumNonZeros_)+
")");
649 template<
class Matrix,
class Vector>
655 #endif // AMESOS2_SHYLUBASKER_DEF_HPP
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:614
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:273
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:510
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_ShyLUBasker_def.hpp:107
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:375
Amesos2 ShyLUBasker declarations.
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_ShyLUBasker_def.hpp:579
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:42
int numericFactorization_impl()
ShyLUBasker specific numeric factorization.
Definition: Amesos2_ShyLUBasker_def.hpp:191
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:101