Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_KLU2_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 
52 #ifndef AMESOS2_KLU2_DEF_HPP
53 #define AMESOS2_KLU2_DEF_HPP
54 
55 #include <Teuchos_Tuple.hpp>
56 #include <Teuchos_ParameterList.hpp>
57 #include <Teuchos_StandardParameterEntryValidators.hpp>
58 
60 #include "Amesos2_KLU2_decl.hpp"
61 
62 namespace Amesos2 {
63 
64 
65 template <class Matrix, class Vector>
67  Teuchos::RCP<const Matrix> A,
68  Teuchos::RCP<Vector> X,
69  Teuchos::RCP<const Vector> B )
70  : SolverCore<Amesos2::KLU2,Matrix,Vector>(A, X, B)
71  , nzvals_() // initialize to empty arrays
72  , rowind_()
73  , colptr_()
74 {
75  ::KLU2::klu_defaults<slu_type, local_ordinal_type> (&(data_.common_)) ;
76  data_.symbolic_ = NULL;
77  data_.numeric_ = NULL;
78 
79  // Override some default options
80  // TODO: use data_ here to init
81 }
82 
83 
84 template <class Matrix, class Vector>
86 {
87  /* Free KLU2 data_types
88  * - Matrices
89  * - Vectors
90  * - Other data
91  */
92  if (data_.symbolic_ != NULL)
93  ::KLU2::klu_free_symbolic<slu_type, local_ordinal_type>
94  (&(data_.symbolic_), &(data_.common_)) ;
95  if (data_.numeric_ != NULL)
96  ::KLU2::klu_free_numeric<slu_type, local_ordinal_type>
97  (&(data_.numeric_), &(data_.common_)) ;
98 
99  // Storage is initialized in numericFactorization_impl()
100  //if ( data_.A.Store != NULL ){
101  // destoy
102  //}
103 
104  // only root allocated these SuperMatrices.
105  //if ( data_.L.Store != NULL ){ // will only be true for this->root_
106  // destroy ..
107  //}
108 }
109 
110 template<class Matrix, class Vector>
111 int
113 {
114  /* TODO: Define what it means for KLU2
115  */
116 #ifdef HAVE_AMESOS2_TIMERS
117  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
118 #endif
119 
120  return(0);
121 }
122 
123 
124 template <class Matrix, class Vector>
125 int
127 {
128  data_.symbolic_ = ::KLU2::klu_analyze<slu_type, local_ordinal_type>
129  ((local_ordinal_type)this->globalNumCols_, colptr_.getRawPtr(),
130  rowind_.getRawPtr(), &(data_.common_)) ;
131 
132  return(0);
133 }
134 
135 
136 template <class Matrix, class Vector>
137 int
139 {
140  using Teuchos::as;
141 
142  // Cleanup old L and U matrices if we are not reusing a symbolic
143  // factorization. Stores and other data will be allocated in gstrf.
144  // Only rank 0 has valid pointers, TODO: for KLU2
145 
146 
147  int info = 0;
148  if ( this->root_ ){
149 
150  { // Do factorization
151 #ifdef HAVE_AMESOS2_TIMERS
152  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
153 #endif
154 
155 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
156  std::cout << "KLU2:: Before numeric factorization" << std::endl;
157  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
158  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
159  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
160 #endif
161 
162  data_.numeric_ = ::KLU2::klu_factor<slu_type, local_ordinal_type>
163  (colptr_.getRawPtr(), rowind_.getRawPtr(), nzvals_.getRawPtr(),
164  data_.symbolic_, &(data_.common_)) ;
165 
166  }
167 
168  }
169 
170  /* All processes should have the same error code */
171  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
172 
173  //global_size_type info_st = as<global_size_type>(info); // unused
174  /* TODO : Proper error messages
175  TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
176  std::runtime_error,
177  "Factorization complete, but matrix is singular. Division by zero eminent");
178  TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
179  std::runtime_error,
180  "Memory allocation failure in KLU2 factorization");*/
181 
182  //data_.options.Fact = SLU::FACTORED;
183  //same_symbolic_ = true;
184 
185  return(info);
186 }
187 
188 
189 template <class Matrix, class Vector>
190 int
192  const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
193  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
194 {
195  using Teuchos::as;
196 
197  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
198  const size_t nrhs = X->getGlobalNumVectors();
199 
200  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
201  Teuchos::Array<slu_type> bValues(val_store_size);
202 
203  { // Get values from RHS B
204 #ifdef HAVE_AMESOS2_TIMERS
205  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
206  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
207 #endif
209  slu_type>::do_get(B, bValues(),
210  as<size_t>(ld_rhs),
211  ROOTED);
212  }
213 
214 
215  int ierr = 0; // returned error code
216 
217  if ( this->root_ ) {
218 
219  //local_ordinal_type i_ld_rhs = as<local_ordinal_type>(ld_rhs);
220 
221  { // Do solve!
222 #ifdef HAVE_AMESOS2_TIMERS
223  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
224 #endif
225  ::KLU2::klu_solve<slu_type, local_ordinal_type>
226  (data_.symbolic_, data_.numeric_,
227  (local_ordinal_type)this->globalNumCols_,
228  (local_ordinal_type)nrhs,
229  bValues.getRawPtr(), &(data_.common_)) ;
230 
231  }
232 
233  }
234 
235  /* All processes should have the same error code */
236  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
237 
238  // global_size_type ierr_st = as<global_size_type>(ierr); // unused
239  // TODO
240  //TEUCHOS_TEST_FOR_EXCEPTION( ierr < 0,
241  //std::invalid_argument,
242  //"Argument " << -ierr << " to KLU2 xgssvx had illegal value" );
243  //TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st <= this->globalNumCols_,
244  //std::runtime_error,
245  //"Factorization complete, but U is exactly singular" );
246  //TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st > this->globalNumCols_ + 1,
247  //std::runtime_error,
248  //"KLU2 allocated " << ierr - this->globalNumCols_ << " bytes of "
249  //"memory before allocation failure occured." );
250 
251  /* Update X's global values */
252  {
253 #ifdef HAVE_AMESOS2_TIMERS
254  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
255 #endif
256 
258  MultiVecAdapter<Vector>,slu_type>::do_put(X, bValues(),
259  as<size_t>(ld_rhs),
260  ROOTED);
261  }
262 
263 
264  return(ierr);
265 }
266 
267 
268 template <class Matrix, class Vector>
269 bool
271 {
272  // The KLU2 factorization routines can handle square as well as
273  // rectangular matrices, but KLU2 can only apply the solve routines to
274  // square matrices, so we check the matrix for squareness.
275  return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
276 }
277 
278 
279 template <class Matrix, class Vector>
280 void
281 KLU2<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
282 {
283  using Teuchos::RCP;
284  using Teuchos::getIntegralValue;
285  using Teuchos::ParameterEntryValidator;
286 
287  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
288 
289  if(parameterList->isParameter("Trans"))
290  {}
291 
292 
293 
294 
295  // The KLU2 transpose option can override the Amesos2 option
296  //if( parameterList->isParameter("Trans") ){
297  //RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("Trans").validator();
298  //parameterList->getEntry("Trans").setValidator(trans_validator);
299 
300  //data_.options.Trans = getIntegralValue<SLU::trans_t>(*parameterList, "Trans");
301  //}
302 
303 }
304 
305 
306 template <class Matrix, class Vector>
307 Teuchos::RCP<const Teuchos::ParameterList>
309 {
310  using Teuchos::ParameterList;
311 
312  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
313 
314  if( is_null(valid_params) ){
315  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
316 
317  pl->set("Equil", true, "Whether to equilibrate the system before solve, does nothing now");
318 
319  valid_params = pl;
320  }
321 
322  return valid_params;
323 }
324 
325 
326 template <class Matrix, class Vector>
327 bool
329 {
330  using Teuchos::as;
331 
332  if(current_phase == SOLVE)return(false);
333 
334 #ifdef HAVE_AMESOS2_TIMERS
335  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
336 #endif
337 
338  // Only the root image needs storage allocated
339  if( this->root_ ){
340  nzvals_.resize(this->globalNumNonZeros_);
341  rowind_.resize(this->globalNumNonZeros_);
342  colptr_.resize(this->globalNumCols_ + 1);
343  }
344 
345  local_ordinal_type nnz_ret = 0;
346  {
347 #ifdef HAVE_AMESOS2_TIMERS
348  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
349 #endif
350 
352  MatrixAdapter<Matrix>,slu_type,local_ordinal_type,local_ordinal_type>
353  ::do_get(this->matrixA_.ptr(), nzvals_(), rowind_(), colptr_(),
354  nnz_ret, ROOTED, ARBITRARY);
355  }
356 
357 
358  if( this->root_ ){
359  TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<local_ordinal_type>(this->globalNumNonZeros_),
360  std::runtime_error,
361  "Did not get the expected number of non-zero vals");
362  }
363 
364  return true;
365 }
366 
367 
368 template<class Matrix, class Vector>
369 const char* KLU2<Matrix,Vector>::name = "KLU2";
370 
371 
372 } // end namespace Amesos2
373 
374 #endif // AMESOS2_KLU2_DEF_HPP
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
Teuchos::Array< slu_type > nzvals_
Stores the values of the nonzero entries for KLU2.
Definition: Amesos2_KLU2_decl.hpp:230
KLU2(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_KLU2_def.hpp:66
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
KLU2 specific solve.
Definition: Amesos2_KLU2_def.hpp:191
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
Amesos2 KLU2 declarations.
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_KLU2_def.hpp:328
bool root_
If true, then this is the root processor.
Definition: Amesos2_SolverCore_decl.hpp:507
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
~KLU2()
Destructor.
Definition: Amesos2_KLU2_def.hpp:85
Definition: Amesos2_TypeDecl.hpp:142
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_KLU2_def.hpp:281
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using KLU2.
Definition: Amesos2_KLU2_def.hpp:126
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
Definition: Amesos2_KLU2_FunctionMap.hpp:67
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_KLU2_def.hpp:112
Teuchos::Array< local_ordinal_type > rowind_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_KLU2_decl.hpp:232
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns a pointer to the Teuchos::Comm communicator with this operator.
Definition: Amesos2_SolverCore_decl.hpp:363
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_KLU2_def.hpp:270
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
int numericFactorization_impl()
KLU2 specific numeric factorization.
Definition: Amesos2_KLU2_def.hpp:138
Teuchos::Array< local_ordinal_type > colptr_
Stores the row indices of the nonzero entries.
Definition: Amesos2_KLU2_decl.hpp:234
Amesos2 interface to the KLU2 package.
Definition: Amesos2_KLU2_decl.hpp:72
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_KLU2_def.hpp:308
Definition: Amesos2_TypeDecl.hpp:127
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:296
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:175
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455