Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Basker_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 
53 #ifndef AMESOS2_BASKER_DEF_HPP
54 #define AMESOS2_BASKER_DEF_HPP
55 
56 #include <Teuchos_Tuple.hpp>
57 #include <Teuchos_ParameterList.hpp>
58 #include <Teuchos_StandardParameterEntryValidators.hpp>
59 
61 #include "Amesos2_Basker_decl.hpp"
62 
63 namespace Amesos2 {
64 
65 
66 template <class Matrix, class Vector>
67 Basker<Matrix,Vector>::Basker(
68  Teuchos::RCP<const Matrix> A,
69  Teuchos::RCP<Vector> X,
70  Teuchos::RCP<const Vector> B )
71  : SolverCore<Amesos2::Basker,Matrix,Vector>(A, X, B)
72  , nzvals_() // initialize to empty arrays
73  , rowind_()
74  , colptr_()
75 {
76 
77  //Nothing
78 
79  // Override some default options
80  // TODO: use data_ here to init
81 
82 
83 
84 #ifdef SHYLUBASKER
85 #ifdef HAVE_AMESOS2_KOKKOS
86 #ifdef KOKKOS_HAVE_OPENMP
87  /*
88  static_assert(std::is_same<kokkos_exe,Kokkos::OpenMP>::value,
89  "Kokkos node type not supported by experimental Basker Amesos2");
90  */
91  typedef Kokkos::OpenMP Exe_Space;
92 #elif defined(KOKKOS_HAVE_SERIAL)
93  typedef Kokkos::Serial Exe_Space;
94 #else
95  TEUCHOS_TEST_FOR_EXCEPTION(1 != 0,
96  std::runtime_error,
97  "Do not have supported Kokkos node type for Basker");
98 #endif
99  //std::cout << "MAKE BASKER" << std::endl;
100  basker = new ::BaskerNS::Basker<local_ordinal_type, slu_type, Exe_Space>();
101  basker->Options.no_pivot = BASKER_TRUE;
102  basker->Options.symmetric = BASKER_FALSE;
103  basker->Options.realloc = BASKER_FALSE;
104  basker->Options.verbose = BASKER_FALSE;
105  basker->Options.matching = BASKER_TRUE;
106  basker->Options.matching_type = 0;
107  basker->Options.btf = BASKER_TRUE;
108  basker->Options.amd_btf = BASKER_TRUE;
109  basker->Options.amd_dom = BASKER_TRUE;
110  basker->Options.transpose = BASKER_FALSE;
111  basker->Options.verbose_matrix_out = BASKER_FALSE;
112  num_threads = Kokkos::OpenMP::max_hardware_threads();
113 #endif
114 #endif
115 
116 }
117 
118 
119 template <class Matrix, class Vector>
120 Basker<Matrix,Vector>::~Basker( )
121 {
122 #ifdef SHYLUBASKER
123 #ifdef HAVE_AMESOS2_KOKKOS
124  //std::cout << "DELETE BASKER" << std::endl;
125  delete basker;
126 #endif
127 #endif
128 
129  /* Basker will cleanup its own internal memory*/
130 }
131 
132 template<class Matrix, class Vector>
133 int
135 {
136  /* TODO: Define what it means for Basker
137  */
138 #ifdef HAVE_AMESOS2_TIMERS
139  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
140 #endif
141 
142  return(0);
143 }
144 
145 
146 template <class Matrix, class Vector>
147 int
149 {
150 
151 #ifdef SHYLUBASKER
152  if(this->root_)
153  {
154  //std::cout << "setting number of threads "
155  // << num_threads << std::endl;
156  basker->SetThreads(num_threads);
157  //std::cout << "Set Threads Done" << std::endl;
158 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
159  std::cout << "Basker:: Before symbolic factorization" << std::endl;
160  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
161  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
162  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
163 #endif
164  int info;
165  info =basker->Symbolic(this->globalNumRows_,
166  this->globalNumCols_,
167  this->globalNumNonZeros_,
168  colptr_.getRawPtr(),
169  rowind_.getRawPtr(),
170  nzvals_.getRawPtr());
171  //std::cout << "Symbolic Factorization Done" << std::endl;
172  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
173  std::runtime_error,
174  "Error in Basker Symbolic");
175 
176  }
177 #endif
178 
179  /*No symbolic factoriztion*/
180  return(0);
181 }
182 
183 
184 template <class Matrix, class Vector>
185 int
187 {
188  using Teuchos::as;
189 
190  int info = 0;
191  if ( this->root_ ){
192  { // Do factorization
193 #ifdef HAVE_AMESOS2_TIMERS
194  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
195 #endif
196 
197  #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
198  std::cout << "Basker:: Before numeric factorization" << std::endl;
199  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
200  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
201  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
202  #endif
203 
204 
205 #ifdef SHYLUBASKER
206  info = basker->Factor(this->globalNumRows_,
207  this->globalNumCols_,
208  this->globalNumNonZeros_,
209  colptr_.getRawPtr(),
210  rowind_.getRawPtr(),
211  nzvals_.getRawPtr());
212  //We need to handle the realloc options
213 
214  //basker->DEBUG_PRINT();
215 
216  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
217  std::runtime_error,
218  "Error Basker Factor");
219 
220 #else
221 
222  info =basker.factor(this->globalNumRows_, this->globalNumCols_, this->globalNumNonZeros_, colptr_.getRawPtr(), rowind_.getRawPtr(), nzvals_.getRawPtr());
223  #endif
224 
225  }
226 
227  }
228 
229  /* All processes should have the same error code */
230  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
231 
232  //global_size_type info_st = as<global_size_type>(info);
233  /* TODO : Proper error messages*/
234  TEUCHOS_TEST_FOR_EXCEPTION( (info == -1) ,
235  std::runtime_error,
236  "Basker: Could not alloc space for L and U");
237  TEUCHOS_TEST_FOR_EXCEPTION( (info == -2),
238  std::runtime_error,
239  "Basker: Could not alloc needed work space");
240  TEUCHOS_TEST_FOR_EXCEPTION( (info == -3) ,
241  std::runtime_error,
242  "Basker: Could not alloc additional memory needed for L and U");
243  TEUCHOS_TEST_FOR_EXCEPTION( (info > 0) ,
244  std::runtime_error,
245  "Basker: Zero pivot found at: " << info );
246 
247  return(info);
248 }
249 
250 
251 template <class Matrix, class Vector>
252 int
254  const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
255  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
256 {
257  using Teuchos::as;
258 
259  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
260  const size_t nrhs = X->getGlobalNumVectors();
261 
262  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
263 
264  xvals_.resize(val_store_size);
265  bvals_.resize(val_store_size);
266 
267  { // Get values from RHS B
268 #ifdef HAVE_AMESOS2_TIMERS
269  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
270  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
271 #endif
273  slu_type>::do_get(B, bvals_(),as<size_t>(ld_rhs),
274  ROOTED);
275  }
276 
277  int ierr = 0; // returned error code
278 
279  if ( this->root_ ) {
280  { // Do solve!
281 #ifdef HAVE_AMESOS2_TIMERS
282  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
283 #endif
284 
285 #ifdef SHYLUBASKER
286  ierr = basker->Solve(nrhs, bvals_.getRawPtr(),
287  xvals_.getRawPtr());
288 #else
289  ierr = basker.solveMultiple(nrhs, bvals_.getRawPtr(),xvals_.getRawPtr());
290 #endif
291  }
292 
293  }
294 
295  /* All processes should have the same error code */
296  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
297 
298  TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
299  std::runtime_error,
300  "Encountered zero diag element at: " << ierr);
301  TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
302  std::runtime_error,
303  "Could not alloc needed working memory for solve" );
304 
305  {
306 #ifdef HAVE_AMESOS2_TIMERS
307  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
308 #endif
309 
311  MultiVecAdapter<Vector>,slu_type>::do_put(X, xvals_(),
312  as<size_t>(ld_rhs),
313  ROOTED);
314  }
315 
316  return(ierr);
317 }
318 
319 
320 template <class Matrix, class Vector>
321 bool
323 {
324  // The Basker can only handle square for right now
325  return( this->globalNumRows_ == this->globalNumCols_ );
326 }
327 
328 
329 template <class Matrix, class Vector>
330 void
331 Basker<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
332 {
333  using Teuchos::RCP;
334  using Teuchos::getIntegralValue;
335  using Teuchos::ParameterEntryValidator;
336 
337  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
338 
339 #ifdef SHYLUBASKER
340  if(parameterList->isParameter("num_threads"))
341  {
342  num_threads = parameterList->get<int>("num_threads");
343  }
344  if(parameterList->isParameter("pivot"))
345  {
346  basker->Options.no_pivot = (!parameterList->get<bool>("pivot"));
347  }
348  if(parameterList->isParameter("pivot_tol"))
349  {
350  basker->Options.pivot_tol = parameterList->get<double>("pivot_tol");
351  }
352  if(parameterList->isParameter("symmetric"))
353  {
354  basker->Options.symmetric = parameterList->get<bool>("symmetric");
355  }
356  if(parameterList->isParameter("realloc"))
357  {
358  basker->Options.realloc = parameterList->get<bool>("realloc");
359  }
360  if(parameterList->isParameter("verbose"))
361  {
362  basker->Options.verbose = parameterList->get<bool>("verbose");
363  }
364  if(parameterList->isParameter("verbose_matrix"))
365  {
366  basker->Options.verbose_matrix_out = parameterList->get<bool>("verbose_matrix");
367  }
368  if(parameterList->isParameter("matching"))
369  {
370  basker->Options.matching = parameterList->get<bool>("matching");
371  }
372  if(parameterList->isParameter("matching_type"))
373  {
374  basker->Options.matching_type =
375  (local_ordinal_type) parameterList->get<int>("matching_type");
376  }
377  if(parameterList->isParameter("btf"))
378  {
379  basker->Options.btf = parameterList->get<bool>("btf");
380  }
381  if(parameterList->isParameter("amd_btf"))
382  {
383  basker->Options.amd_btf = parameterList->get<bool>("amd_btf");
384  }
385  if(parameterList->isParameter("amd_dom"))
386  {
387  basker->Options.amd_dom = parameterList->get<bool>("amd_dom");
388  }
389  if(parameterList->isParameter("transpose"))
390  {
391  basker->Options.transpose = parameterList->get<bool>("transpose");
392  }
393 #endif
394 
395 }
396 
397 template <class Matrix, class Vector>
398 Teuchos::RCP<const Teuchos::ParameterList>
400 {
401  using Teuchos::ParameterList;
402 
403  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
404 
405 
406 #ifdef SHYLUBASKER
407  if( is_null(valid_params) )
408  {
409  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
410  pl->set("num_threads", 1,
411  "Number of threads");
412  pl->set("pivot", false,
413  "Should not pivot");
414  pl->set("pivot_tol", .0001,
415  "Tolerance before pivot, currently not used");
416  pl->set("symmetric", false,
417  "Should Symbolic assume symmetric nonzero pattern");
418  pl->set("realloc" , false,
419  "Should realloc space if not enough");
420  pl->set("verbose", false,
421  "Information about factoring");
422  pl->set("verbose_matrix", false,
423  "Give Permuted Matrices");
424  pl->set("matching", true,
425  "Use WC matching (Not Supported)");
426  pl->set("matching_type", 0,
427  "Type of WC matching (Not Supported)");
428  pl->set("btf", true,
429  "Use BTF ordering");
430  pl->set("amd_btf", true,
431  "Use AMD on BTF blocks (Not Supported)");
432  pl->set("amd_dom", true,
433  "Use CAMD on ND blocks (Not Supported)");
434  pl->set("transpose", false,
435  "Solve the transpose A");
436  valid_params = pl;
437  }
438  return valid_params;
439 #else
440  if( is_null(valid_params) ){
441  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
442 
443  pl->set("alnnz", 2, "Approx number of nonzeros in L, default is 2*nnz(A)");
444  pl->set("aunnx", 2, "Approx number of nonzeros in I, default is 2*nnz(U)");
445  valid_params = pl;
446  }
447 #endif
448  return valid_params;
449 }
450 
451 
452 template <class Matrix, class Vector>
453 bool
455 {
456  using Teuchos::as;
457 
458  if(current_phase == SOLVE) return (false);
459 
460 #ifdef HAVE_AMESOS2_TIMERS
461  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
462 #endif
463 
464 
465 
466  // Only the root image needs storage allocated
467  if( this->root_ ){
468  nzvals_.resize(this->globalNumNonZeros_);
469  rowind_.resize(this->globalNumNonZeros_);
470  colptr_.resize(this->globalNumCols_ + 1);
471  }
472 
473  local_ordinal_type nnz_ret = 0;
474  {
475 #ifdef HAVE_AMESOS2_TIMERS
476  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
477 #endif
478 
480  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
481  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
482  nnz_ret, ROOTED, ARBITRARY);
483  }
484 
485 
486  if( this->root_ ){
487  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
488  std::runtime_error,
489  "Did not get the expected number of non-zero vals");
490  }
491 
492  return true;
493 }
494 
495 
496 template<class Matrix, class Vector>
497 const char* Basker<Matrix,Vector>::name = "Basker";
498 
499 
500 } // end namespace Amesos2
501 
502 #endif // AMESOS2_Basker_DEF_HPP
Definition: basker.cpp:35
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Basker_def.hpp:399
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Basker_def.hpp:454
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
Definition: Amesos2_TypeDecl.hpp:142
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
Basker specific solve.
Definition: Amesos2_Basker_def.hpp:253
Amesos2 interface to the Baker package.
Definition: Amesos2_Basker_decl.hpp:74
Amesos2 Basker declarations.
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Basker_def.hpp:134
Definition: Amesos2_TypeDecl.hpp:127
int numericFactorization_impl()
Basker specific numeric factorization.
Definition: Amesos2_Basker_def.hpp:186
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:296
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:175
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Basker_def.hpp:322