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