53 #ifndef AMESOS2_CHOLMOD_DEF_HPP
54 #define AMESOS2_CHOLMOD_DEF_HPP
56 #include <Teuchos_Tuple.hpp>
57 #include <Teuchos_ParameterList.hpp>
58 #include <Teuchos_StandardParameterEntryValidators.hpp>
67 template <
class Matrix,
class Vector>
69 Teuchos::RCP<const Matrix> A,
70 Teuchos::RCP<Vector> X,
71 Teuchos::RCP<const Vector> B )
78 , is_contiguous_(true)
80 CHOL::cholmod_start(&data_.c);
81 (&data_.c)->nmethods=9;
82 (&data_.c)->quick_return_if_not_posdef=1;
91 template <
class Matrix,
class Vector>
94 CHOL::cholmod_free_factor (&(data_.L), &(data_.c));
96 CHOL::cholmod_free_dense(&(data_.Y), &data_.c);
97 CHOL::cholmod_free_dense(&(data_.E), &data_.c);
99 CHOL::cholmod_finish(&(data_.c));
102 template<
class Matrix,
class Vector>
106 #ifdef HAVE_AMESOS2_TIMERS
107 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
109 data_.L = CHOL::cholmod_analyze (&data_.A, &(data_.c));
119 template <
class Matrix,
class Vector>
123 if ( !skip_symfact ){
124 #ifdef HAVE_AMESOS2_TIMERS
125 Teuchos::TimeMonitor symFactTimer(this->timers_.symFactTime_);
127 CHOL::cholmod_resymbol (&data_.A, NULL, 0,
true, data_.L, &(data_.c));
134 skip_symfact =
false;
141 template <
class Matrix,
class Vector>
149 #ifdef HAVE_AMESOS2_DEBUG
150 TEUCHOS_TEST_FOR_EXCEPTION(data_.A.ncol != as<size_t>(this->globalNumCols_),
152 "Error in converting to cholmod_sparse: wrong number of global columns." );
153 TEUCHOS_TEST_FOR_EXCEPTION(data_.A.nrow != as<size_t>(this->globalNumRows_),
155 "Error in converting to cholmod_sparse: wrong number of global rows." );
159 #ifdef HAVE_AMESOS2_TIMERS
160 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
163 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG
164 std::cout <<
"Cholmod:: Before numeric factorization" << std::endl;
165 std::cout <<
"nzvals_ : " << nzvals_.toString() << std::endl;
166 std::cout <<
"rowind_ : " << rowind_.toString() << std::endl;
167 std::cout <<
"colptr_ : " << colptr_.toString() << std::endl;
170 CHOL::cholmod_factorize(&data_.A, data_.L, &(data_.c));
172 info = (&data_.c)->status;
179 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
181 TEUCHOS_TEST_FOR_EXCEPTION(info == 2,
183 "Memory allocation failure in Cholmod factorization");
189 template <
class Matrix,
class Vector>
196 const global_size_type ld_rhs = X->getGlobalLength();
197 const size_t nrhs = X->getGlobalNumVectors();
199 Teuchos::RCP<
const Tpetra::Map<local_ordinal_type,
200 global_ordinal_type, node_type> > map2;
202 const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
203 Teuchos::Array<chol_type> xValues(val_store_size);
204 Teuchos::Array<chol_type> bValues(val_store_size);
207 #ifdef HAVE_AMESOS2_TIMERS
208 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
209 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
211 if ( is_contiguous_ ==
true ) {
213 chol_type>::do_get(B, bValues(),
215 ROOTED, this->rowIndexBase_);
219 chol_type>::do_get(B, bValues(),
221 CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
230 #ifdef HAVE_AMESOS2_TIMERS
231 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
234 function_map::cholmod_init_dense(as<int>(this->globalNumRows_),
235 as<int>(nrhs), as<int>(ld_rhs),
236 bValues.getRawPtr(), &data_.b);
237 function_map::cholmod_init_dense(as<int>(this->globalNumRows_),
238 as<int>(nrhs), as<int>(ld_rhs),
239 xValues.getRawPtr(), &data_.x);
244 #ifdef HAVE_AMESOS2_TIMERS
245 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
248 CHOL::cholmod_dense *xtemp = &(data_.x);
250 CHOL::cholmod_solve2(CHOLMOD_A, data_.L, &data_.b, NULL,
251 &xtemp, NULL, &data_.Y, &data_.E, &data_.c);
253 ierr = (&data_.c)->status;
259 Teuchos::broadcast(*(this->getComm()), 0, &ierr);
261 TEUCHOS_TEST_FOR_EXCEPTION(ierr == -2, std::runtime_error,
"Ran out of memory" );
265 #ifdef HAVE_AMESOS2_TIMERS
266 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
269 if ( is_contiguous_ ==
true ) {
273 ROOTED, this->rowIndexBase_);
279 CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
289 template <
class Matrix,
class Vector>
293 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
297 template <
class Matrix,
class Vector>
302 using Teuchos::getIntegralValue;
303 using Teuchos::ParameterEntryValidator;
305 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
308 if( parameterList->isParameter(
"Trans") ){
309 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"Trans").validator();
310 parameterList->getEntry(
"Trans").setValidator(trans_validator);
314 if( parameterList->isParameter(
"IterRefine") ){
315 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry(
"IterRefine").validator();
316 parameterList->getEntry(
"IterRefine").setValidator(refine_validator);
321 if( parameterList->isParameter(
"ColPerm") ){
322 RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry(
"ColPerm").validator();
323 parameterList->getEntry(
"ColPerm").setValidator(colperm_validator);
325 if( parameterList->isParameter(
"IsContiguous") ){
326 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
331 (&data_.c)->dbound = parameterList->get<
double>(
"dbound", 0.0);
333 bool prefer_upper = parameterList->get<
bool>(
"PreferUpper",
true);
335 (&data_.c)->prefer_upper = prefer_upper ? 1 : 0;
337 (&data_.c)->print = parameterList->get<
int>(
"print",3);
339 (&data_.c)->nmethods = parameterList->get<
int>(
"nmethods",0);
344 template <
class Matrix,
class Vector>
345 Teuchos::RCP<const Teuchos::ParameterList>
349 using Teuchos::tuple;
350 using Teuchos::ParameterList;
351 using Teuchos::EnhancedNumberValidator;
352 using Teuchos::setStringToIntegralParameter;
353 using Teuchos::stringToIntegralParameterEntryValidator;
355 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
357 if( is_null(valid_params) ){
358 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
361 Teuchos::RCP<EnhancedNumberValidator<int> > print_validator
362 = Teuchos::rcp(
new EnhancedNumberValidator<int>(0,5));
364 Teuchos::RCP<EnhancedNumberValidator<int> > nmethods_validator
365 = Teuchos::rcp(
new EnhancedNumberValidator<int>(0,9));
367 pl->set(
"nmethods", 0,
"Specifies the number of different ordering methods to try", nmethods_validator);
369 pl->set(
"print", 3,
"Specifies the verbosity of the print statements", print_validator);
371 pl->set(
"dbound", 0.0,
372 "Specifies the smallest absolute value on the diagonal D for the LDL' factorization");
375 pl->set(
"Equil",
true,
"Whether to equilibrate the system before solve");
377 pl->set(
"PreferUpper",
true,
378 "Specifies whether the matrix will be "
379 "stored in upper triangular form.");
381 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
390 template <
class Matrix,
class Vector>
396 #ifdef HAVE_AMESOS2_TIMERS
397 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
402 nzvals_.resize(this->globalNumNonZeros_);
403 rowind_.resize(this->globalNumNonZeros_);
404 colptr_.resize(this->globalNumCols_ + 1);
409 #ifdef HAVE_AMESOS2_TIMERS
410 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
413 TEUCHOS_TEST_FOR_EXCEPTION(this->rowIndexBase_ != this->columnIndexBase_,
415 "Row and column maps have different indexbase ");
416 if ( is_contiguous_ ==
true ) {
419 nzvals_(), rowind_(),
420 colptr_(), nnz_ret, ROOTED,
422 this->rowIndexBase_);
427 nzvals_(), rowind_(),
428 colptr_(), nnz_ret, CONTIGUOUS_AND_ROOTED,
430 this->rowIndexBase_);
435 TEUCHOS_TEST_FOR_EXCEPTION(nnz_ret != as<int>(this->globalNumNonZeros_),
437 "Did not get the expected number of non-zero vals");
439 function_map::cholmod_init_sparse(as<size_t>(this->globalNumRows_),
440 as<size_t>(this->globalNumCols_),
441 as<size_t>(this->globalNumNonZeros_),
453 template<
class Matrix,
class Vector>
459 #endif // AMESOS2_CHOLMOD_DEF_HPP
~Cholmod()
Destructor.
Definition: Amesos2_Cholmod_def.hpp:92
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Cholmod_def.hpp:104
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
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Cholmod_def.hpp:346
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:266
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Cholmod_def.hpp:291
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
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:654
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:191
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_Cholmod_def.hpp:299
Amesos2 CHOLMOD declarations.
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Cholmod_def.hpp:392
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using CHOLMOD.
Definition: Amesos2_Cholmod_def.hpp:121
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:322
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
int numericFactorization_impl()
CHOLMOD specific numeric factorization.
Definition: Amesos2_Cholmod_def.hpp:143