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
55namespace Amesos2 {
56
57template <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
69template <class Matrix, class Vector>
71{
72 if ( this->root_ ) {
73 data_.solver.release();
74 }
75}
76
77template <class Matrix, class Vector>
78std::string
80{
81 std::ostringstream oss;
82 oss << "Tacho solver interface";
83 return oss.str();
84}
85
86template<class Matrix, class Vector>
87int
92
93template <class Matrix, class Vector>
94int
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
122template <class Matrix, class Vector>
123int
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
140template <class Matrix, class Vector>
141int
142TachoSolver<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
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
220template <class Matrix, class Vector>
221bool
223{
224 // Tacho can only apply the solve routines to square matrices
225 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
226}
227
228
229template <class Matrix, class Vector>
230void
231TachoSolver<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
256template <class Matrix, class Vector>
257Teuchos::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
278template <class Matrix, class Vector>
279bool
281 return (this->root_ && (this->matrixA_->getComm()->getSize() == 1));
282}
283
284template <class Matrix, class Vector>
285bool
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,
331 this->columnIndexBase_);
332 }
333 }
334
335 return true;
336}
337
338
339template<class Matrix, class Vector>
340const char* TachoSolver<Matrix,Vector>::name = "Tacho";
341
342
343} // end namespace Amesos2
344
345#endif // AMESOS2_TACHO_DEF_HPP
@ ROOTED
Definition Amesos2_TypeDecl.hpp:127
@ ARBITRARY
Definition Amesos2_TypeDecl.hpp:143
Utility functions for Amesos2.
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition Amesos2_SolverCore_decl.hpp:106
Amesos2 interface to the Tacho package.
Definition Amesos2_Tacho_decl.hpp:68
~TachoSolver()
Destructor.
Definition Amesos2_Tacho_def.hpp:70
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
EPhase
Used to indicate a phase in the direct solution.
Definition Amesos2_TypeDecl.hpp:65
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition Amesos2_Util.hpp:663