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 
64 namespace Amesos2 {
65 
66 
67 template <class Matrix, class Vector>
68 ShyLUBasker<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  , nzvals_() // initialize to empty arrays
74  , rowind_()
75  , colptr_()
76  , is_contiguous_(true)
77 {
78 
79  //Nothing
80 
81  // Override some default options
82  // TODO: use data_ here to init
83 #if defined(HAVE_AMESOS2_KOKKOS) && defined(KOKKOS_ENABLE_OPENMP)
84  /*
85  static_assert(std::is_same<kokkos_exe,Kokkos::OpenMP>::value,
86  "Kokkos node type not supported by experimental ShyLUBasker Amesos2");
87  */
88  typedef Kokkos::OpenMP Exe_Space;
89 
90  ShyLUbasker = new ::BaskerNS::BaskerTrilinosInterface<local_ordinal_type, slu_type, Exe_Space>();
91  ShyLUbasker->Options.no_pivot = BASKER_FALSE;
92  ShyLUbasker->Options.static_delayed_pivot = 0;
93  ShyLUbasker->Options.symmetric = BASKER_FALSE;
94  ShyLUbasker->Options.realloc = BASKER_TRUE;
95  ShyLUbasker->Options.verbose = BASKER_FALSE;
96  ShyLUbasker->Options.prune = BASKER_TRUE;
97  ShyLUbasker->Options.btf_matching = 2; // use cardinary matching from Trilinos, globally
98  ShyLUbasker->Options.blk_matching = 1; // use max-weight matching from Basker on each diagonal block
99  ShyLUbasker->Options.min_block_size = 0; // no merging small blocks
100  ShyLUbasker->Options.amd_dom = BASKER_TRUE; // use block-wise AMD
101  ShyLUbasker->Options.use_metis = BASKER_TRUE; // use scotch/metis for ND (TODO: should METIS optional?)
102  ShyLUbasker->Options.run_nd_on_leaves = BASKER_FALSE; // run ND on the final leaf-nodes
103  ShyLUbasker->Options.transpose = BASKER_FALSE;
104  ShyLUbasker->Options.replace_tiny_pivot = BASKER_TRUE;
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, slu_type, 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.min_block_size = 0; // no merging small blocks
125  ShyLUbaskerTr->Options.amd_dom = BASKER_TRUE; // use block-wise AMD
126  ShyLUbaskerTr->Options.use_metis = BASKER_TRUE; // use scotch/metis for ND (TODO: should METIS optional?)
127  ShyLUbaskerTr->Options.run_nd_on_leaves = BASKER_FALSE; // run ND on the final leaf-nodes
128  ShyLUbaskerTr->Options.transpose = BASKER_TRUE;
129  ShyLUbaskerTr->Options.replace_tiny_pivot = BASKER_TRUE;
130  ShyLUbaskerTr->Options.verbose_matrix_out = BASKER_FALSE;
131 
132  ShyLUbaskerTr->Options.user_fill = (double)BASKER_FILL_USER;
133  ShyLUbaskerTr->Options.use_sequential_diag_facto = BASKER_FALSE;
134 #else
135  TEUCHOS_TEST_FOR_EXCEPTION(1 != 0,
136  std::runtime_error,
137  "Amesos2_ShyLUBasker Exception: Do not have supported Kokkos node type (OpenMP) enabled for ShyLUBasker");
138 #endif
139 }
140 
141 
142 template <class Matrix, class Vector>
143 ShyLUBasker<Matrix,Vector>::~ShyLUBasker( )
144 {
145  /* ShyLUBasker will cleanup its own internal memory*/
146 #if defined(HAVE_AMESOS2_KOKKOS) && defined(KOKKOS_ENABLE_OPENMP)
147  ShyLUbasker->Finalize();
148  ShyLUbaskerTr->Finalize();
149  delete ShyLUbasker;
150  delete ShyLUbaskerTr;
151 #endif
152 }
153 
154 template <class Matrix, class Vector>
155 bool
157  return (this->root_ && (this->matrixA_->getComm()->getSize() == 1) && is_contiguous_);
158 }
159 
160 template<class Matrix, class Vector>
161 int
163 {
164  /* TODO: Define what it means for ShyLUBasker
165  */
166 #ifdef HAVE_AMESOS2_TIMERS
167  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
168 #endif
169 
170  return(0);
171 }
172 
173 
174 template <class Matrix, class Vector>
175 int
177 {
178 
179  int info = 0;
180  if(this->root_)
181  {
182  ShyLUbasker->SetThreads(num_threads);
183  ShyLUbaskerTr->SetThreads(num_threads);
184 
185 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
186  std::cout << "ShyLUBasker:: Before symbolic factorization" << std::endl;
187  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
188  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
189  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
190 #endif
191 
192  // NDE: Special case
193  // Rather than going through the Amesos2 machinery to convert the matrixA_ CRS pointer data to CCS and store in Teuchos::Arrays,
194  // in this special case we pass the CRS raw pointers directly to ShyLUBasker which copies+transposes+sorts the data for CCS format
195  // loadA_impl is essentially an empty function in this case, as the raw pointers are handled here and similarly in Symbolic
196 
197  if ( single_proc_optimization() ) {
198 
199  // this needs to be checked during loadA_impl...
200  auto sp_rowptr = this->matrixA_->returnRowPtr();
201  TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr == nullptr,
202  std::runtime_error, "Amesos2 Runtime Error: sp_rowptr returned null ");
203  auto sp_colind = this->matrixA_->returnColInd();
204  TEUCHOS_TEST_FOR_EXCEPTION(sp_colind == nullptr,
205  std::runtime_error, "Amesos2 Runtime Error: sp_colind returned null ");
206 #ifndef HAVE_TEUCHOS_COMPLEX
207  auto sp_values = this->matrixA_->returnValues();
208 #else
209  // NDE: 09/25/2017
210  // Cannot convert Kokkos::complex<T>* to std::complex<T>*; in this case, use reinterpret_cast
211  using complex_type = typename Util::getStdCplxType< magnitude_type, typename matrix_adapter_type::spmtx_vals_t >::type;
212  complex_type * sp_values = nullptr;
213  sp_values = reinterpret_cast< complex_type * > ( this->matrixA_->returnValues() );
214 #endif
215  TEUCHOS_TEST_FOR_EXCEPTION(sp_values == nullptr,
216  std::runtime_error, "Amesos2 Runtime Error: sp_values returned null ");
217 
218  // In this case, colptr_, rowind_, nzvals_ are invalid
219  info = ShyLUbasker->Symbolic(this->globalNumRows_,
220  this->globalNumCols_,
221  this->globalNumNonZeros_,
222  sp_rowptr,
223  sp_colind,
224  sp_values,
225  true);
226 
227  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
228  std::runtime_error, "Error in ShyLUBasker Symbolic");
229 
230  if (info == BASKER_SUCCESS) {
231  info = ShyLUbaskerTr->Symbolic(this->globalNumRows_,
232  this->globalNumCols_,
233  this->globalNumNonZeros_,
234  sp_rowptr,
235  sp_colind,
236  sp_values,
237  true);
238 
239  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
240  std::runtime_error, "Error in ShyLUBaskerTr Symbolic");
241  }
242  }
243  else
244  { //follow original code path if conditions not met
245  // In this case, loadA_impl updates colptr_, rowind_, nzvals_
246  info = ShyLUbasker->Symbolic(this->globalNumRows_,
247  this->globalNumCols_,
248  this->globalNumNonZeros_,
249  colptr_.getRawPtr(),
250  rowind_.getRawPtr(),
251  nzvals_.getRawPtr());
252 
253  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
254  std::runtime_error, "Error in ShyLUBasker Symbolic");
255 
256  if (info == BASKER_SUCCESS) {
257  info = ShyLUbaskerTr->Symbolic(this->globalNumRows_,
258  this->globalNumCols_,
259  this->globalNumNonZeros_,
260  colptr_.getRawPtr(),
261  rowind_.getRawPtr(),
262  nzvals_.getRawPtr(),
263  false);
264 
265  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
266  std::runtime_error, "Error in ShyLUBaskerTr Symbolic");
267  }
268  }
269  } // end if (this->root_)
270  /*No symbolic factoriztion*/
271 
272  /* All processes should have the same error code */
273  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
274  return(info);
275 }
276 
277 
278 template <class Matrix, class Vector>
279 int
281 {
282  using Teuchos::as;
283 
284  int info = 0;
285  if ( this->root_ ){
286  { // Do factorization
287 #ifdef HAVE_AMESOS2_TIMERS
288  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
289 #endif
290 
291 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
292  std::cout << "ShyLUBasker:: Before numeric factorization" << std::endl;
293  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
294  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
295  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
296 #endif
297 
298  // NDE: Special case
299  // Rather than going through the Amesos2 machinery to convert the matrixA_ CRS pointer data to CCS and store in Teuchos::Arrays,
300  // in this special case we pass the CRS raw pointers directly to ShyLUBasker which copies+transposes+sorts the data for CCS format
301  // loadA_impl is essentially an empty function in this case, as the raw pointers are handled here and similarly in Symbolic
302 
303  if ( single_proc_optimization() ) {
304 
305  auto sp_rowptr = this->matrixA_->returnRowPtr();
306  TEUCHOS_TEST_FOR_EXCEPTION(sp_rowptr == nullptr,
307  std::runtime_error, "Amesos2 Runtime Error: sp_rowptr returned null ");
308  auto sp_colind = this->matrixA_->returnColInd();
309  TEUCHOS_TEST_FOR_EXCEPTION(sp_colind == nullptr,
310  std::runtime_error, "Amesos2 Runtime Error: sp_colind returned null ");
311 #ifndef HAVE_TEUCHOS_COMPLEX
312  auto sp_values = this->matrixA_->returnValues();
313 #else
314  // NDE: 09/25/2017
315  // Cannot convert Kokkos::complex<T>* to std::complex<T>*; in this case, use reinterpret_cast
316  using complex_type = typename Util::getStdCplxType< magnitude_type, typename matrix_adapter_type::spmtx_vals_t >::type;
317  complex_type * sp_values = nullptr;
318  sp_values = reinterpret_cast< complex_type * > ( this->matrixA_->returnValues() );
319 #endif
320  TEUCHOS_TEST_FOR_EXCEPTION(sp_values == nullptr,
321  std::runtime_error, "Amesos2 Runtime Error: sp_values returned null ");
322 
323  // In this case, colptr_, rowind_, nzvals_ are invalid
324  info = ShyLUbasker->Factor( this->globalNumRows_,
325  this->globalNumCols_,
326  this->globalNumNonZeros_,
327  sp_rowptr,
328  sp_colind,
329  sp_values);
330 
331  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
332  std::runtime_error, "Error ShyLUBasker Factor");
333 
334  if (info == 0) {
335  info = ShyLUbaskerTr->Factor( this->globalNumRows_,
336  this->globalNumCols_,
337  this->globalNumNonZeros_,
338  sp_rowptr,
339  sp_colind,
340  sp_values);
341 
342  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
343  std::runtime_error, "Error ShyLUBaskerTr Factor");
344  }
345  }
346  else
347  {
348  // In this case, loadA_impl updates colptr_, rowind_, nzvals_
349  info = ShyLUbasker->Factor(this->globalNumRows_,
350  this->globalNumCols_,
351  this->globalNumNonZeros_,
352  colptr_.getRawPtr(),
353  rowind_.getRawPtr(),
354  nzvals_.getRawPtr());
355 
356  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
357  std::runtime_error, "Error ShyLUBasker Factor");
358 
359  if (info == 0) {
360  info = ShyLUbaskerTr->Factor(this->globalNumRows_,
361  this->globalNumCols_,
362  this->globalNumNonZeros_,
363  colptr_.getRawPtr(),
364  rowind_.getRawPtr(),
365  nzvals_.getRawPtr());
366 
367  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
368  std::runtime_error, "Error ShyLUBaskerTr Factor");
369  }
370  //We need to handle the realloc options
371  }
372 
373  //ShyLUbasker->DEBUG_PRINT();
374 
375  local_ordinal_type blnnz = local_ordinal_type(0);
376  local_ordinal_type bunnz = local_ordinal_type(0);
377  ShyLUbasker->GetLnnz(blnnz); // Add exception handling?
378  ShyLUbasker->GetUnnz(bunnz);
379 
380  local_ordinal_type Trblnnz = local_ordinal_type(0);
381  local_ordinal_type Trbunnz = local_ordinal_type(0);
382  ShyLUbaskerTr->GetLnnz(Trblnnz); // Add exception handling?
383  ShyLUbaskerTr->GetUnnz(Trbunnz);
384 
385  // This is set after numeric factorization complete as pivoting can be used;
386  // In this case, a discrepancy between symbolic and numeric nnz total can occur.
387  this->setNnzLU( as<size_t>( blnnz + bunnz ) );
388 
389  } // end scope for timer
390  } // end if (this->root_)
391 
392  /* All processes should have the same error code */
393  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
394 
395  //global_size_type info_st = as<global_size_type>(info);
396  /* TODO : Proper error messages*/
397  TEUCHOS_TEST_FOR_EXCEPTION(info == -1,
398  std::runtime_error,
399  "ShyLUBasker: Could not alloc space for L and U");
400  TEUCHOS_TEST_FOR_EXCEPTION(info == -2,
401  std::runtime_error,
402  "ShyLUBasker: Could not alloc needed work space");
403  TEUCHOS_TEST_FOR_EXCEPTION(info == -3,
404  std::runtime_error,
405  "ShyLUBasker: Could not alloc additional memory needed for L and U");
406  TEUCHOS_TEST_FOR_EXCEPTION(info > 0,
407  std::runtime_error,
408  "ShyLUBasker: Zero pivot found at: " << info );
409 
410  return(info);
411 }
412 
413 
414 template <class Matrix, class Vector>
415 int
417  const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
418  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
419 {
420  int ierr = 0; // returned error code
421 
422  using Teuchos::as;
423 
424  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
425  const size_t nrhs = X->getGlobalNumVectors();
426 
427  bool ShyluBaskerTransposeRequest = this->control_.useTranspose_;
428 
429  if ( single_proc_optimization() && nrhs == 1 ) {
430 
431 #ifdef HAVE_AMESOS2_TIMERS
432  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
433 #endif
434 
435 #ifndef HAVE_TEUCHOS_COMPLEX
436  auto b_vector = Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( B );
437  auto x_vector = Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( X );
438 #else
439  // NDE: 09/25/2017
440  // Cannot convert Kokkos::complex<T>* to std::complex<T>*; in this case, use reinterpret_cast
441  using complex_type = typename Util::getStdCplxType< magnitude_type, typename matrix_adapter_type::spmtx_vals_t >::type;
442  complex_type * b_vector = reinterpret_cast< complex_type * >( Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( B ) );
443  complex_type * x_vector = reinterpret_cast< complex_type * >( Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( X ) );
444 #endif
445  TEUCHOS_TEST_FOR_EXCEPTION(b_vector == nullptr,
446  std::runtime_error, "Amesos2 Runtime Error: b_vector returned null ");
447 
448  TEUCHOS_TEST_FOR_EXCEPTION(x_vector == nullptr,
449  std::runtime_error, "Amesos2 Runtime Error: x_vector returned null ");
450 
451  if ( this->root_ ) {
452  { // Do solve!
453 #ifdef HAVE_AMESOS2_TIMERS
454  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
455 #endif
456  if (!ShyluBaskerTransposeRequest)
457  ierr = ShyLUbasker->Solve(nrhs, b_vector, x_vector);
458  else
459  ierr = ShyLUbaskerTr->Solve(nrhs, b_vector, x_vector);
460  } // end scope for timer
461 
462  /* All processes should have the same error code */
463  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
464 
465  TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
466  std::runtime_error,
467  "Encountered zero diag element at: " << ierr);
468  TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
469  std::runtime_error,
470  "Could not alloc needed working memory for solve" );
471  } //end if (this->root_)
472  } // end if ( single_proc_optimization() && nrhs == 1 )
473  else
474  {
475  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
476 
477  xvals_.resize(val_store_size);
478  bvals_.resize(val_store_size);
479 
480  { // Get values from RHS B
481 #ifdef HAVE_AMESOS2_TIMERS
482  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
483  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
484 #endif
485 
486  if ( is_contiguous_ == true ) {
488  slu_type>::do_get(B, bvals_(), as<size_t>(ld_rhs), ROOTED, this->rowIndexBase_);
489  }
490  else {
492  slu_type>::do_get(B, bvals_(), as<size_t>(ld_rhs), CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
493  }
494  }
495 
496  if ( this->root_ ) {
497  { // Do solve!
498 #ifdef HAVE_AMESOS2_TIMERS
499  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
500 #endif
501  if (!ShyluBaskerTransposeRequest)
502  ierr = ShyLUbasker->Solve(nrhs, bvals_.getRawPtr(), xvals_.getRawPtr());
503  else
504  ierr = ShyLUbaskerTr->Solve(nrhs, bvals_.getRawPtr(), xvals_.getRawPtr());
505  } // end scope for timer
506  } // end if (this->root_)
507 
508  /* All processes should have the same error code */
509  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
510 
511  TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
512  std::runtime_error,
513  "Encountered zero diag element at: " << ierr);
514  TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
515  std::runtime_error,
516  "Could not alloc needed working memory for solve" );
517 
518  {
519 #ifdef HAVE_AMESOS2_TIMERS
520  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
521 #endif
522 
523  if ( is_contiguous_ == true ) {
525  MultiVecAdapter<Vector>,slu_type>::do_put(X, xvals_(),
526  as<size_t>(ld_rhs),
527  ROOTED);
528  }
529  else {
531  MultiVecAdapter<Vector>,slu_type>::do_put(X, xvals_(),
532  as<size_t>(ld_rhs),
534  }
535  } // end scope for timer
536  } // end else
537 
538  return(ierr);
539 }
540 
541 
542 template <class Matrix, class Vector>
543 bool
545 {
546  // The ShyLUBasker can only handle square for right now
547  return( this->globalNumRows_ == this->globalNumCols_ );
548 }
549 
550 
551 template <class Matrix, class Vector>
552 void
553 ShyLUBasker<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
554 {
555  using Teuchos::RCP;
556  using Teuchos::getIntegralValue;
557  using Teuchos::ParameterEntryValidator;
558 
559  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
560 
561  if(parameterList->isParameter("IsContiguous"))
562  {
563  is_contiguous_ = parameterList->get<bool>("IsContiguous");
564  }
565 
566  if(parameterList->isParameter("num_threads"))
567  {
568  num_threads = parameterList->get<int>("num_threads");
569  }
570  if(parameterList->isParameter("pivot"))
571  {
572  ShyLUbasker->Options.no_pivot = (!parameterList->get<bool>("pivot"));
573  ShyLUbaskerTr->Options.no_pivot = (!parameterList->get<bool>("pivot"));
574  }
575  if(parameterList->isParameter("delayed pivot"))
576  {
577  ShyLUbasker->Options.static_delayed_pivot = (parameterList->get<int>("delayed pivot"));
578  ShyLUbaskerTr->Options.static_delayed_pivot = (parameterList->get<int>("delayed pivot"));
579  }
580  if(parameterList->isParameter("pivot_tol"))
581  {
582  ShyLUbasker->Options.pivot_tol = parameterList->get<double>("pivot_tol");
583  ShyLUbaskerTr->Options.pivot_tol = parameterList->get<double>("pivot_tol");
584  }
585  if(parameterList->isParameter("symmetric"))
586  {
587  ShyLUbasker->Options.symmetric = parameterList->get<bool>("symmetric");
588  ShyLUbaskerTr->Options.symmetric = parameterList->get<bool>("symmetric");
589  }
590  if(parameterList->isParameter("realloc"))
591  {
592  ShyLUbasker->Options.realloc = parameterList->get<bool>("realloc");
593  ShyLUbaskerTr->Options.realloc = parameterList->get<bool>("realloc");
594  }
595  if(parameterList->isParameter("verbose"))
596  {
597  ShyLUbasker->Options.verbose = parameterList->get<bool>("verbose");
598  ShyLUbaskerTr->Options.verbose = parameterList->get<bool>("verbose");
599  }
600  if(parameterList->isParameter("verbose_matrix"))
601  {
602  ShyLUbasker->Options.verbose_matrix_out = parameterList->get<bool>("verbose_matrix");
603  ShyLUbaskerTr->Options.verbose_matrix_out = parameterList->get<bool>("verbose_matrix");
604  }
605  if(parameterList->isParameter("btf"))
606  {
607  ShyLUbasker->Options.btf = parameterList->get<bool>("btf");
608  ShyLUbaskerTr->Options.btf = parameterList->get<bool>("btf");
609  }
610  if(parameterList->isParameter("use_metis"))
611  {
612  ShyLUbasker->Options.use_metis = parameterList->get<bool>("use_metis");
613  ShyLUbaskerTr->Options.use_metis = parameterList->get<bool>("use_metis");
614  }
615  if(parameterList->isParameter("run_nd_on_leaves"))
616  {
617  ShyLUbasker->Options.run_nd_on_leaves = parameterList->get<bool>("run_nd_on_leaves");
618  ShyLUbaskerTr->Options.run_nd_on_leaves = parameterList->get<bool>("run_nd_on_leaves");
619  }
620  if(parameterList->isParameter("transpose"))
621  {
622  // NDE: set transpose vs non-transpose mode as bool; track separate shylubasker objects
623  const auto transpose = parameterList->get<bool>("transpose");
624  if (transpose == true)
625  this->control_.useTranspose_ = true;
626  //ShyLUbasker->Options.transpose = parameterList->get<bool>("transpose");
627  //ShyLUbaskerTr->Options.transpose = parameterList->get<bool>("transpose");
628  }
629  if(parameterList->isParameter("use_sequential_diag_facto"))
630  {
631  ShyLUbasker->Options.use_sequential_diag_facto = parameterList->get<bool>("use_sequential_diag_facto");
632  ShyLUbaskerTr->Options.use_sequential_diag_facto = parameterList->get<bool>("use_sequential_diag_facto");
633  }
634  if(parameterList->isParameter("user_fill"))
635  {
636  ShyLUbasker->Options.user_fill = parameterList->get<double>("user_fill");
637  ShyLUbaskerTr->Options.user_fill = parameterList->get<double>("user_fill");
638  }
639  if(parameterList->isParameter("prune"))
640  {
641  ShyLUbasker->Options.prune = parameterList->get<bool>("prune");
642  ShyLUbaskerTr->Options.prune = parameterList->get<bool>("prune");
643  }
644  if(parameterList->isParameter("replace_tiny_pivot"))
645  {
646  ShyLUbasker->Options.prune = parameterList->get<bool>("replace_tiny_pivot");
647  ShyLUbaskerTr->Options.prune = parameterList->get<bool>("replace_tiny_pivot");
648  }
649  if(parameterList->isParameter("btf_matching"))
650  {
651  ShyLUbasker->Options.btf_matching = parameterList->get<int>("btf_matching");
652  ShyLUbaskerTr->Options.btf_matching = parameterList->get<int>("btf_matching");
653  if (ShyLUbasker->Options.btf_matching == 1 || ShyLUbasker->Options.btf_matching == 2) {
654  ShyLUbasker->Options.matching = true;
655  ShyLUbaskerTr->Options.matching = true;
656  } else {
657  ShyLUbasker->Options.matching = false;
658  ShyLUbaskerTr->Options.matching = false;
659  }
660  }
661  if(parameterList->isParameter("blk_matching"))
662  {
663  ShyLUbasker->Options.blk_matching = parameterList->get<int>("blk_matching");
664  ShyLUbaskerTr->Options.blk_matching = parameterList->get<int>("blk_matching");
665  }
666  if(parameterList->isParameter("min_block_size"))
667  {
668  ShyLUbasker->Options.min_block_size = parameterList->get<int>("min_block_size");
669  ShyLUbaskerTr->Options.min_block_size = parameterList->get<int>("min_block_size");
670  }
671 }
672 
673 template <class Matrix, class Vector>
674 Teuchos::RCP<const Teuchos::ParameterList>
676 {
677  using Teuchos::ParameterList;
678 
679  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
680 
681  if( is_null(valid_params) )
682  {
683  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
684  pl->set("IsContiguous", true,
685  "Are GIDs contiguous");
686  pl->set("num_threads", 1,
687  "Number of threads");
688  pl->set("pivot", false,
689  "Should not pivot");
690  pl->set("delayed pivot", 0,
691  "Apply static delayed pivot on a big block");
692  pl->set("pivot_tol", .0001,
693  "Tolerance before pivot, currently not used");
694  pl->set("symmetric", false,
695  "Should Symbolic assume symmetric nonzero pattern");
696  pl->set("realloc" , false,
697  "Should realloc space if not enough");
698  pl->set("verbose", false,
699  "Information about factoring");
700  pl->set("verbose_matrix", false,
701  "Give Permuted Matrices");
702  pl->set("btf", true,
703  "Use BTF ordering");
704  pl->set("prune", false,
705  "Use prune on BTF blocks (Not Supported)");
706  pl->set("btf_matching", 2,
707  "Matching option for BTF: 0 = none, 1 = Basker, 2 = Trilinos (default), (3 = MC64 if enabled)");
708  pl->set("blk_matching", 1,
709  "Matching optioon for block: 0 = none, 1 or anything else = Basker (default), (2 = MC64 if enabled)");
710  pl->set("min_block_size", 0,
711  "Size of the minimum diagonal blocks");
712  pl->set("replace_tiny_pivot", true,
713  "Replace tiny pivots during the numerical factorization");
714  pl->set("use_metis", true,
715  "Use METIS for ND");
716  pl->set("run_nd_on_leaves", false,
717  "Run ND on the final leaf-nodes");
718  pl->set("transpose", false,
719  "Solve the transpose A");
720  pl->set("use_sequential_diag_facto", false,
721  "Use sequential algorithm to factor each diagonal block");
722  pl->set("user_fill", (double)BASKER_FILL_USER,
723  "User-provided padding for the fill ratio");
724  valid_params = pl;
725  }
726  return valid_params;
727 }
728 
729 
730 template <class Matrix, class Vector>
731 bool
733 {
734  using Teuchos::as;
735  if(current_phase == SOLVE) return (false);
736 
737  #ifdef HAVE_AMESOS2_TIMERS
738  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
739  #endif
740 
741 
742  // NDE: Can clean up duplicated code with the #ifdef guards
743  if ( single_proc_optimization() ) {
744  // NDE: Nothing is done in this special case - CRS raw pointers are passed to SHYLUBASKER and transpose of copies handled there
745  // In this case, colptr_, rowind_, nzvals_ are invalid
746  }
747  else
748  {
749 
750  // Only the root image needs storage allocated
751  if( this->root_ ){
752  nzvals_.resize(this->globalNumNonZeros_);
753  rowind_.resize(this->globalNumNonZeros_);
754  colptr_.resize(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...
755  }
756 
757  local_ordinal_type nnz_ret = 0;
758  {
759  #ifdef HAVE_AMESOS2_TIMERS
760  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
761  #endif
762 
763  if ( is_contiguous_ == true ) {
764  Util::get_ccs_helper<
765  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
766  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
767  nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_); // copies from matrixA_ to ShyLUBasker ConcreteSolver cp, ri, nzval members
768  }
769  else {
770  Util::get_ccs_helper<
771  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
772  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
773  nnz_ret, CONTIGUOUS_AND_ROOTED, ARBITRARY, this->rowIndexBase_); // copies from matrixA_ to ShyLUBasker ConcreteSolver cp, ri, nzval members
774  }
775  }
776 
777  if( this->root_ ){
778  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
779  std::runtime_error,
780  "Amesos2_ShyLUBasker loadA_impl: Did not get the expected number of non-zero vals");
781  }
782 
783  } //end alternative path
784  return true;
785 }
786 
787 
788 template<class Matrix, class Vector>
789 const char* ShyLUBasker<Matrix,Vector>::name = "ShyLUBasker";
790 
791 
792 } // end namespace Amesos2
793 
794 #endif // AMESOS2_SHYLUBASKER_DEF_HPP
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:156
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
Amesos2 interface to the Baker package.
Definition: Amesos2_ShyLUBasker_decl.hpp:76
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:266
Definition: Amesos2_TypeDecl.hpp:143
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_ShyLUBasker_def.hpp:162
void transpose(ArrayView< Scalar > vals, ArrayView< GlobalOrdinal > indices, ArrayView< GlobalSizeT > ptr, ArrayView< Scalar > trans_vals, ArrayView< GlobalOrdinal > trans_indices, ArrayView< GlobalSizeT > trans_ptr)
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
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
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:416
Amesos2 ShyLUBasker declarations.
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_ShyLUBasker_def.hpp:732
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
int numericFactorization_impl()
ShyLUBasker specific numeric factorization.
Definition: Amesos2_ShyLUBasker_def.hpp:280
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_ShyLUBasker_def.hpp:544
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_ShyLUBasker_def.hpp:675
Definition: Amesos2_TypeDecl.hpp:127
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:372
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
Definition: Amesos2_TypeDecl.hpp:128