Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Tacho_def.hpp
1 // @HEADER
2 // *****************************************************************************
3 // Amesos2: Templated Direct Sparse Solver Package
4 //
5 // Copyright 2011 NTESS and the Amesos2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef AMESOS2_TACHO_DEF_HPP
11 #define AMESOS2_TACHO_DEF_HPP
12 
13 #include <Teuchos_Tuple.hpp>
14 #include <Teuchos_ParameterList.hpp>
15 #include <Teuchos_StandardParameterEntryValidators.hpp>
16 
18 #include "Amesos2_Tacho_decl.hpp"
19 #include "Amesos2_Util.hpp"
20 
21 namespace Amesos2 {
22 
23 template <class Matrix, class Vector>
25  Teuchos::RCP<const Matrix> A,
26  Teuchos::RCP<Vector> X,
27  Teuchos::RCP<const Vector> B )
28  : SolverCore<Amesos2::TachoSolver,Matrix,Vector>(A, X, B)
29 {
30  data_.method = 1; // Cholesky
31  data_.variant = 2; // solver variant
32  data_.streams = 1; // # of streams
33  data_.dofs_per_node = 1; // DoFs / node
34  data_.pivot_pert = false; // Diagonal pertubation
35  data_.verbose = false; // verbose
36 }
37 
38 
39 template <class Matrix, class Vector>
41 {
42  if ( this->root_ ) {
43  data_.solver.release();
44  }
45 }
46 
47 template <class Matrix, class Vector>
48 std::string
50 {
51  std::ostringstream oss;
52  oss << "Tacho solver interface";
53  return oss.str();
54 }
55 
56 template<class Matrix, class Vector>
57 int
59 {
60  return(0);
61 }
62 
63 template <class Matrix, class Vector>
64 int
66 {
67 #ifdef HAVE_AMESOS2_TIMERS
68  Teuchos::TimeMonitor symFactTime( this->timers_.symFactTime_ );
69 #endif
70 
71  int status = 0;
72  if ( this->root_ ) {
73  if(do_optimization()) {
74  this->matrixA_->returnRowPtr_kokkos_view(host_row_ptr_view_);
75  this->matrixA_->returnColInd_kokkos_view(host_cols_view_);
76  }
77 
78  data_.solver.setSolutionMethod(data_.method);
79  data_.solver.setLevelSetOptionAlgorithmVariant(data_.variant);
80  data_.solver.setSmallProblemThresholdsize(data_.small_problem_threshold_size);
81  data_.solver.setVerbose(data_.verbose);
82  data_.solver.setLevelSetOptionNumStreams(data_.streams);
83  // TODO: Confirm param options
84  // data_.solver.setMaxNumberOfSuperblocks(data_.max_num_superblocks);
85 
86  // Symbolic factorization currently must be done on host
87  if (data_.dofs_per_node > 1) {
88  data_.solver.analyze(this->globalNumCols_, data_.dofs_per_node, host_row_ptr_view_, host_cols_view_);
89  } else {
90  data_.solver.analyze(this->globalNumCols_, host_row_ptr_view_, host_cols_view_);
91  }
92  data_.solver.initialize();
93  }
94  return status;
95 }
96 
97 
98 template <class Matrix, class Vector>
99 int
101 {
102 #ifdef HAVE_AMESOS2_TIMERS
103  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
104 #endif
105 
106  int status = 0;
107  if ( this->root_ ) {
108  if(do_optimization()) {
109  this->matrixA_->returnValues_kokkos_view(device_nzvals_view_);
110  }
111  if (data_.pivot_pert) {
112  data_.solver.useDefaultPivotTolerance();
113  } else {
114  data_.solver.useNoPivotTolerance();
115  }
116  data_.solver.factorize(device_nzvals_view_);
117  }
118  return status;
119 }
120 
121 template <class Matrix, class Vector>
122 int
124  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
125 {
126  using Teuchos::as;
127 
128  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
129  const size_t nrhs = X->getGlobalNumVectors();
130 
131  // don't allocate b since it's handled by the copy manager and might just be
132  // be assigned, not copied anyways.
133  // also don't allocate x since we will also use do_get to allocate this if
134  // necessary. When a copy is not necessary we'll solve directly to the x
135  // values in the MV.
136  bool bDidAssignX;
137  { // Get values from RHS B
138 #ifdef HAVE_AMESOS2_TIMERS
139  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
140 #endif
141  const bool initialize_data = true;
142  const bool do_not_initialize_data = false;
143  Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
144  device_solve_array_t>::do_get(initialize_data, B, this->bValues_,
145  as<size_t>(ld_rhs),
146  ROOTED, this->rowIndexBase_);
147  bDidAssignX = Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
148  device_solve_array_t>::do_get(do_not_initialize_data, X, this->xValues_,
149  as<size_t>(ld_rhs),
150  ROOTED, this->rowIndexBase_);
151  }
152 
153  int ierr = 0; // returned error code
154 
155  if ( this->root_ ) { // Do solve!
156  // Bump up the workspace size if needed
157 #ifdef HAVE_AMESOS2_TIMERS
158  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
159 #endif
160  if (workspace_.extent(0) < this->globalNumRows_ || workspace_.extent(1) < nrhs) {
161  workspace_ = device_solve_array_t(
162  Kokkos::ViewAllocateWithoutInitializing("t"), this->globalNumRows_, nrhs);
163  }
164 
165  data_.solver.solve(xValues_, bValues_, workspace_);
166 
167  int status = 0; // TODO: determine what error handling will be
168  if(status != 0) {
169  ierr = status;
170  }
171  }
172 
173  /* All processes should have the same error code */
174  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
175 
176  TEUCHOS_TEST_FOR_EXCEPTION( ierr != 0, std::runtime_error,
177  "tacho_solve has error code: " << ierr );
178 
179  /* Update X's global values */
180 
181  // if bDidAssignX, then we solved straight to the adapter's X memory space without
182  // requiring additional memory allocation, so the x data is already in place.
183  if(!bDidAssignX) {
184 #ifdef HAVE_AMESOS2_TIMERS
185  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
186 #endif
187 
188  // This will do nothing is if the target view matches the src view, which
189  // can be the case if the memory spaces match. See comments above for do_get.
190  Util::template put_1d_data_helper_kokkos_view<
191  MultiVecAdapter<Vector>,device_solve_array_t>::do_put(X, xValues_,
192  as<size_t>(ld_rhs),
193  ROOTED, this->rowIndexBase_);
194  }
195 
196  return(ierr);
197 }
198 
199 
200 template <class Matrix, class Vector>
201 bool
203 {
204  // Tacho can only apply the solve routines to square matrices
205  return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
206 }
207 
208 
209 template <class Matrix, class Vector>
210 void
211 TachoSolver<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
212 {
213  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
214 
215  // TODO: Confirm param options
216 
217  // factorization type
218  auto method_name = parameterList->get<std::string> ("method", "chol");
219  if (method_name == "chol")
220  data_.method = 1;
221  else if (method_name == "ldl")
222  data_.method = 2;
223  else if (method_name == "lu")
224  data_.method = 3;
225  else {
226  std::cout << "Error: not supported solution method\n";
227  }
228  // solver type
229  data_.variant = parameterList->get<int> ("variant", 2);
230  // small problem threshold
231  data_.small_problem_threshold_size = parameterList->get<int> ("small problem threshold size", 1024);
232  // verbosity
233  data_.verbose = parameterList->get<bool> ("verbose", false);
234  // # of streams
235  data_.streams = parameterList->get<int> ("num-streams", 1);
236  // DoFs / node
237  data_.dofs_per_node = parameterList->get<int> ("dofs-per-node", 1);
238  // Perturb tiny pivots
239  data_.pivot_pert = parameterList->get<bool> ("perturb-pivot", false);
240  // TODO: Confirm param options
241  // data_.num_kokkos_threads = parameterList->get<int>("kokkos-threads", 1);
242  // data_.max_num_superblocks = parameterList->get<int>("max-num-superblocks", 4);
243 }
244 
245 
246 template <class Matrix, class Vector>
247 Teuchos::RCP<const Teuchos::ParameterList>
249 {
250  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
251 
252  if( is_null(valid_params) ){
253  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
254 
255  pl->set("method", "chol", "Type of factorization, chol, ldl, or lu");
256  pl->set("variant", 2, "Type of solver variant, 0, 1, or 2");
257  pl->set("small problem threshold size", 1024, "Problem size threshold below with Tacho uses LAPACK.");
258  pl->set("verbose", false, "Verbosity");
259  pl->set("num-streams", 1, "Number of GPU streams");
260  pl->set("dofs-per-node", 1, "DoFs per node");
261  pl->set("perturb-pivot", false, "Perturb tiny pivots");
262 
263  // TODO: Confirm param options
264  // pl->set("kokkos-threads", 1, "Number of threads");
265  // pl->set("max-num-superblocks", 4, "Max number of superblocks");
266 
267  valid_params = pl;
268  }
269 
270  return valid_params;
271 }
272 
273 template <class Matrix, class Vector>
274 bool
276  return (this->root_ && (this->matrixA_->getComm()->getSize() == 1));
277 }
278 
279 template <class Matrix, class Vector>
280 bool
282 {
283 
284  if(current_phase == SOLVE) {
285  return(false);
286  }
287 
288  if(!do_optimization()) {
289 #ifdef HAVE_AMESOS2_TIMERS
290  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
291 #endif
292 
293  // Note views are allocated but eventually we should remove this.
294  // The internal copy manager will decide if we can assign or deep_copy
295  // and then allocate if necessary. However the GPU solvers are serial right
296  // now so I didn't complete refactoring the matrix code for the parallel
297  // case. If we added that later, we should have it hooked up to the copy
298  // manager and then these allocations can go away.
299  if( this->root_ ) {
300  device_nzvals_view_ = device_value_type_array(
301  Kokkos::ViewAllocateWithoutInitializing("nzvals"), this->globalNumNonZeros_);
302  host_cols_view_ = host_ordinal_type_array(
303  Kokkos::ViewAllocateWithoutInitializing("colind"), this->globalNumNonZeros_);
304  host_row_ptr_view_ = host_size_type_array(
305  Kokkos::ViewAllocateWithoutInitializing("rowptr"), this->globalNumRows_ + 1);
306  }
307 
308  typename host_size_type_array::value_type nnz_ret = 0;
309  {
310  #ifdef HAVE_AMESOS2_TIMERS
311  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
312  #endif
313 
314  TEUCHOS_TEST_FOR_EXCEPTION( this->rowIndexBase_ != this->columnIndexBase_,
315  std::runtime_error,
316  "Row and column maps have different indexbase ");
317 
319  device_value_type_array, host_ordinal_type_array, host_size_type_array>::do_get(
320  this->matrixA_.ptr(),
321  device_nzvals_view_,
322  host_cols_view_,
323  host_row_ptr_view_,
324  nnz_ret,
325  ROOTED, ARBITRARY,
326  this->columnIndexBase_);
327  }
328  }
329 
330  return true;
331 }
332 
333 
334 template<class Matrix, class Vector>
335 const char* TachoSolver<Matrix,Vector>::name = "Tacho";
336 
337 
338 } // end namespace Amesos2
339 
340 #endif // AMESOS2_TACHO_DEF_HPP
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:71
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:31
Amesos2 interface to the Tacho package.
Definition: Amesos2_Tacho_decl.hpp:33
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using Tacho.
Definition: Amesos2_Tacho_def.hpp:65
Utility functions for Amesos2.
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Tacho_def.hpp:202
int numericFactorization_impl()
Tacho specific numeric factorization.
Definition: Amesos2_Tacho_def.hpp:100
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Tacho_def.hpp:58
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Tacho_def.hpp:248
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:625
std::string description() const override
Returns a short description of this Solver.
Definition: Amesos2_Tacho_def.hpp:49
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
Tacho specific solve.
Definition: Amesos2_Tacho_def.hpp:123
TachoSolver(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Tacho_def.hpp:24
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Tacho_def.hpp:281
bool do_optimization() const
can we optimize size_type and ordinal_type for straight pass through
Definition: Amesos2_Tacho_def.hpp:275
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:142
~TachoSolver()
Destructor.
Definition: Amesos2_Tacho_def.hpp:40