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