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  {
73  // Set default parameters
74  Teuchos::RCP<Teuchos::ParameterList> default_params
75  = Teuchos::parameterList( *(this->getValidParameters()) );
76  this->setParameters(default_params);
77  }
78 
79 
80  template <class Matrix, class Vector>
82  {
83  /*
84  * Free any memory allocated by the Lapack library functions (i.e. none)
85  */
86  }
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  return(0);
102  }
103 
104 
105  template <class Matrix, class Vector>
106  int
108  {
109  int factor_ierr = 0;
110 
111  if( this->root_ ){
112  // Set here so that solver_ can refresh it's internal state
113  solver_.setMatrix( Teuchos::rcpFromRef(lu_) );
114 
115  {
116 #ifdef HAVE_AMESOS2_TIMERS
117  Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
118 #endif
119  factor_ierr = solver_.factor();
120  }
121  }
122 
123  Teuchos::broadcast(*(this->getComm()), 0, &factor_ierr);
124  TEUCHOS_TEST_FOR_EXCEPTION( factor_ierr != 0,
125  std::runtime_error,
126  "Lapack factor routine returned error code "
127  << factor_ierr );
128  return( 0 );
129  }
130 
131 
132  template <class Matrix, class Vector>
133  int
135  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
136  {
137  using Teuchos::as;
138 
139  // Convert X and B to SerialDenseMatrix's
140  const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
141  const size_t nrhs = X->getGlobalNumVectors();
142 
143  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
144  if( this->root_ ){
145  rhsvals_.resize(val_store_size);
146  }
147 
148  { // Get values from RHS B
149 #ifdef HAVE_AMESOS2_TIMERS
150  Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
151  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
152 #endif
154  scalar_type> copy_helper;
155  copy_helper::do_get(B, rhsvals_(), as<size_t>(ld_rhs), ROOTED);
156  }
157 
158  int solve_ierr = 0;
159  // int unequilibrate_ierr = 0; // unused
160 
161  if( this->root_ ){
162 #ifdef HAVE_AMESOS2_TIMERS
163  Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
164 #endif
165 
166  using Teuchos::rcpFromRef;
167  typedef Teuchos::SerialDenseMatrix<int,scalar_type> DenseMat;
168 
169  DenseMat rhs_dense_mat(Teuchos::View, rhsvals_.getRawPtr(),
170  as<int>(ld_rhs), as<int>(ld_rhs), as<int>(nrhs));
171 
172  solver_.setVectors( rcpFromRef(rhs_dense_mat),
173  rcpFromRef(rhs_dense_mat) );
174 
175  solve_ierr = solver_.solve();
176 
177  // Solution is found in rhsvals_
178  }
179 
180  // Consolidate and check error codes
181  Teuchos::broadcast(*(this->getComm()), 0, &solve_ierr);
182  TEUCHOS_TEST_FOR_EXCEPTION( solve_ierr != 0,
183  std::runtime_error,
184  "Lapack solver solve method returned with error code "
185  << solve_ierr );
186 
187  /* Update X's global values */
188  {
189 #ifdef HAVE_AMESOS2_TIMERS
190  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
191 #endif
192 
194  MultiVecAdapter<Vector>,scalar_type>::do_put(X, rhsvals_(),
195  as<size_t>(ld_rhs),
196  ROOTED);
197  }
198 
199  return( 0 );
200  }
201 
202 
203  template <class Matrix, class Vector>
204  bool
206  {
207  // Factorization of rectangular matrices is supported, but not
208  // their solution. For solution we can have square matrices.
209 
210  return( this->globalNumCols_ == this->globalNumRows_ );
211  }
212 
213 
214  template <class Matrix, class Vector>
215  void
216  Lapack<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
217  {
218  solver_.solveWithTranspose( parameterList->get<bool>("Transpose",
219  this->control_.useTranspose_) );
220 
221  solver_.factorWithEquilibration( parameterList->get<bool>("Equilibrate", true) );
222 
223  // solver_.solveToRefinedSolution( parameterList->get<bool>("Refine", false) );
224  }
225 
226  template <class Matrix, class Vector>
227  Teuchos::RCP<const Teuchos::ParameterList>
229  {
230  using Teuchos::ParameterList;
231 
232  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
233 
234  if( is_null(valid_params) ){
235  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
236 
237  pl->set("Equilibrate", true, "Whether to equilibrate the input matrix");
238 
239  // TODO: Refinement does not seem to be working with the SerialDenseSolver. Not sure why.
240  // pl->set("Refine", false, "Whether to apply iterative refinement");
241 
242  valid_params = pl;
243  }
244 
245  return valid_params;
246  }
247 
248  template <class Matrix, class Vector>
249  bool
251  {
252 #ifdef HAVE_AMESOS2_TIMERS
253  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
254 #endif
255 
256  // We only load the matrix when numericFactorization is called
257  if( current_phase < NUMFACT ) return( false );
258 
259  if( this->root_ ){
260  nzvals_.resize(this->globalNumNonZeros_);
261  rowind_.resize(this->globalNumNonZeros_);
262  colptr_.resize(this->globalNumCols_ + 1);
263  }
264 
265  // global_size_type nnz_ret = 0;
266  int nnz_ret = 0;
267  {
268 #ifdef HAVE_AMESOS2_TIMERS
269  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
270 #endif
271 
272  // typedef Util::get_ccs_helper<MatrixAdapter<Matrix>,
273  // scalar_type, global_ordinal_type, global_size_type> ccs_helper;
275  scalar_type, int, int> ccs_helper;
276  ccs_helper::do_get(this->matrixA_.ptr(),
277  nzvals_(), rowind_(), colptr_(),
278  nnz_ret, ROOTED, SORTED_INDICES);
279  }
280 
281  if( this->root_ ){
282  // entries are initialized to zero in here:
283  lu_.shape(this->globalNumRows_, this->globalNumCols_);
284 
285  // Put entries of ccs representation into the dense matrix
286  global_size_type end_col = this->globalNumCols_;
287  for( global_size_type col = 0; col < end_col; ++col ){
288  global_ordinal_type ptr = colptr_[col];
289  global_ordinal_type end_ptr = colptr_[col+1];
290  for( ; ptr < end_ptr; ++ptr ){
291  lu_(rowind_[ptr], col) = nzvals_[ptr];
292  }
293  }
294 
295  // lu_.print(std::cout);
296  }
297 
298  return( true );
299 }
300 
301 
302  template<class Matrix, class Vector>
303  const char* Lapack<Matrix,Vector>::name = "LAPACK";
304 
305 
306 } // end namespace Amesos2
307 
308 #endif // AMESOS2_LAPACK_DEF_HPP
Teuchos::Array< int > colptr_
Stores the location in Ai_ and Aval_ that starts col j.
Definition: Amesos2_Lapack_decl.hpp:208
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
Definition: Amesos2_TypeDecl.hpp:141
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:99
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
Teuchos::SerialDenseMatrix< int, scalar_type > lu_
L and U storage.
Definition: Amesos2_Lapack_decl.hpp:214
bool root_
If true, then this is the root processor.
Definition: Amesos2_SolverCore_decl.hpp:507
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
Teuchos::SerialDenseSolver< int, scalar_type > solver_
The serial solver.
Definition: Amesos2_Lapack_decl.hpp:218
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_Lapack_def.hpp:216
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a const parameter list of all of the valid parameters that this->setParameterList(...) will accept.
Definition: Amesos2_SolverCore_def.hpp:307
Teuchos::Array< int > rowind_
Stores the column indices of the nonzero entries.
Definition: Amesos2_Lapack_decl.hpp:205
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Lapack_def.hpp:205
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
Teuchos::Array< scalar_type > nzvals_
Stores the values of the nonzero entries.
Definition: Amesos2_Lapack_decl.hpp:202
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
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:91
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Lapack_def.hpp:228
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
Lapack solve.
Definition: Amesos2_Lapack_def.hpp:134
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:275
Declarations for the Amesos2 interface to LAPACK.
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Lapack_def.hpp:250
Definition: Amesos2_TypeDecl.hpp:127
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
~Lapack()
Destructor.
Definition: Amesos2_Lapack_def.hpp:81
Teuchos::Array< scalar_type > rhsvals_
Store for RHS and Solution values.
Definition: Amesos2_Lapack_decl.hpp:210
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
int numericFactorization_impl()
Perform numeric factorization using LAPACK.
Definition: Amesos2_Lapack_def.hpp:107