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