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  , is_contiguous_(true)
76 // , basker()
77 {
78  //Nothing
79 }
80 
81 
82 template <class Matrix, class Vector>
83 Basker<Matrix,Vector>::~Basker( )
84 {}
85 
86 template <class Matrix, class Vector>
87 bool
89  return (this->root_ && (this->matrixA_->getComm()->getSize() == 1) && is_contiguous_);
90 }
91 
92 template<class Matrix, class Vector>
93 int
95 {
96  /* TODO: Define what it means for Basker
97  */
98 #ifdef HAVE_AMESOS2_TIMERS
99  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
100 #endif
101 
102  return(0);
103 }
104 
105 
106 template <class Matrix, class Vector>
107 int
109 {
110  /*No symbolic factoriztion*/
111  return(0);
112 }
113 
114 
115 template <class Matrix, class Vector>
116 int
118 {
119  using Teuchos::as;
120 
121  int info = 0;
122  if ( this->root_ ){
123  { // Do factorization
124  #ifdef HAVE_AMESOS2_TIMERS
125  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
126  #endif
127 
128  #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
129  std::cout << "Basker:: Before numeric factorization" << std::endl;
130  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
131  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
132  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
133  #endif
134 
135  info = basker.factor(this->globalNumRows_, this->globalNumCols_, this->globalNumNonZeros_, colptr_.getRawPtr(), rowind_.getRawPtr(), nzvals_.getRawPtr());
136 
137  // This is set after numeric factorization complete as pivoting can be used;
138  // In this case, a discrepancy between symbolic and numeric nnz total can occur.
139  this->setNnzLU( as<size_t>(basker.get_NnzLU() ) ) ;
140  }
141 
142  }
143 
144  /* All processes should have the same error code */
145  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
146 
147  //global_size_type info_st = as<global_size_type>(info);
148  /* TODO : Proper error messages*/
149  TEUCHOS_TEST_FOR_EXCEPTION( (info == -1) ,
150  std::runtime_error,
151  "Basker: Could not alloc space for L and U");
152  TEUCHOS_TEST_FOR_EXCEPTION( (info == -2),
153  std::runtime_error,
154  "Basker: Could not alloc needed work space");
155  TEUCHOS_TEST_FOR_EXCEPTION( (info == -3) ,
156  std::runtime_error,
157  "Basker: Could not alloc additional memory needed for L and U");
158  TEUCHOS_TEST_FOR_EXCEPTION( (info > 0) ,
159  std::runtime_error,
160  "Basker: Zero pivot found at: " << info );
161 
162  return(info);
163 }
164 
165 
166 template <class Matrix, class Vector>
167 int
169  const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
170  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
171 {
172  int ierr = 0; // returned error code
173 
174  using Teuchos::as;
175 
176  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
177  const size_t nrhs = X->getGlobalNumVectors();
178 
179  if ( single_proc_optimization() && nrhs == 1 ) {
180 
181 #ifdef HAVE_AMESOS2_TIMERS
182  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
183 #endif
184 
185 #ifndef HAVE_TEUCHOS_COMPLEX
186  auto b_vector = Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( B );
187  auto x_vector = Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( X );
188 #else
189  // NDE: 09/25/2017
190  // Cannot convert Kokkos::complex<T>* to std::complex<T>*; in this case, use reinterpret_cast
191  using complex_type = typename Util::getStdCplxType< magnitude_type, typename matrix_adapter_type::spmtx_vals_t >::type;
192  complex_type * b_vector = reinterpret_cast< complex_type * >( Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( B ) );
193  complex_type * x_vector = reinterpret_cast< complex_type * >( Util::vector_pointer_helper< MultiVecAdapter<Vector>, Vector >::get_pointer_to_vector( X ) );
194 #endif
195  TEUCHOS_TEST_FOR_EXCEPTION(b_vector == nullptr,
196  std::runtime_error, "Amesos2 Runtime Error: b_vector returned null ");
197 
198  TEUCHOS_TEST_FOR_EXCEPTION(x_vector == nullptr,
199  std::runtime_error, "Amesos2 Runtime Error: x_vector returned null ");
200 
201  if ( this->root_ ) {
202  { // Do solve!
203 #ifdef HAVE_AMESOS2_TIMERS
204  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
205 #endif
206  ierr = basker.solveMultiple(nrhs, b_vector, x_vector);
207  }
208 
209  /* All processes should have the same error code */
210  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
211 
212  TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
213  std::runtime_error,
214  "Encountered zero diag element at: " << ierr);
215  TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
216  std::runtime_error,
217  "Could not alloc needed working memory for solve" );
218  }
219  }
220  else
221  {
222  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
223 
224  xvals_.resize(val_store_size);
225  bvals_.resize(val_store_size);
226 
227  { // Get values from RHS B
228 #ifdef HAVE_AMESOS2_TIMERS
229  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
230  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
231 #endif
232 
233  if ( is_contiguous_ == true ) {
235  slu_type>::do_get(B, bvals_(), as<size_t>(ld_rhs), ROOTED, this->rowIndexBase_);
236  }
237  else {
239  slu_type>::do_get(B, bvals_(), as<size_t>(ld_rhs), CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
240  }
241  }
242 
243  if ( this->root_ ) {
244  { // Do solve!
245 #ifdef HAVE_AMESOS2_TIMERS
246  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
247 #endif
248 
249  ierr = basker.solveMultiple(nrhs, bvals_.getRawPtr(),xvals_.getRawPtr());
250  }
251 
252  }
253 
254  /* All processes should have the same error code */
255  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
256 
257  TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
258  std::runtime_error,
259  "Encountered zero diag element at: " << ierr);
260  TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
261  std::runtime_error,
262  "Could not alloc needed working memory for solve" );
263 
264  {
265 #ifdef HAVE_AMESOS2_TIMERS
266  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
267 #endif
268 
269  if ( is_contiguous_ == true ) {
271  MultiVecAdapter<Vector>,slu_type>::do_put(X, xvals_(),
272  as<size_t>(ld_rhs),
273  ROOTED);
274  }
275  else {
277  MultiVecAdapter<Vector>,slu_type>::do_put(X, xvals_(),
278  as<size_t>(ld_rhs),
279  CONTIGUOUS_AND_ROOTED);
280  }
281  }
282  }
283 
284  return(ierr);
285 }
286 
287 
288 template <class Matrix, class Vector>
289 bool
291 {
292  // The Basker can only handle square for right now
293  return( this->globalNumRows_ == this->globalNumCols_ );
294 }
295 
296 
297 template <class Matrix, class Vector>
298 void
299 Basker<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
300 {
301  using Teuchos::RCP;
302  using Teuchos::getIntegralValue;
303  using Teuchos::ParameterEntryValidator;
304 
305  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
306 
307  if(parameterList->isParameter("IsContiguous"))
308  {
309  is_contiguous_ = parameterList->get<bool>("IsContiguous");
310  }
311 
312 }
313 
314 template <class Matrix, class Vector>
315 Teuchos::RCP<const Teuchos::ParameterList>
317 {
318  using Teuchos::ParameterList;
319 
320  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
321 
322  if( is_null(valid_params) ){
323  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
324 
325  pl->set("IsContiguous", true, "Are GIDs contiguous");
326  pl->set("alnnz", 2, "Approx number of nonzeros in L, default is 2*nnz(A)");
327  pl->set("aunnx", 2, "Approx number of nonzeros in I, default is 2*nnz(U)");
328  valid_params = pl;
329  }
330  return valid_params;
331 }
332 
333 
334 template <class Matrix, class Vector>
335 bool
337 {
338  using Teuchos::as;
339  if(current_phase == SOLVE) return (false);
340 
341  #ifdef HAVE_AMESOS2_TIMERS
342  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
343  #endif
344 
345  // Only the root image needs storage allocated
346  if( this->root_ ){
347  nzvals_.resize(this->globalNumNonZeros_);
348  rowind_.resize(this->globalNumNonZeros_);
349  colptr_.resize(this->globalNumCols_ + 1);
350  }
351 
352  local_ordinal_type nnz_ret = 0;
353  {
354  #ifdef HAVE_AMESOS2_TIMERS
355  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
356  #endif
357 
358  if ( is_contiguous_ == true ) {
360  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
361  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
362  nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_);
363  }
364  else {
366  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
367  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
368  nnz_ret, CONTIGUOUS_AND_ROOTED, ARBITRARY, this->rowIndexBase_);
369  }
370  }
371 
372  if( this->root_ ){
373  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
374  std::runtime_error,
375  "Amesos2_Basker loadA_impl: Did not get the expected number of non-zero vals");
376  }
377  return true;
378 }
379 
380 
381 template<class Matrix, class Vector>
382 const char* Basker<Matrix,Vector>::name = "Basker";
383 
384 
385 } // end namespace Amesos2
386 
387 #endif // AMESOS2_Basker_DEF_HPP
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Basker_def.hpp:336
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:266
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Basker_def.hpp:316
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
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:654
Amesos2 interface to the Baker package.
Definition: Amesos2_Basker_decl.hpp:73
Amesos2 Basker declarations.
bool single_proc_optimization() const
can we optimize size_type and ordinal_type for straight pass through, also check that is_contiguous_ ...
Definition: Amesos2_Basker_def.hpp:88
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Basker_def.hpp:290
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Basker_def.hpp:94
int numericFactorization_impl()
Basker specific numeric factorization.
Definition: Amesos2_Basker_def.hpp:117
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:322
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
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:168