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 #ifdef SHYLUBASKER
64 #include "Kokkos_Core.hpp"
65 #endif
66 
67 namespace Amesos2 {
68 
69 
70 template <class Matrix, class Vector>
71 Basker<Matrix,Vector>::Basker(
72  Teuchos::RCP<const Matrix> A,
73  Teuchos::RCP<Vector> X,
74  Teuchos::RCP<const Vector> B )
75  : SolverCore<Amesos2::Basker,Matrix,Vector>(A, X, B)
76  , nzvals_() // initialize to empty arrays
77  , rowind_()
78  , colptr_()
79 {
80 
81  //Nothing
82 
83  // Override some default options
84  // TODO: use data_ here to init
85 
86 
87 
88 #ifdef SHYLUBASKER
89 #ifdef HAVE_AMESOS2_KOKKOS
90 
91  typedef Kokkos::OpenMP Exe_Space;
92 
93  basker = new ::BaskerNS::Basker<local_ordinal_type, slu_type, Exe_Space>();
94 
95  printf("Constructor Called \n");
96 
97  basker->Options.no_pivot = true;
98  basker->Options.symmetric = false;
99  basker->Options.realloc = false;
100  basker->Options.verbose = false;
101  basker->Options.btf = true;
102 
103  num_threads = Kokkos::OpenMP::max_hardware_threads();
104 
105 
106 #endif
107 #endif
108 
109 }
110 
111 
112 template <class Matrix, class Vector>
113 Basker<Matrix,Vector>::~Basker( )
114 {
115  printf("Destructor Call \n");
116 #ifdef SHYLUBASKER
117 #ifdef HAVE_AMESOS2_KOKKOS
118 
119  delete basker;
120 #endif
121 #endif
122  /* Basker will cleanup its own internal memory*/
123 }
124 
125 template<class Matrix, class Vector>
126 int
128 {
129  /* TODO: Define what it means for Basker
130  */
131 #ifdef HAVE_AMESOS2_TIMERS
132  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
133 #endif
134 
135  return(0);
136 }
137 
138 
139 template <class Matrix, class Vector>
140 int
142 {
143 
144 #ifdef SHYLUBASKER
145  //std::cout << "shylubasker sfactor" << std::endl;
146  if(this->root_)
147  {
148 
149 
150  basker->SetThreads(num_threads);
151 
152 
153  //std::cout << "Set Threads Done" << std::endl;
154 
155 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
156  std::cout << "Basker:: Before symbolic factorization" << std::endl;
157  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
158  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
159  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
160 #endif
161 
162 
163  int info;
164  info =basker->Symbolic(this->globalNumRows_,
165  this->globalNumCols_,
166  this->globalNumNonZeros_,
167  colptr_.getRawPtr(),
168  rowind_.getRawPtr(),
169  nzvals_.getRawPtr());
170 
171 
172  //std::cout << "Symbolic Factorization Done" << std::endl;
173 
174  }
175 
176 #endif
177  /*No symbolic factoriztion*/
178  return(0);
179 }
180 
181 
182 template <class Matrix, class Vector>
183 int
185 {
186  using Teuchos::as;
187 
188  int info = 0;
189  if ( this->root_ ){
190  { // Do factorization
191 #ifdef HAVE_AMESOS2_TIMERS
192  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
193 #endif
194 
195 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
196  std::cout << "Basker:: Before numeric factorization" << std::endl;
197  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
198  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
199  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
200 #endif
201 
202 #ifdef SHYLUBASKER
203  std::cout << "SHYLUBASKER FACTOR " << std::endl;
204 
205 
206  info = basker->Factor(this->globalNumRows_,
207  this->globalNumCols_,
208  this->globalNumNonZeros_,
209  colptr_.getRawPtr(),
210  rowind_.getRawPtr(),
211  nzvals_.getRawPtr());
212 
213 
214 #else
215  info =basker.factor(this->globalNumRows_, this->globalNumCols_, this->globalNumNonZeros_, colptr_.getRawPtr(), rowind_.getRawPtr(), nzvals_.getRawPtr());
216 #endif
217 
218  }
219 
220  }
221 
222  /* All processes should have the same error code */
223  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
224 
225  //global_size_type info_st = as<global_size_type>(info);
226  /* TODO : Proper error messages*/
227  TEUCHOS_TEST_FOR_EXCEPTION( (info == -1) ,
228  std::runtime_error,
229  "Basker: Could not alloc space for L and U");
230  TEUCHOS_TEST_FOR_EXCEPTION( (info == -2),
231  std::runtime_error,
232  "Basker: Could not alloc needed work space");
233  TEUCHOS_TEST_FOR_EXCEPTION( (info == -3) ,
234  std::runtime_error,
235  "Basker: Could not alloc additional memory needed for L and U");
236  TEUCHOS_TEST_FOR_EXCEPTION( (info > 0) ,
237  std::runtime_error,
238  "Basker: Zero pivot found at: " << info );
239 
240  return(info);
241 }
242 
243 
244 template <class Matrix, class Vector>
245 int
247  const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
248  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
249 {
250  using Teuchos::as;
251 
252  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
253  const size_t nrhs = X->getGlobalNumVectors();
254 
255  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
256 
257  xvals_.resize(val_store_size);
258  bvals_.resize(val_store_size);
259 
260  { // Get values from RHS B
261 #ifdef HAVE_AMESOS2_TIMERS
262  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
263  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
264 #endif
266  slu_type>::do_get(B, bvals_(),as<size_t>(ld_rhs),
267  ROOTED);
268  }
269 
270  int ierr = 0; // returned error code
271 
272  if ( this->root_ ) {
273  { // Do solve!
274 #ifdef HAVE_AMESOS2_TIMERS
275  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
276 #endif
277 
278 #ifdef SHYLUBASKER
279 
280  std::cout << "SHYLUBASKER Only handles 1 solve right now" << std::endl;
281  ierr = basker->Solve(nrhs, bvals_.getRawPtr(), xvals_.getRawPtr());
282 
283 #else
284  ierr = basker.solveMultiple(nrhs, bvals_.getRawPtr(),xvals_.getRawPtr());
285 #endif
286  }
287 
288  }
289 
290  /* All processes should have the same error code */
291  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
292 
293  TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0,
294  std::runtime_error,
295  "Encountered zero diag element at: " << ierr);
296  TEUCHOS_TEST_FOR_EXCEPTION( ierr == -1,
297  std::runtime_error,
298  "Could not alloc needed working memory for solve" );
299 
300  {
301 #ifdef HAVE_AMESOS2_TIMERS
302  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
303 #endif
304 
306  MultiVecAdapter<Vector>,slu_type>::do_put(X, xvals_(),
307  as<size_t>(ld_rhs),
308  ROOTED);
309  }
310 
311  return(ierr);
312 }
313 
314 
315 template <class Matrix, class Vector>
316 bool
318 {
319  // The Basker can only handle square for right now
320  return( this->globalNumRows_ == this->globalNumCols_ );
321 }
322 
323 
324 template <class Matrix, class Vector>
325 void
326 Basker<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
327 {
328  using Teuchos::RCP;
329  using Teuchos::getIntegralValue;
330  using Teuchos::ParameterEntryValidator;
331 
332  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
333  /*To Do --- add support for parameters */
334  if(parameterList->isParameter("num_threads"))
335  {
336  num_threads = parameterList->get<int>("num_threads");
337  }
338 }
339 
340 template <class Matrix, class Vector>
341 Teuchos::RCP<const Teuchos::ParameterList>
343 {
344  using Teuchos::ParameterList;
345 
346  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
347 
348 
349 #ifdef SHYLUBASKER
350  if( is_null(valid_params) )
351  {
352  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
353  pl->set("num_threads", 1, "Number of threads");
354  //This is all I am supporting right this second in experimental
355  }
356  return valid_params;
357 #else
358  if( is_null(valid_params) ){
359  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
360 
361  pl->set("alnnz", 2, "Approx number of nonzeros in L, default is 2*nnz(A)");
362  pl->set("aunnx", 2, "Approx number of nonzeros in I, default is 2*nnz(U)");
363  valid_params = pl;
364  }
365 #endif
366  return valid_params;
367 }
368 
369 
370 template <class Matrix, class Vector>
371 bool
373 {
374  using Teuchos::as;
375 
376  if(current_phase == SOLVE) return (false);
377 
378 #ifdef HAVE_AMESOS2_TIMERS
379  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
380 #endif
381 
382 
383 
384  // Only the root image needs storage allocated
385  if( this->root_ ){
386  nzvals_.resize(this->globalNumNonZeros_);
387  rowind_.resize(this->globalNumNonZeros_);
388  colptr_.resize(this->globalNumCols_ + 1);
389  }
390 
391  local_ordinal_type nnz_ret = 0;
392  {
393 #ifdef HAVE_AMESOS2_TIMERS
394  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
395 #endif
396 
398  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
399  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
400  nnz_ret, ROOTED, ARBITRARY);
401  }
402 
403 
404  if( this->root_ ){
405  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
406  std::runtime_error,
407  "Did not get the expected number of non-zero vals");
408  }
409 
410  return true;
411 }
412 
413 
414 template<class Matrix, class Vector>
415 const char* Basker<Matrix,Vector>::name = "Basker";
416 
417 
418 } // end namespace Amesos2
419 
420 #endif // AMESOS2_Basker_DEF_HPP
Definition: basker.cpp:35
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
bool root_
If true, then this is the root processor.
Definition: Amesos2_SolverCore_decl.hpp:507
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
Teuchos::Array< slu_type > bvals_
Persisting 1D store for B.
Definition: Amesos2_Basker_decl.hpp:192
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Basker_def.hpp:372
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
Definition: Amesos2_TypeDecl.hpp:142
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Basker_def.hpp:342
Teuchos::Array< local_ordinal_type > colptr_
Stores the row indices of the nonzero entries.
Definition: Amesos2_Basker_decl.hpp:187
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns a pointer to the Teuchos::Comm communicator with this operator.
Definition: Amesos2_SolverCore_decl.hpp:363
Amesos2 interface to the Baker package.
Definition: Amesos2_Basker_decl.hpp:75
Amesos2 Basker declarations.
Teuchos::Array< local_ordinal_type > rowind_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_Basker_decl.hpp:185
Teuchos::Array< slu_type > xvals_
Persisting 1D store for X.
Definition: Amesos2_Basker_decl.hpp:190
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
Teuchos::Array< slu_type > nzvals_
Stores the values of the nonzero entries for Basker.
Definition: Amesos2_Basker_decl.hpp:183
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Basker_def.hpp:317
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Basker_def.hpp:127
Definition: Amesos2_TypeDecl.hpp:127
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
int numericFactorization_impl()
Basker specific numeric factorization.
Definition: Amesos2_Basker_def.hpp:184
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
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455
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:246