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
62
63namespace Amesos2 {
64
65
66template <class Matrix, class Vector>
67Basker<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
79template <class Matrix, class Vector>
80Basker<Matrix,Vector>::~Basker( )
81{}
82
83template <class Matrix, class Vector>
84bool
86 return (this->root_ && (this->matrixA_->getComm()->getSize() == 1) && is_contiguous_);
87}
88
89template<class Matrix, class Vector>
90int
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
103template <class Matrix, class Vector>
104int
106{
107 /*No symbolic factoriztion*/
108 return(0);
109}
110
111
112template <class Matrix, class Vector>
113int
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
164template <class Matrix, class Vector>
165int
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
260template <class Matrix, class Vector>
261bool
263{
264 // The Basker can only handle square for right now
265 return( this->globalNumRows_ == this->globalNumCols_ );
266}
267
268
269template <class Matrix, class Vector>
270void
271Basker<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
286template <class Matrix, class Vector>
287Teuchos::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
306template <class Matrix, class Vector>
307bool
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
356template<class Matrix, class Vector>
357const char* Basker<Matrix,Vector>::name = "Basker";
358
359
360} // end namespace Amesos2
361
362#endif // AMESOS2_Basker_DEF_HPP
Amesos2 Basker declarations.
@ ROOTED
Definition: Amesos2_TypeDecl.hpp:127
@ CONTIGUOUS_AND_ROOTED
Definition: Amesos2_TypeDecl.hpp:128
@ ARBITRARY
Definition: Amesos2_TypeDecl.hpp:143
Amesos2 interface to the Baker package.
Definition: Amesos2_Basker_decl.hpp:74
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:652