Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Lapack_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 
44 
53 #ifndef AMESOS2_LAPACK_DEF_HPP
54 #define AMESOS2_LAPACK_DEF_HPP
55 
56 #include <Teuchos_RCP.hpp>
57 
59 #include "Amesos2_Lapack_decl.hpp"
60 
61 namespace Amesos2 {
62 
63 
64  template <class Matrix, class Vector>
65  Lapack<Matrix,Vector>::Lapack(Teuchos::RCP<const Matrix> A,
66  Teuchos::RCP<Vector> X,
67  Teuchos::RCP<const Vector> B)
68  : SolverCore<Amesos2::Lapack,Matrix,Vector>(A, X, B) // instantiate superclass
69  , nzvals_()
70  , rowind_()
71  , colptr_()
72  , is_contiguous_(true)
73  {
74  // Set default parameters
75  Teuchos::RCP<Teuchos::ParameterList> default_params
76  = Teuchos::parameterList( *(this->getValidParameters()) );
77  this->setParameters(default_params);
78  }
79 
80 
81  template <class Matrix, class Vector>
83  {
84  /*
85  * Free any memory allocated by the Lapack library functions (i.e. none)
86  */
87  }
88 
89 
90  template<class Matrix, class Vector>
91  int
93  {
94  return(0);
95  }
96 
97 
98  template <class Matrix, class Vector>
99  int
101  {
102  return(0);
103  }
104 
105 
106  template <class Matrix, class Vector>
107  int
109  {
110  int factor_ierr = 0;
111 
112  if( this->root_ ){
113  // Set here so that solver_ can refresh it's internal state
114  solver_.setMatrix( Teuchos::rcpFromRef(lu_) );
115 
116  {
117 #ifdef HAVE_AMESOS2_TIMERS
118  Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
119 #endif
120  factor_ierr = solver_.factor();
121  }
122  }
123 
124  Teuchos::broadcast(*(this->getComm()), 0, &factor_ierr);
125  TEUCHOS_TEST_FOR_EXCEPTION( factor_ierr != 0,
126  std::runtime_error,
127  "Lapack factor routine returned error code "
128  << factor_ierr );
129  return( 0 );
130  }
131 
132 
133  template <class Matrix, class Vector>
134  int
136  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
137  {
138  using Teuchos::as;
139 
140  // Convert X and B to SerialDenseMatrix's
141  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
142  const size_t nrhs = X->getGlobalNumVectors();
143 
144  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
145  if( this->root_ ){
146  rhsvals_.resize(val_store_size);
147  }
148 
149  { // Get values from RHS B
150 #ifdef HAVE_AMESOS2_TIMERS
151  Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
152  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
153 #endif
155  scalar_type> copy_helper;
156  if ( is_contiguous_ == true ) {
157  copy_helper::do_get(B, rhsvals_(), as<size_t>(ld_rhs), ROOTED, 0);
158  }
159  else {
160  copy_helper::do_get(B, rhsvals_(), as<size_t>(ld_rhs), CONTIGUOUS_AND_ROOTED, 0);
161  }
162  }
163 
164  int solve_ierr = 0;
165  // int unequilibrate_ierr = 0; // unused
166 
167  if( this->root_ ){
168 #ifdef HAVE_AMESOS2_TIMERS
169  Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
170 #endif
171 
172  using Teuchos::rcpFromRef;
173  typedef Teuchos::SerialDenseMatrix<int,scalar_type> DenseMat;
174 
175  DenseMat rhs_dense_mat(Teuchos::View, rhsvals_.getRawPtr(),
176  as<int>(ld_rhs), as<int>(ld_rhs), as<int>(nrhs));
177 
178  solver_.setVectors( rcpFromRef(rhs_dense_mat),
179  rcpFromRef(rhs_dense_mat) );
180 
181  solve_ierr = solver_.solve();
182 
183  // Solution is found in rhsvals_
184  }
185 
186  // Consolidate and check error codes
187  Teuchos::broadcast(*(this->getComm()), 0, &solve_ierr);
188  TEUCHOS_TEST_FOR_EXCEPTION( solve_ierr != 0,
189  std::runtime_error,
190  "Lapack solver solve method returned with error code "
191  << solve_ierr );
192 
193  /* Update X's global values */
194  {
195 #ifdef HAVE_AMESOS2_TIMERS
196  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
197 #endif
198 
199  if ( is_contiguous_ == true ) {
201  MultiVecAdapter<Vector>,scalar_type>::do_put(X, rhsvals_(),
202  as<size_t>(ld_rhs),
203  ROOTED);
204  }
205  else {
207  MultiVecAdapter<Vector>,scalar_type>::do_put(X, rhsvals_(),
208  as<size_t>(ld_rhs),
209  CONTIGUOUS_AND_ROOTED);
210  }
211  }
212 
213  return( 0 );
214  }
215 
216 
217  template <class Matrix, class Vector>
218  bool
220  {
221  // Factorization of rectangular matrices is supported, but not
222  // their solution. For solution we can have square matrices.
223 
224  return( this->globalNumCols_ == this->globalNumRows_ );
225  }
226 
227 
228  template <class Matrix, class Vector>
229  void
230  Lapack<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
231  {
232  solver_.solveWithTranspose( parameterList->get<bool>("Transpose",
233  this->control_.useTranspose_) );
234 
235  solver_.factorWithEquilibration( parameterList->get<bool>("Equilibrate", true) );
236 
237  if( parameterList->isParameter("IsContiguous") ){
238  is_contiguous_ = parameterList->get<bool>("IsContiguous");
239  }
240 
241  // solver_.solveToRefinedSolution( parameterList->get<bool>("Refine", false) );
242  }
243 
244  template <class Matrix, class Vector>
245  Teuchos::RCP<const Teuchos::ParameterList>
247  {
248  using Teuchos::ParameterList;
249 
250  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
251 
252  if( is_null(valid_params) ){
253  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
254 
255  pl->set("Equilibrate", true, "Whether to equilibrate the input matrix");
256 
257  pl->set("IsContiguous", true, "Whether GIDs contiguous");
258 
259  // TODO: Refinement does not seem to be working with the SerialDenseSolver. Not sure why.
260  // pl->set("Refine", false, "Whether to apply iterative refinement");
261 
262  valid_params = pl;
263  }
264 
265  return valid_params;
266  }
267 
268  template <class Matrix, class Vector>
269  bool
271  {
272 #ifdef HAVE_AMESOS2_TIMERS
273  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
274 #endif
275 
276  // We only load the matrix when numericFactorization is called
277  if( current_phase < NUMFACT ) return( false );
278 
279  if( this->root_ ){
280  nzvals_.resize(this->globalNumNonZeros_);
281  rowind_.resize(this->globalNumNonZeros_);
282  colptr_.resize(this->globalNumCols_ + 1);
283  }
284 
285  // global_size_type nnz_ret = 0;
286  int nnz_ret = 0;
287  {
288 #ifdef HAVE_AMESOS2_TIMERS
289  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
290 #endif
291 
292  // typedef Util::get_ccs_helper<MatrixAdapter<Matrix>,
293  // scalar_type, global_ordinal_type, global_size_type> ccs_helper;
295  scalar_type, int, int> ccs_helper;
296  if ( is_contiguous_ == true ) {
297  ccs_helper::do_get(this->matrixA_.ptr(),
298  nzvals_(), rowind_(), colptr_(),
299  nnz_ret, ROOTED, SORTED_INDICES, 0);
300  }
301  else {
302  ccs_helper::do_get(this->matrixA_.ptr(),
303  nzvals_(), rowind_(), colptr_(),
304  nnz_ret, CONTIGUOUS_AND_ROOTED, SORTED_INDICES, 0);
305  }
306  }
307 
308  if( this->root_ ){
309  // entries are initialized to zero in here:
310  lu_.shape(this->globalNumRows_, this->globalNumCols_);
311 
312  // Put entries of ccs representation into the dense matrix
313  global_size_type end_col = this->globalNumCols_;
314  for( global_size_type col = 0; col < end_col; ++col ){
315  global_ordinal_type ptr = colptr_[col];
316  global_ordinal_type end_ptr = colptr_[col+1];
317  for( ; ptr < end_ptr; ++ptr ){
318  lu_(rowind_[ptr], col) = nzvals_[ptr];
319  }
320  }
321 
322  // lu_.print(std::cout);
323  }
324 
325  return( true );
326 }
327 
328 
329  template<class Matrix, class Vector>
330  const char* Lapack<Matrix,Vector>::name = "LAPACK";
331 
332 
333 } // end namespace Amesos2
334 
335 #endif // AMESOS2_LAPACK_DEF_HPP
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
Amesos2 interface to the LAPACK.
Definition: Amesos2_Lapack_decl.hpp:80
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
int symbolicFactorization_impl()
No-op.
Definition: Amesos2_Lapack_def.hpp:100
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:266
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_Lapack_def.hpp:230
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a const parameter list of all of the valid parameters that this-&gt;setParameterList(...) will accept.
Definition: Amesos2_SolverCore_def.hpp:314
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Lapack_def.hpp:219
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:765
Lapack(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Lapack_def.hpp:65
int preOrdering_impl()
No-op.
Definition: Amesos2_Lapack_def.hpp:92
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Lapack_def.hpp:246
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
Lapack solve.
Definition: Amesos2_Lapack_def.hpp:135
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:282
Declarations for the Amesos2 interface to LAPACK.
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Lapack_def.hpp:270
~Lapack()
Destructor.
Definition: Amesos2_Lapack_def.hpp:82
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:344
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
int numericFactorization_impl()
Perform numeric factorization using LAPACK.
Definition: Amesos2_Lapack_def.hpp:108