Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Umfpack_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_UMFPACK_DEF_HPP
45#define AMESOS2_UMFPACK_DEF_HPP
46
47#include <Teuchos_Tuple.hpp>
48#include <Teuchos_ParameterList.hpp>
49#include <Teuchos_StandardParameterEntryValidators.hpp>
50
52#include "Amesos2_Umfpack_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::Umfpack,Matrix,Vector>(A, X, B)
63 , is_contiguous_(true)
64{
65 data_.Symbolic = NULL;
66 data_.Numeric = NULL;
67}
68
69
70template <class Matrix, class Vector>
72{
73 if (data_.Symbolic) function_map::umfpack_free_symbolic (&data_.Symbolic);
74 if (data_.Numeric) function_map::umfpack_free_numeric (&data_.Numeric);
75}
76
77template <class Matrix, class Vector>
78std::string
80{
81 std::ostringstream oss;
82 oss << "Umfpack solver interface";
83 return oss.str();
84}
85
86template<class Matrix, class Vector>
87int
92
93
94template <class Matrix, class Vector>
95int
97{
98 int status = 0;
99 if ( this->root_ ) {
100 if (data_.Symbolic) {
101 function_map::umfpack_free_symbolic(&(data_.Symbolic));
102 }
103
104 function_map::umfpack_defaults(data_.Control);
105
106 status = function_map::umfpack_symbolic(
107 this->globalNumRows_,this->globalNumCols_,
108 &(this->colptr_view_[0]), &(this->rowind_view_[0]),
109 &(this->nzvals_view_[0]), &(data_.Symbolic), data_.Control, data_.Info);
110 }
111
112 return status;
113}
114
115
116template <class Matrix, class Vector>
117int
119{
120 int status = 0;
121 if ( this->root_ ) {
122 if(!data_.Symbolic) {
123 symbolicFactorization_impl();
124 }
125
126 function_map::umfpack_defaults(data_.Control);
127
128 if (data_.Numeric) {
129 function_map::umfpack_free_numeric(&(data_.Numeric));
130 }
131
132 status = function_map::umfpack_numeric(
133 &(this->colptr_view_[0]),
134 &(this->rowind_view_[0]), &(this->nzvals_view_[0]), data_.Symbolic,
135 &(data_.Numeric), data_.Control, data_.Info);
136 }
137 return status;
138}
139
140template <class Matrix, class Vector>
141int
142Umfpack<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 const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
151 Teuchos::Array<umfpack_type> xValues(val_store_size);
152 Teuchos::Array<umfpack_type> bValues(val_store_size);
153
154 { // Get values from RHS B
155#ifdef HAVE_AMESOS2_TIMERS
156 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
157 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
158#endif
159 if ( is_contiguous_ == true ) {
161 umfpack_type>::do_get(B, bValues(),
162 as<size_t>(ld_rhs),
163 ROOTED,
164 this->rowIndexBase_);
165 }
166 else {
168 umfpack_type>::do_get(B, bValues(),
169 as<size_t>(ld_rhs),
171 this->rowIndexBase_);
172 }
173 }
174
175 int UmfpackRequest = this->control_.useTranspose_ ? UMFPACK_At : UMFPACK_A;
176
177 int ierr = 0; // returned error code
178
179 if ( this->root_ ) {
180 { // Do solve!
181#ifdef HAVE_AMESOS2_TIMER
182 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
183#endif
184 if (data_.Symbolic) {
185 function_map::umfpack_free_symbolic(&(data_.Symbolic));
186 }
187
188 // validate
189 int i_ld_rhs = as<int>(ld_rhs);
190
191 for(size_t j = 0 ; j < nrhs; j++) {
192 int status = function_map::umfpack_solve(
193 UmfpackRequest,
194 &(this->colptr_view_[0]), &(this->rowind_view_[0]), &(this->nzvals_view_[0]),
195 &xValues.getRawPtr()[j*i_ld_rhs],
196 &bValues.getRawPtr()[j*i_ld_rhs],
197 data_.Numeric, data_.Control, data_.Info);
198
199 if(status != 0) {
200 ierr = status;
201 break; // need to verify best ways to handle error in this loop
202 }
203 }
204 }
205 }
206
207 /* All processes should have the same error code */
208 Teuchos::broadcast(*(this->getComm()), 0, &ierr);
209
210 TEUCHOS_TEST_FOR_EXCEPTION( ierr != 0, std::runtime_error,
211 "umfpack_solve has error code: " << ierr );
212
213 /* Update X's global values */
214 {
215#ifdef HAVE_AMESOS2_TIMERS
216 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
217#endif
218
219 if ( is_contiguous_ == true ) {
221 MultiVecAdapter<Vector>,umfpack_type>::do_put(X, xValues(),
222 as<size_t>(ld_rhs),
223 ROOTED,
224 this->rowIndexBase_);
225 }
226 else {
228 MultiVecAdapter<Vector>,umfpack_type>::do_put(X, xValues(),
229 as<size_t>(ld_rhs),
231 this->rowIndexBase_);
232 }
233 }
234
235 return(ierr);
236}
237
238
239template <class Matrix, class Vector>
240bool
242{
243 // The Umfpack factorization routines can handle square as well as
244 // rectangular matrices, but Umfpack can only apply the solve routines to
245 // square matrices, so we check the matrix for squareness.
246 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
247}
248
249
250template <class Matrix, class Vector>
251void
252Umfpack<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
253{
254 using Teuchos::RCP;
255 using Teuchos::getIntegralValue;
256 using Teuchos::ParameterEntryValidator;
257
258 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
259
260 if( parameterList->isParameter("IsContiguous") ){
261 is_contiguous_ = parameterList->get<bool>("IsContiguous");
262 }
263}
264
265
266template <class Matrix, class Vector>
267Teuchos::RCP<const Teuchos::ParameterList>
269{
270 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
271
272 if( is_null(valid_params) ){
273 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
274
275 pl->set("IsContiguous", true, "Whether GIDs contiguous");
276
277 valid_params = pl;
278 }
279
280 return valid_params;
281}
282
283
284template <class Matrix, class Vector>
285bool
287{
288 if(current_phase == SOLVE) {
289 return(false);
290 }
291
292#ifdef HAVE_AMESOS2_TIMERS
293 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
294#endif
295
296 // Only the root image needs storage allocated
297 if( this->root_ ){
298 Kokkos::resize(nzvals_view_,this->globalNumNonZeros_);
299 Kokkos::resize(rowind_view_,this->globalNumNonZeros_);
300 Kokkos::resize(colptr_view_,this->globalNumCols_ + 1);
301 }
302
303 int nnz_ret = 0;
304 {
305#ifdef HAVE_AMESOS2_TIMERS
306 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
307#endif
308
309 TEUCHOS_TEST_FOR_EXCEPTION( this->rowIndexBase_ != this->columnIndexBase_,
310 std::runtime_error,
311 "Row and column maps have different indexbase ");
312 if ( is_contiguous_ == true ) {
314 host_value_type_array,host_ordinal_type_array, host_size_type_array>::do_get(this->matrixA_.ptr(),
315 nzvals_view_, rowind_view_,
316 colptr_view_, nnz_ret,
317 ROOTED,
318 ARBITRARY,
319 this->rowIndexBase_);
320 } else {
322 host_value_type_array,host_ordinal_type_array, host_size_type_array>::do_get(this->matrixA_.ptr(),
323 nzvals_view_, rowind_view_,
324 colptr_view_, nnz_ret,
326 ARBITRARY,
327 this->rowIndexBase_);
328 }
329 }
330
331 return true;
332}
333
334
335template<class Matrix, class Vector>
336const char* Umfpack<Matrix,Vector>::name = "Umfpack";
337
338
339} // end namespace Amesos2
340
341#endif // AMESOS2_UMFPACK_DEF_HPP
@ ROOTED
Definition Amesos2_TypeDecl.hpp:127
@ CONTIGUOUS_AND_ROOTED
Definition Amesos2_TypeDecl.hpp:128
@ 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 Umfpack package.
Definition Amesos2_Umfpack_decl.hpp:64
~Umfpack()
Destructor.
Definition Amesos2_Umfpack_def.hpp:71
Umfpack(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition Amesos2_Umfpack_def.hpp:58
EPhase
Used to indicate a phase in the direct solution.
Definition Amesos2_TypeDecl.hpp:65
Helper class for getting 1-D copies of multivectors.
Definition Amesos2_MultiVecAdapter_decl.hpp:267
A generic helper class for getting a CCS representation of a Matrix.
Definition Amesos2_Util.hpp:652
Helper class for putting 1-D data arrays into multivectors.
Definition Amesos2_MultiVecAdapter_decl.hpp:373