Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Cholmod_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 
53 #ifndef AMESOS2_CHOLMOD_DEF_HPP
54 #define AMESOS2_CHOLMOD_DEF_HPP
55 
56 #include <Teuchos_Tuple.hpp>
57 #include <Teuchos_ParameterList.hpp>
58 #include <Teuchos_StandardParameterEntryValidators.hpp>
59 
61 #include "Amesos2_Cholmod_decl.hpp"
62 
63 
64 namespace Amesos2 {
65 
66 
67 template <class Matrix, class Vector>
69  Teuchos::RCP<const Matrix> A,
70  Teuchos::RCP<Vector> X,
71  Teuchos::RCP<const Vector> B )
72  : SolverCore<Amesos2::Cholmod,Matrix,Vector>(A, X, B)
73  , nzvals_() // initialize to empty arrays
74  , rowind_()
75  , colptr_()
76  , firstsolve(true)
77  , map()
78 {
79  CHOL::cholmod_start(&data_.c);
80  (&data_.c)->nmethods=9;
81  (&data_.c)->quick_return_if_not_posdef=1;
82 
83  //(&data_.c)->method[0].ordering=CHOLMOD_NATURAL;
84  //(&data_.c)->print=5;
85  data_.Y = NULL;
86  data_.E = NULL;
87 }
88 
89 
90 template <class Matrix, class Vector>
92 {
93  CHOL::cholmod_free_factor (&(data_.L), &(data_.c));
94 
95  CHOL::cholmod_free_dense(&(data_.Y), &data_.c);
96  CHOL::cholmod_free_dense(&(data_.E), &data_.c);
97 
98  CHOL::cholmod_finish(&(data_.c));
99 }
100 
101 template<class Matrix, class Vector>
102 int
104 {
105  #ifdef HAVE_AMESOS2_TIMERS
106  Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
107 #endif
108  data_.L = CHOL::cholmod_analyze (&data_.A, &(data_.c));
109  //data_.L->is_super=1;
110  data_.L->is_ll=1;
111 
112  skip_symfact = true;
113 
114  return(0);
115 }
116 
117 
118 template <class Matrix, class Vector>
119 int
121 {
122  if ( !skip_symfact ){
123 #ifdef HAVE_AMESOS2_TIMERS
124  Teuchos::TimeMonitor symFactTimer(this->timers_.symFactTime_);
125 #endif
126  CHOL::cholmod_resymbol (&data_.A, NULL, 0, true, data_.L, &(data_.c));
127  } else {
128  /*
129  * Symbolic factorization has already occured in preOrdering_impl,
130  * but if the user calls this routine directly later, we need to
131  * redo the symbolic factorization.
132  */
133  skip_symfact = false;
134  }
135 
136  return(0);
137 }
138 
139 
140 template <class Matrix, class Vector>
141 int
143 {
144  using Teuchos::as;
145 
146  int info = 0;
147 
148 #ifdef HAVE_AMESOS2_DEBUG
149  TEUCHOS_TEST_FOR_EXCEPTION(data_.A.ncol != as<size_t>(this->globalNumCols_),
150  std::runtime_error,
151  "Error in converting to cholmod_sparse: wrong number of global columns." );
152  TEUCHOS_TEST_FOR_EXCEPTION(data_.A.nrow != as<size_t>(this->globalNumRows_),
153  std::runtime_error,
154  "Error in converting to cholmod_sparse: wrong number of global rows." );
155 #endif
156 
157  { // Do factorization
158 #ifdef HAVE_AMESOS2_TIMERS
159  Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
160 #endif
161 
162 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
163  std::cout << "Cholmod:: Before numeric factorization" << std::endl;
164  std::cout << "nzvals_ : " << nzvals_.toString() << std::endl;
165  std::cout << "rowind_ : " << rowind_.toString() << std::endl;
166  std::cout << "colptr_ : " << colptr_.toString() << std::endl;
167 #endif
168  //cholmod_print_sparse(data_.A,"A",&(data_.c));
169  CHOL::cholmod_factorize(&data_.A, data_.L, &(data_.c));
170  //cholmod_print_factor(data_.L,"L",&(data_.c));
171  info = (&data_.c)->status;
172 
173  }
174 
175 
176 
177  /* All processes should have the same error code */
178  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
179 
180  TEUCHOS_TEST_FOR_EXCEPTION(info == 2,
181  std::runtime_error,
182  "Memory allocation failure in Cholmod factorization");
183 
184  return(info);
185 }
186 
187 
188 template <class Matrix, class Vector>
189 int
191  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
192 {
193  using Teuchos::as;
194 
195  const global_size_type ld_rhs = X->getGlobalLength();
196  const size_t nrhs = X->getGlobalNumVectors();
197 
198  Teuchos::RCP<const Tpetra::Map<local_ordinal_type,
199  global_ordinal_type, node_type> > map2;
200 
201  const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
202  Teuchos::Array<chol_type> xValues(val_store_size);
203  Teuchos::Array<chol_type> bValues(val_store_size);
204 
205  { // Get values from RHS B
206 #ifdef HAVE_AMESOS2_TIMERS
207  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
208  Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
209 #endif
211  chol_type>::do_get(B, bValues(),
212  as<size_t>(ld_rhs),
213  ROOTED, this->rowIndexBase_);
214 
215  }
216 
217 
218  int ierr = 0; // returned error code
219 
220  {
221 #ifdef HAVE_AMESOS2_TIMERS
222  Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
223 #endif
224 
225  function_map::cholmod_init_dense(as<int>(this->globalNumRows_),
226  as<int>(nrhs), as<int>(ld_rhs),
227  bValues.getRawPtr(), &data_.b);
228  function_map::cholmod_init_dense(as<int>(this->globalNumRows_),
229  as<int>(nrhs), as<int>(ld_rhs),
230  xValues.getRawPtr(), &data_.x);
231  }
232 
233 
234  { // Do solve!
235 #ifdef HAVE_AMESOS2_TIMERS
236  Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
237 #endif
238 
239  CHOL::cholmod_dense *xtemp = &(data_.x);
240 
241  CHOL::cholmod_solve2(CHOLMOD_A, data_.L, &data_.b, NULL,
242  &xtemp, NULL, &data_.Y, &data_.E, &data_.c);
243 
244  ierr = (&data_.c)->status;
245  }
246 
247 
248 
249  /* All processes should have the same error code */
250  Teuchos::broadcast(*(this->getComm()), 0, &ierr);
251 
252  TEUCHOS_TEST_FOR_EXCEPTION(ierr == -2, std::runtime_error, "Ran out of memory" );
253 
254  /* Update X's global values */
255  {
256 #ifdef HAVE_AMESOS2_TIMERS
257  Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
258 #endif
259 
261  MultiVecAdapter<Vector>,chol_type>::do_put(X, xValues(),
262  as<size_t>(ld_rhs),
263  ROOTED, this->rowIndexBase_);
264 
265  }
266 
267 
268  return(ierr);
269 }
270 
271 
272 template <class Matrix, class Vector>
273 bool
275 {
276  return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
277 }
278 
279 
280 template <class Matrix, class Vector>
281 void
282 Cholmod<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
283 {
284  using Teuchos::RCP;
285  using Teuchos::getIntegralValue;
286  using Teuchos::ParameterEntryValidator;
287 
288  RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
289 
290 
291  if( parameterList->isParameter("Trans") ){
292  RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("Trans").validator();
293  parameterList->getEntry("Trans").setValidator(trans_validator);
294 
295  }
296 
297  if( parameterList->isParameter("IterRefine") ){
298  RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry("IterRefine").validator();
299  parameterList->getEntry("IterRefine").setValidator(refine_validator);
300 
301 
302  }
303 
304  if( parameterList->isParameter("ColPerm") ){
305  RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry("ColPerm").validator();
306  parameterList->getEntry("ColPerm").setValidator(colperm_validator);
307 
308 
309  }
310 
311  (&data_.c)->dbound = parameterList->get<double>("dbound", 0.0);
312 
313  bool prefer_upper = parameterList->get<bool>("PreferUpper", true);
314 
315  (&data_.c)->prefer_upper = prefer_upper ? 1 : 0;
316 
317  (&data_.c)->print = parameterList->get<int>("print",3);
318 
319  (&data_.c)->nmethods = parameterList->get<int>("nmethods",0);
320 
321 }
322 
323 
324 template <class Matrix, class Vector>
325 Teuchos::RCP<const Teuchos::ParameterList>
327 {
328  using std::string;
329  using Teuchos::tuple;
330  using Teuchos::ParameterList;
331  using Teuchos::EnhancedNumberValidator;
332  using Teuchos::setStringToIntegralParameter;
333  using Teuchos::stringToIntegralParameterEntryValidator;
334 
335  static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
336 
337  if( is_null(valid_params) ){
338  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
339 
340 
341  Teuchos::RCP<EnhancedNumberValidator<int> > print_validator
342  = Teuchos::rcp( new EnhancedNumberValidator<int>(0,5));
343 
344  Teuchos::RCP<EnhancedNumberValidator<int> > nmethods_validator
345  = Teuchos::rcp( new EnhancedNumberValidator<int>(0,9));
346 
347  pl->set("nmethods", 0, "Specifies the number of different ordering methods to try", nmethods_validator);
348 
349  pl->set("print", 3, "Specifies the verbosity of the print statements", print_validator);
350 
351  pl->set("dbound", 0.0,
352  "Specifies the smallest absolute value on the diagonal D for the LDL' factorization");
353 
354 
355  pl->set("Equil", true, "Whether to equilibrate the system before solve");
356 
357  pl->set("PreferUpper", true,
358  "Specifies whether the matrix will be "
359  "stored in upper triangular form.");
360 
361  valid_params = pl;
362  }
363 
364  return valid_params;
365 }
366 
367 
368 template <class Matrix, class Vector>
369 bool
371 {
372  using Teuchos::as;
373 
374 #ifdef HAVE_AMESOS2_TIMERS
375  Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
376 #endif
377 
378  // Only the root image needs storage allocated
379 
380  nzvals_.resize(this->globalNumNonZeros_);
381  rowind_.resize(this->globalNumNonZeros_);
382  colptr_.resize(this->globalNumCols_ + 1);
383 
384 
385  int nnz_ret = 0;
386  {
387 #ifdef HAVE_AMESOS2_TIMERS
388  Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
389 #endif
390 
391  TEUCHOS_TEST_FOR_EXCEPTION(this->rowIndexBase_ != this->columnIndexBase_,
392  std::runtime_error,
393  "Row and column maps have different indexbase ");
395  MatrixAdapter<Matrix>,chol_type,int,int>::do_get(this->matrixA_.ptr(),
396  nzvals_(), rowind_(),
397  colptr_(), nnz_ret, ROOTED,
398  ARBITRARY,
399  this->rowIndexBase_);
400  }
401 
402 
403  TEUCHOS_TEST_FOR_EXCEPTION(nnz_ret != as<int>(this->globalNumNonZeros_),
404  std::runtime_error,
405  "Did not get the expected number of non-zero vals");
406 
407  function_map::cholmod_init_sparse(as<size_t>(this->globalNumRows_),
408  as<size_t>(this->globalNumCols_),
409  as<size_t>(this->globalNumNonZeros_),
410  0,
411  colptr_.getRawPtr(),
412  nzvals_.getRawPtr(),
413  rowind_.getRawPtr(),
414  &(data_.A));
415 
416 
417  return true;
418 }
419 
420 
421 template<class Matrix, class Vector>
422 const char* Cholmod<Matrix,Vector>::name = "Cholmod";
423 
424 
425 } // end namespace Amesos2
426 
427 #endif // AMESOS2_CHOLMOD_DEF_HPP
~Cholmod()
Destructor.
Definition: Amesos2_Cholmod_def.hpp:91
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Cholmod_def.hpp:103
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
CHOLMOD specific solve.
Definition: Amesos2_Cholmod_def.hpp:190
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
Definition: Amesos2_TypeDecl.hpp:142
Cholmod(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Cholmod_def.hpp:68
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_Cholmod_def.hpp:282
Amesos2 CHOLMOD declarations.
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Cholmod_def.hpp:326
Definition: Amesos2_TypeDecl.hpp:127
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Cholmod_def.hpp:370
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using CHOLMOD.
Definition: Amesos2_Cholmod_def.hpp:120
Amesos2 interface to the SuperLU package.
Definition: Amesos2_Cholmod_decl.hpp:73
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
int numericFactorization_impl()
CHOLMOD specific numeric factorization.
Definition: Amesos2_Cholmod_def.hpp:142
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Cholmod_def.hpp:274