Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_ShyLUBasker_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Amesos2: Templated Direct Sparse Solver Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//
42// @HEADER
43
54#ifndef AMESOS2_SHYLUBASKER_DEF_HPP
55#define AMESOS2_SHYLUBASKER_DEF_HPP
56
57#include <Teuchos_Tuple.hpp>
58#include <Teuchos_ParameterList.hpp>
59#include <Teuchos_StandardParameterEntryValidators.hpp>
60
63
64namespace Amesos2 {
65
66
67template <class Matrix, class Vector>
68ShyLUBasker<Matrix,Vector>::ShyLUBasker(
69 Teuchos::RCP<const Matrix> A,
70 Teuchos::RCP<Vector> X,
71 Teuchos::RCP<const Vector> B )
72 : SolverCore<Amesos2::ShyLUBasker,Matrix,Vector>(A, X, B)
73 , is_contiguous_(true)
74{
75
76 //Nothing
77
78 // Override some default options
79 // TODO: use data_ here to init
80#if defined(HAVE_AMESOS2_KOKKOS) && defined(KOKKOS_ENABLE_OPENMP)
81 /*
82 static_assert(std::is_same<kokkos_exe,Kokkos::OpenMP>::value,
83 "Kokkos node type not supported by experimental ShyLUBasker Amesos2");
84 */
85 typedef Kokkos::OpenMP Exe_Space;
86
87 ShyLUbasker = new ::BaskerNS::BaskerTrilinosInterface<local_ordinal_type, shylubasker_dtype, Exe_Space>();
88 ShyLUbasker->Options.no_pivot = BASKER_FALSE;
89 ShyLUbasker->Options.static_delayed_pivot = 0;
90 ShyLUbasker->Options.symmetric = BASKER_FALSE;
91 ShyLUbasker->Options.realloc = BASKER_TRUE;
92 ShyLUbasker->Options.verbose = BASKER_FALSE;
93 ShyLUbasker->Options.prune = BASKER_TRUE;
94 ShyLUbasker->Options.btf_matching = 2; // use cardinary matching from Trilinos, globally
95 ShyLUbasker->Options.blk_matching = 1; // use max-weight matching from Basker on each diagonal block
96 ShyLUbasker->Options.matrix_scaling = 0; // use matrix scaling on a big A block
97 ShyLUbasker->Options.min_block_size = 0; // no merging small blocks
98 ShyLUbasker->Options.amd_dom = BASKER_TRUE; // use block-wise AMD
99 ShyLUbasker->Options.use_metis = BASKER_TRUE; // use scotch/metis for ND (TODO: should METIS optional?)
100 ShyLUbasker->Options.use_nodeNDP = BASKER_TRUE; // use nodeNDP to compute ND partition
101 ShyLUbasker->Options.run_nd_on_leaves = BASKER_TRUE; // run ND on the final leaf-nodes
102 ShyLUbasker->Options.run_amd_on_leaves = BASKER_FALSE; // run AMD on the final leaf-nodes
103 ShyLUbasker->Options.transpose = BASKER_FALSE;
104 ShyLUbasker->Options.replace_tiny_pivot = BASKER_FALSE;
105 ShyLUbasker->Options.verbose_matrix_out = BASKER_FALSE;
106
107 ShyLUbasker->Options.user_fill = (double)BASKER_FILL_USER;
108 ShyLUbasker->Options.use_sequential_diag_facto = BASKER_FALSE;
109#ifdef KOKKOS_ENABLE_DEPRECATED_CODE
110 num_threads = Kokkos::OpenMP::max_hardware_threads();
111#else
112 num_threads = Kokkos::OpenMP::impl_max_hardware_threads();
113#endif
114
115 ShyLUbaskerTr = new ::BaskerNS::BaskerTrilinosInterface<local_ordinal_type, shylubasker_dtype, Exe_Space>();
116 ShyLUbaskerTr->Options.no_pivot = BASKER_FALSE;
117 ShyLUbaskerTr->Options.static_delayed_pivot = 0;
118 ShyLUbaskerTr->Options.symmetric = BASKER_FALSE;
119 ShyLUbaskerTr->Options.realloc = BASKER_TRUE;
120 ShyLUbaskerTr->Options.verbose = BASKER_FALSE;
121 ShyLUbaskerTr->Options.prune = BASKER_TRUE;
122 ShyLUbaskerTr->Options.btf_matching = 2; // use cardinary matching from Trilinos, globally
123 ShyLUbaskerTr->Options.blk_matching = 1; // use max-weight matching from Basker on each diagonal block
124 ShyLUbaskerTr->Options.matrix_scaling = 0; // use matrix scaling on a big A block
125 ShyLUbaskerTr->Options.min_block_size = 0; // no merging small blocks
126 ShyLUbaskerTr->Options.amd_dom = BASKER_TRUE; // use block-wise AMD
127 ShyLUbaskerTr->Options.use_metis = BASKER_TRUE; // use scotch/metis for ND (TODO: should METIS optional?)
128 ShyLUbaskerTr->Options.use_nodeNDP = BASKER_TRUE; // use nodeNDP to compute ND partition
129 ShyLUbaskerTr->Options.run_nd_on_leaves = BASKER_TRUE; // run ND on the final leaf-nodes
130 ShyLUbaskerTr->Options.run_amd_on_leaves = BASKER_FALSE; // run ND on the final leaf-nodes
131 ShyLUbaskerTr->Options.transpose = BASKER_TRUE;
132 ShyLUbaskerTr->Options.replace_tiny_pivot = BASKER_FALSE;
133 ShyLUbaskerTr->Options.verbose_matrix_out = BASKER_FALSE;
134
135 ShyLUbaskerTr->Options.user_fill = (double)BASKER_FILL_USER;
136 ShyLUbaskerTr->Options.use_sequential_diag_facto = BASKER_FALSE;
137#else
138 TEUCHOS_TEST_FOR_EXCEPTION(1 != 0,
139 std::runtime_error,
140 "Amesos2_ShyLUBasker Exception: Do not have supported Kokkos node type (OpenMP) enabled for ShyLUBasker");
141#endif
142}
143
144
145template <class Matrix, class Vector>
146ShyLUBasker<Matrix,Vector>::~ShyLUBasker( )
147{
148 /* ShyLUBasker will cleanup its own internal memory*/
149#if defined(HAVE_AMESOS2_KOKKOS) && defined(KOKKOS_ENABLE_OPENMP)
150 ShyLUbasker->Finalize();
151 ShyLUbaskerTr->Finalize();
152 delete ShyLUbasker;
153 delete ShyLUbaskerTr;
154#endif
155}
156
157template <class Matrix, class Vector>
158bool
160 return (this->root_ && (this->matrixA_->getComm()->getSize() == 1) && is_contiguous_);
161}
162
163template<class Matrix, class Vector>
164int
166{
167 /* TODO: Define what it means for ShyLUBasker
168 */
169#ifdef HAVE_AMESOS2_TIMERS
170 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
171#endif
172
173 return(0);
174}
175
176
177template <class Matrix, class Vector>
178int
180{
181
182 int info = 0;
183 if(this->root_)
184 {
185 ShyLUbasker->SetThreads(num_threads);
186 ShyLUbaskerTr->SetThreads(num_threads);
187
188
189 // NDE: Special case
190 // Rather than going through the Amesos2 machinery to convert the matrixA_ CRS pointer data to CCS and store in Teuchos::Arrays,
191 // in this special case we pass the CRS raw pointers directly to ShyLUBasker which copies+transposes+sorts the data for CCS format
192 // loadA_impl is essentially an empty function in this case, as the raw pointers are handled here and similarly in Symbolic
193
194 if ( single_proc_optimization() ) {
195
196 host_ordinal_type_array sp_rowptr;
197 host_ordinal_type_array sp_colind;
198 // this needs to be checked during loadA_impl...
199 this->matrixA_->returnRowPtr_kokkos_view(sp_rowptr);
200 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr.data() == nullptr,
201 std::runtime_error, "Amesos2 Runtime Error: sp_rowptr returned null ");
202 this->matrixA_->returnColInd_kokkos_view(sp_colind);
203 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind.data() == nullptr,
204 std::runtime_error, "Amesos2 Runtime Error: sp_colind returned null ");
205
206 host_value_type_array hsp_values;
207 this->matrixA_->returnValues_kokkos_view(hsp_values);
208 shylubasker_dtype * sp_values = function_map::convert_scalar(hsp_values.data());
209 //shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
210 TEUCHOS_TEST_FOR_EXCEPTION(sp_values == nullptr,
211 std::runtime_error, "Amesos2 Runtime Error: sp_values returned null ");
212
213 // In this case, colptr_, rowind_, nzvals_ are invalid
214 info = ShyLUbasker->Symbolic(this->globalNumRows_,
215 this->globalNumCols_,
216 this->globalNumNonZeros_,
217 sp_rowptr.data(),
218 sp_colind.data(),
219 sp_values,
220 true);
221
222 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
223 std::runtime_error, "Error in ShyLUBasker Symbolic");
224
225 if (info == BASKER_SUCCESS) {
226 info = ShyLUbaskerTr->Symbolic(this->globalNumRows_,
227 this->globalNumCols_,
228 this->globalNumNonZeros_,
229 sp_rowptr.data(),
230 sp_colind.data(),
231 sp_values,
232 true);
233
234 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
235 std::runtime_error, "Error in ShyLUBaskerTr Symbolic");
236 }
237 }
238 else
239 { //follow original code path if conditions not met
240 // In this case, loadA_impl updates colptr_, rowind_, nzvals_
241 shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
242 info = ShyLUbasker->Symbolic(this->globalNumRows_,
243 this->globalNumCols_,
244 this->globalNumNonZeros_,
245 colptr_view_.data(),
246 rowind_view_.data(),
247 sp_values);
248
249 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
250 std::runtime_error, "Error in ShyLUBasker Symbolic");
251
252 if (info == BASKER_SUCCESS) {
253 info = ShyLUbaskerTr->Symbolic(this->globalNumRows_,
254 this->globalNumCols_,
255 this->globalNumNonZeros_,
256 colptr_view_.data(),
257 rowind_view_.data(),
258 sp_values,
259 false);
260
261 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
262 std::runtime_error, "Error in ShyLUBaskerTr Symbolic");
263 }
264 }
265 } // end if (this->root_)
266 /*No symbolic factoriztion*/
267
268 /* All processes should have the same error code */
269 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
270 return(info);
271}
272
273
274template <class Matrix, class Vector>
275int
277{
278 using Teuchos::as;
279
280 int info = 0;
281 if ( this->root_ ){
282 { // Do factorization
283#ifdef HAVE_AMESOS2_TIMERS
284 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
285#endif
286
287
288 // NDE: Special case
289 // Rather than going through the Amesos2 machinery to convert the matrixA_ CRS pointer data to CCS and store in Teuchos::Arrays,
290 // in this special case we pass the CRS raw pointers directly to ShyLUBasker which copies+transposes+sorts the data for CCS format
291 // loadA_impl is essentially an empty function in this case, as the raw pointers are handled here and similarly in Symbolic
292
293 if ( single_proc_optimization() ) {
294
295 host_ordinal_type_array sp_rowptr;
296 host_ordinal_type_array sp_colind;
297 this->matrixA_->returnRowPtr_kokkos_view(sp_rowptr);
298 TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr.data() == nullptr,
299 std::runtime_error, "Amesos2 Runtime Error: sp_rowptr returned null ");
300 this->matrixA_->returnColInd_kokkos_view(sp_colind);
301 TEUCHOS_TEST_FOR_EXCEPTION(sp_colind.data() == nullptr,
302 std::runtime_error, "Amesos2 Runtime Error: sp_colind returned null ");
303
304 host_value_type_array hsp_values;
305 this->matrixA_->returnValues_kokkos_view(hsp_values);
306 shylubasker_dtype * sp_values = function_map::convert_scalar(hsp_values.data());
307 //shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
308
309 TEUCHOS_TEST_FOR_EXCEPTION(sp_values == nullptr,
310 std::runtime_error, "Amesos2 Runtime Error: sp_values returned null ");
311
312 // In this case, colptr_, rowind_, nzvals_ are invalid
313 info = ShyLUbasker->Factor( this->globalNumRows_,
314 this->globalNumCols_,
315 this->globalNumNonZeros_,
316 sp_rowptr.data(),
317 sp_colind.data(),
318 sp_values);
319
320 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
321 std::runtime_error, "Error ShyLUBasker Factor");
322
323 if (info == 0) {
324 info = ShyLUbaskerTr->Factor( this->globalNumRows_,
325 this->globalNumCols_,
326 this->globalNumNonZeros_,
327 sp_rowptr.data(),
328 sp_colind.data(),
329 sp_values);
330
331 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
332 std::runtime_error, "Error ShyLUBaskerTr Factor");
333 }
334 }
335 else
336 {
337 // In this case, loadA_impl updates colptr_, rowind_, nzvals_
338 shylubasker_dtype * sp_values = function_map::convert_scalar(nzvals_view_.data());
339 info = ShyLUbasker->Factor(this->globalNumRows_,
340 this->globalNumCols_,
341 this->globalNumNonZeros_,
342 colptr_view_.data(),
343 rowind_view_.data(),
344 sp_values);
345
346 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
347 std::runtime_error, "Error ShyLUBasker Factor");
348
349 if (info == 0) {
350 info = ShyLUbaskerTr->Factor(this->globalNumRows_,
351 this->globalNumCols_,
352 this->globalNumNonZeros_,
353 colptr_view_.data(),
354 rowind_view_.data(),
355 sp_values);
356
357 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
358 std::runtime_error, "Error ShyLUBaskerTr Factor");
359 }
360 //We need to handle the realloc options
361 }
362
363 //ShyLUbasker->DEBUG_PRINT();
364
365 local_ordinal_type blnnz = local_ordinal_type(0);
366 local_ordinal_type bunnz = local_ordinal_type(0);
367 ShyLUbasker->GetLnnz(blnnz); // Add exception handling?
368 ShyLUbasker->GetUnnz(bunnz);
369
370 local_ordinal_type Trblnnz = local_ordinal_type(0);
371 local_ordinal_type Trbunnz = local_ordinal_type(0);
372 ShyLUbaskerTr->GetLnnz(Trblnnz); // Add exception handling?
373 ShyLUbaskerTr->GetUnnz(Trbunnz);
374
375 // This is set after numeric factorization complete as pivoting can be used;
376 // In this case, a discrepancy between symbolic and numeric nnz total can occur.
377 this->setNnzLU( as<size_t>( blnnz + bunnz ) );
378
379 } // end scope for timer
380 } // end if (this->root_)
381
382 /* All processes should have the same error code */
383 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
384
385 //global_size_type info_st = as<global_size_type>(info);
386 /* TODO : Proper error messages*/
387 TEUCHOS_TEST_FOR_EXCEPTION(info == -1,
388 std::runtime_error,
389 "ShyLUBasker: Could not alloc space for L and U");
390 TEUCHOS_TEST_FOR_EXCEPTION(info == -2,
391 std::runtime_error,
392 "ShyLUBasker: Could not alloc needed work space");
393 TEUCHOS_TEST_FOR_EXCEPTION(info == -3,
394 std::runtime_error,
395 "ShyLUBasker: Could not alloc additional memory needed for L and U");
396 TEUCHOS_TEST_FOR_EXCEPTION(info > 0,
397 std::runtime_error,
398 "ShyLUBasker: Zero pivot found at: " << info );
399
400 return(info);
401}
402
403
404template <class Matrix, class Vector>
405int
407 const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
408 const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
409{
410#ifdef HAVE_AMESOS2_TIMERS
411 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
412#endif
413
414 int ierr = 0; // returned error code
415
416 using Teuchos::as;
417
418 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
419 const size_t nrhs = X->getGlobalNumVectors();
420
421 bool ShyluBaskerTransposeRequest = this->control_.useTranspose_;
422 const bool initialize_data = true;
423 const bool do_not_initialize_data = false;
424
425 if ( single_proc_optimization() && nrhs == 1 ) {
426
427 // no msp creation
428 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
429 host_solve_array_t>::do_get(initialize_data, B, bValues_, as<size_t>(ld_rhs));
430
431 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
432 host_solve_array_t>::do_get(do_not_initialize_data, X, xValues_, as<size_t>(ld_rhs));
433
434 } // end if ( single_proc_optimization() && nrhs == 1 )
435 else {
436
437 if ( is_contiguous_ == true ) {
438 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
439 host_solve_array_t>::do_get(initialize_data, B, bValues_, as<size_t>(ld_rhs), ROOTED, this->rowIndexBase_);
440 }
441 else {
442 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
443 host_solve_array_t>::do_get(initialize_data, B, bValues_, as<size_t>(ld_rhs), CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
444 }
445 // See Amesos2_Tacho_def.hpp for notes on why we 'get' x here.
446 if ( is_contiguous_ == true ) {
447 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
448 host_solve_array_t>::do_get(do_not_initialize_data, X, xValues_, as<size_t>(ld_rhs), ROOTED, this->rowIndexBase_);
449 }
450 else {
451 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
452 host_solve_array_t>::do_get(do_not_initialize_data, X, xValues_, as<size_t>(ld_rhs), CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
453 }
454 }
455
456 if ( this->root_ ) { // do solve
457 shylubasker_dtype * pxValues = function_map::convert_scalar(xValues_.data());
458 shylubasker_dtype * pbValues = function_map::convert_scalar(bValues_.data());
459 if (!ShyluBaskerTransposeRequest)
460 ierr = ShyLUbasker->Solve(nrhs, pbValues, pxValues);
461 else
462 ierr = ShyLUbaskerTr->Solve(nrhs, pbValues, pxValues);
463 }
464
465 /* All processes should have the same error code */
466 Teuchos::broadcast(*(this->getComm()), 0, &ierr);
467
468 TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
469 std::runtime_error,
470 "Encountered zero diag element at: " << ierr);
471 TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
472 std::runtime_error,
473 "Could not alloc needed working memory for solve" );
474
475 if ( is_contiguous_ == true ) {
476 Util::put_1d_data_helper_kokkos_view<
477 MultiVecAdapter<Vector>,host_solve_array_t>::do_put(X, xValues_,
478 as<size_t>(ld_rhs),
479 ROOTED);
480 }
481 else {
482 Util::put_1d_data_helper_kokkos_view<
483 MultiVecAdapter<Vector>,host_solve_array_t>::do_put(X, xValues_,
484 as<size_t>(ld_rhs),
486 }
487
488 return(ierr);
489}
490
491
492template <class Matrix, class Vector>
493bool
495{
496 // The ShyLUBasker can only handle square for right now
497 return( this->globalNumRows_ == this->globalNumCols_ );
498}
499
500
501template <class Matrix, class Vector>
502void
503ShyLUBasker<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
504{
505 using Teuchos::RCP;
506 using Teuchos::getIntegralValue;
507 using Teuchos::ParameterEntryValidator;
508
509 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
510
511 if(parameterList->isParameter("IsContiguous"))
512 {
513 is_contiguous_ = parameterList->get<bool>("IsContiguous");
514 }
515
516 if(parameterList->isParameter("num_threads"))
517 {
518 num_threads = parameterList->get<int>("num_threads");
519 }
520 if(parameterList->isParameter("pivot"))
521 {
522 ShyLUbasker->Options.no_pivot = (!parameterList->get<bool>("pivot"));
523 ShyLUbaskerTr->Options.no_pivot = (!parameterList->get<bool>("pivot"));
524 }
525 if(parameterList->isParameter("delayed pivot"))
526 {
527 ShyLUbasker->Options.static_delayed_pivot = (parameterList->get<int>("delayed pivot"));
528 ShyLUbaskerTr->Options.static_delayed_pivot = (parameterList->get<int>("delayed pivot"));
529 }
530 if(parameterList->isParameter("pivot_tol"))
531 {
532 ShyLUbasker->Options.pivot_tol = parameterList->get<double>("pivot_tol");
533 ShyLUbaskerTr->Options.pivot_tol = parameterList->get<double>("pivot_tol");
534 }
535 if(parameterList->isParameter("symmetric"))
536 {
537 ShyLUbasker->Options.symmetric = parameterList->get<bool>("symmetric");
538 ShyLUbaskerTr->Options.symmetric = parameterList->get<bool>("symmetric");
539 }
540 if(parameterList->isParameter("realloc"))
541 {
542 ShyLUbasker->Options.realloc = parameterList->get<bool>("realloc");
543 ShyLUbaskerTr->Options.realloc = parameterList->get<bool>("realloc");
544 }
545 if(parameterList->isParameter("verbose"))
546 {
547 ShyLUbasker->Options.verbose = parameterList->get<bool>("verbose");
548 ShyLUbaskerTr->Options.verbose = parameterList->get<bool>("verbose");
549 }
550 if(parameterList->isParameter("verbose_matrix"))
551 {
552 ShyLUbasker->Options.verbose_matrix_out = parameterList->get<bool>("verbose_matrix");
553 ShyLUbaskerTr->Options.verbose_matrix_out = parameterList->get<bool>("verbose_matrix");
554 }
555 if(parameterList->isParameter("btf"))
556 {
557 ShyLUbasker->Options.btf = parameterList->get<bool>("btf");
558 ShyLUbaskerTr->Options.btf = parameterList->get<bool>("btf");
559 }
560 if(parameterList->isParameter("use_metis"))
561 {
562 ShyLUbasker->Options.use_metis = parameterList->get<bool>("use_metis");
563 ShyLUbaskerTr->Options.use_metis = parameterList->get<bool>("use_metis");
564 }
565 if(parameterList->isParameter("use_nodeNDP"))
566 {
567 ShyLUbasker->Options.use_nodeNDP = parameterList->get<bool>("use_nodeNDP");
568 ShyLUbaskerTr->Options.use_nodeNDP = parameterList->get<bool>("use_nodeNDP");
569 }
570 if(parameterList->isParameter("run_nd_on_leaves"))
571 {
572 ShyLUbasker->Options.run_nd_on_leaves = parameterList->get<bool>("run_nd_on_leaves");
573 ShyLUbaskerTr->Options.run_nd_on_leaves = parameterList->get<bool>("run_nd_on_leaves");
574 }
575 if(parameterList->isParameter("run_amd_on_leaves"))
576 {
577 ShyLUbasker->Options.run_amd_on_leaves = parameterList->get<bool>("run_amd_on_leaves");
578 ShyLUbaskerTr->Options.run_amd_on_leaves = parameterList->get<bool>("run_amd_on_leaves");
579 }
580 if(parameterList->isParameter("amd_on_blocks"))
581 {
582 ShyLUbasker->Options.amd_dom = parameterList->get<bool>("amd_on_blocks");
583 ShyLUbaskerTr->Options.amd_dom = parameterList->get<bool>("amd_on_blocks");
584 }
585 if(parameterList->isParameter("transpose"))
586 {
587 // NDE: set transpose vs non-transpose mode as bool; track separate shylubasker objects
588 const auto transpose = parameterList->get<bool>("transpose");
589 if (transpose == true)
590 this->control_.useTranspose_ = true;
591 //ShyLUbasker->Options.transpose = parameterList->get<bool>("transpose");
592 //ShyLUbaskerTr->Options.transpose = parameterList->get<bool>("transpose");
593 }
594 if(parameterList->isParameter("use_sequential_diag_facto"))
595 {
596 ShyLUbasker->Options.use_sequential_diag_facto = parameterList->get<bool>("use_sequential_diag_facto");
597 ShyLUbaskerTr->Options.use_sequential_diag_facto = parameterList->get<bool>("use_sequential_diag_facto");
598 }
599 if(parameterList->isParameter("user_fill"))
600 {
601 ShyLUbasker->Options.user_fill = parameterList->get<double>("user_fill");
602 ShyLUbaskerTr->Options.user_fill = parameterList->get<double>("user_fill");
603 }
604 if(parameterList->isParameter("prune"))
605 {
606 ShyLUbasker->Options.prune = parameterList->get<bool>("prune");
607 ShyLUbaskerTr->Options.prune = parameterList->get<bool>("prune");
608 }
609 if(parameterList->isParameter("replace_tiny_pivot"))
610 {
611 ShyLUbasker->Options.replace_tiny_pivot = parameterList->get<bool>("replace_tiny_pivot");
612 ShyLUbaskerTr->Options.replace_tiny_pivot = parameterList->get<bool>("replace_tiny_pivot");
613 }
614 if(parameterList->isParameter("btf_matching"))
615 {
616 ShyLUbasker->Options.btf_matching = parameterList->get<int>("btf_matching");
617 ShyLUbaskerTr->Options.btf_matching = parameterList->get<int>("btf_matching");
618 if (ShyLUbasker->Options.btf_matching == 1 || ShyLUbasker->Options.btf_matching == 2) {
619 ShyLUbasker->Options.matching = true;
620 ShyLUbaskerTr->Options.matching = true;
621 } else {
622 ShyLUbasker->Options.matching = false;
623 ShyLUbaskerTr->Options.matching = false;
624 }
625 }
626 if(parameterList->isParameter("blk_matching"))
627 {
628 ShyLUbasker->Options.blk_matching = parameterList->get<int>("blk_matching");
629 ShyLUbaskerTr->Options.blk_matching = parameterList->get<int>("blk_matching");
630 }
631 if(parameterList->isParameter("matrix_scaling"))
632 {
633 ShyLUbasker->Options.matrix_scaling = parameterList->get<int>("matrix_scaling");
634 ShyLUbaskerTr->Options.matrix_scaling = parameterList->get<int>("matrix_scaling");
635 }
636 if(parameterList->isParameter("min_block_size"))
637 {
638 ShyLUbasker->Options.min_block_size = parameterList->get<int>("min_block_size");
639 ShyLUbaskerTr->Options.min_block_size = parameterList->get<int>("min_block_size");
640 }
641}
642
643template <class Matrix, class Vector>
644Teuchos::RCP<const Teuchos::ParameterList>
646{
647 using Teuchos::ParameterList;
648
649 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
650
651 if( is_null(valid_params) )
652 {
653 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
654 pl->set("IsContiguous", true,
655 "Are GIDs contiguous");
656 pl->set("num_threads", 1,
657 "Number of threads");
658 pl->set("pivot", false,
659 "Should not pivot");
660 pl->set("delayed pivot", 0,
661 "Apply static delayed pivot on a big block");
662 pl->set("pivot_tol", .0001,
663 "Tolerance before pivot, currently not used");
664 pl->set("symmetric", false,
665 "Should Symbolic assume symmetric nonzero pattern");
666 pl->set("realloc" , false,
667 "Should realloc space if not enough");
668 pl->set("verbose", false,
669 "Information about factoring");
670 pl->set("verbose_matrix", false,
671 "Give Permuted Matrices");
672 pl->set("btf", true,
673 "Use BTF ordering");
674 pl->set("prune", false,
675 "Use prune on BTF blocks (Not Supported)");
676 pl->set("btf_matching", 2,
677 "Matching option for BTF: 0 = none, 1 = Basker, 2 = Trilinos (default), (3 = MC64 if enabled)");
678 pl->set("blk_matching", 1,
679 "Matching optioon for block: 0 = none, 1 or anything else = Basker (default), (2 = MC64 if enabled)");
680 pl->set("matrix_scaling", 0,
681 "Use matrix scaling to biig A BTF block: 0 = no-scaling, 1 = symmetric diagonal scaling, 2 = row-max, and then col-max scaling");
682 pl->set("min_block_size", 0,
683 "Size of the minimum diagonal blocks");
684 pl->set("replace_tiny_pivot", false,
685 "Replace tiny pivots during the numerical factorization");
686 pl->set("use_metis", true,
687 "Use METIS for ND");
688 pl->set("use_nodeNDP", true,
689 "Use nodeND to compute ND partition");
690 pl->set("run_nd_on_leaves", false,
691 "Run ND on the final leaf-nodes for ND factorization");
692 pl->set("run_amd_on_leaves", false,
693 "Run AMD on the final leaf-nodes for ND factorization");
694 pl->set("amd_on_blocks", true,
695 "Run AMD on each diagonal blocks");
696 pl->set("transpose", false,
697 "Solve the transpose A");
698 pl->set("use_sequential_diag_facto", false,
699 "Use sequential algorithm to factor each diagonal block");
700 pl->set("user_fill", (double)BASKER_FILL_USER,
701 "User-provided padding for the fill ratio");
702 valid_params = pl;
703 }
704 return valid_params;
705}
706
707
708template <class Matrix, class Vector>
709bool
711{
712 using Teuchos::as;
713 if(current_phase == SOLVE) return (false);
714
715 #ifdef HAVE_AMESOS2_TIMERS
716 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
717 #endif
718
719
720 // NDE: Can clean up duplicated code with the #ifdef guards
721 if ( single_proc_optimization() ) {
722 // NDE: Nothing is done in this special case - CRS raw pointers are passed to SHYLUBASKER and transpose of copies handled there
723 // In this case, colptr_, rowind_, nzvals_ are invalid
724 }
725 else
726 {
727
728 // Only the root image needs storage allocated
729 if( this->root_ ){
730 Kokkos::resize(nzvals_view_, this->globalNumNonZeros_);
731 Kokkos::resize(rowind_view_, this->globalNumNonZeros_);
732 Kokkos::resize(colptr_view_, this->globalNumCols_ + 1); //this will be wrong for case of gapped col ids, e.g. 0,2,4,9; num_cols = 10 ([0,10)) but num GIDs = 4...
733 }
734
735 local_ordinal_type nnz_ret = 0;
736 {
737 #ifdef HAVE_AMESOS2_TIMERS
738 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
739 #endif
740
741 if ( is_contiguous_ == true ) {
743 MatrixAdapter<Matrix>, host_value_type_array, host_ordinal_type_array, host_ordinal_type_array>
744 ::do_get(this->matrixA_.ptr(), nzvals_view_, rowind_view_, colptr_view_,
745 nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_); // copies from matrixA_ to ShyLUBasker ConcreteSolver cp, ri, nzval members
746 }
747 else {
749 MatrixAdapter<Matrix>, host_value_type_array, host_ordinal_type_array, host_ordinal_type_array>
750 ::do_get(this->matrixA_.ptr(), nzvals_view_, rowind_view_, colptr_view_,
751 nnz_ret, CONTIGUOUS_AND_ROOTED, ARBITRARY, this->rowIndexBase_); // copies from matrixA_ to ShyLUBasker ConcreteSolver cp, ri, nzval members
752 }
753 }
754
755 if( this->root_ ){
756 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
757 std::runtime_error,
758 "Amesos2_ShyLUBasker loadA_impl: Did not get the expected number of non-zero vals");
759 }
760
761 } //end alternative path
762 return true;
763}
764
765
766template<class Matrix, class Vector>
767const char* ShyLUBasker<Matrix,Vector>::name = "ShyLUBasker";
768
769
770} // end namespace Amesos2
771
772#endif // AMESOS2_SHYLUBASKER_DEF_HPP
Amesos2 ShyLUBasker declarations.
@ ROOTED
Definition Amesos2_TypeDecl.hpp:127
@ CONTIGUOUS_AND_ROOTED
Definition Amesos2_TypeDecl.hpp:128
@ ARBITRARY
Definition Amesos2_TypeDecl.hpp:143
Amesos2 interface to the Baker package.
Definition Amesos2_ShyLUBasker_decl.hpp:77
EPhase
Used to indicate a phase in the direct solution.
Definition Amesos2_TypeDecl.hpp:65
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
A generic helper class for getting a CCS representation of a Matrix.
Definition Amesos2_Util.hpp:652