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