Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_CssMKL_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Amesos2: Templated Direct Sparse Solver Package
4//
5// Copyright 2011 NTESS and the Amesos2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10
19#ifndef AMESOS2_CSSMKL_DEF_HPP
20#define AMESOS2_CSSMKL_DEF_HPP
21
22#include <map>
23
24#include <Teuchos_Tuple.hpp>
25#include <Teuchos_toString.hpp>
26#include <Teuchos_StandardParameterEntryValidators.hpp>
27
28#include "Amesos2_SolverCore_def.hpp"
29#include "Amesos2_CssMKL_decl.hpp"
30
31
32namespace Amesos2 {
33
34 namespace PMKL {
35# include <mkl.h>
36# include <mkl_pardiso.h>
37 }
38
39 template <class Matrix, class Vector>
40 CssMKL<Matrix,Vector>::CssMKL(Teuchos::RCP<const Matrix> A,
41 Teuchos::RCP<Vector> X,
42 Teuchos::RCP<const Vector> B)
43 : SolverCore<Amesos2::CssMKL,Matrix,Vector>(A, X, B) // instantiate superclass
44 , n_(Teuchos::as<int_t>(this->globalNumRows_))
45 , perm_(this->globalNumRows_)
46 , nrhs_(0)
47 , css_initialized_(false)
48 , is_contiguous_(true)
49 , msglvl_(0)
50 {
51 // Matrix info
52 Teuchos::RCP<const Teuchos::Comm<int> > matComm = this->matrixA_->getComm ();
53 const global_ordinal_type indexBase = this->matrixA_->getRowMap ()->getIndexBase ();
54 const local_ordinal_type nrows = this->matrixA_->getLocalNumRows();
55
56 // rowmap for loadA (to have locally contiguous)
57 css_rowmap_ =
58 Teuchos::rcp (new map_type (this->globalNumRows_, nrows, indexBase, matComm));
59 css_contig_rowmap_ = Teuchos::rcp (new map_type (0, 0, indexBase, matComm));
60 css_contig_colmap_ = Teuchos::rcp (new map_type (0, 0, indexBase, matComm));
61
62 // set the default matrix type
63 set_css_mkl_matrix_type();
64 set_css_mkl_default_parameters(pt_, iparm_);
65
66 // index base
67 iparm_[34] = (indexBase == 0 ? 1 : 0); /* Use one or zero-based indexing */
68 // 1D block-row distribution (using Contiguous map)
69 auto frow = css_rowmap_->getMinGlobalIndex();
70 iparm_[39] = 2; /* Matrix input format. */
71 iparm_[40] = frow; /* > Beginning of input domain. */
72 iparm_[41] = frow+nrows-1; /* > End of input domain. */
73
74 // get MPI Comm
75 TEUCHOS_TEST_FOR_EXCEPTION(
76 matComm.is_null (), std::logic_error, "Amesos2::CssMKL "
77 "constructor: The matrix's communicator is null!");
78 Teuchos::RCP<const Teuchos::MpiComm<int> > matMpiComm =
79 Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> > (matComm);
80 TEUCHOS_TEST_FOR_EXCEPTION(
81 matMpiComm.is_null (), std::logic_error, "Amesos2::CssMKL "
82 "constructor: The matrix's communicator is not an MpiComm!");
83 TEUCHOS_TEST_FOR_EXCEPTION(
84 matMpiComm->getRawMpiComm ().is_null (), std::logic_error, "Amesos2::"
85 "CssMKL constructor: The matrix's communicator claims to be a "
86 "Teuchos::MpiComm<int>, but its getRawPtrComm() method returns "
87 "Teuchos::null! This means that the underlying MPI_Comm doesn't even "
88 "exist, which likely implies that the Teuchos::MpiComm was constructed "
89 "incorrectly. It means something different than if the MPI_Comm were "
90 "MPI_COMM_NULL.");
91 MPI_Comm CssComm = *(matMpiComm->getRawMpiComm ());
92 CssComm_ = MPI_Comm_c2f(CssComm);
93 }
94
95
96 template <class Matrix, class Vector>
97 CssMKL<Matrix,Vector>::~CssMKL( )
98 {
99 /*
100 * Free any memory allocated by the CssMKL library functions
101 */
102 int_t error = 0;
103 if (css_initialized_)
104 {
105 int_t phase = -1; // release all internal solver memory
106 void *bdummy, *xdummy;
107 const MPI_Fint CssComm = CssComm_;
108 function_map::cluster_sparse_solver( pt_, const_cast<int_t*>(&maxfct_),
109 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
110 nzvals_view_.data(), rowptr_view_.data(),
111 colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
112 const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &CssComm, &error );
113 css_initialized_ = false;
114 }
115 check_css_mkl_error(Amesos2::CLEAN, error);
116 }
117
118
119 template<class Matrix, class Vector>
120 int
121 CssMKL<Matrix,Vector>::preOrdering_impl()
122 {
123 // preOrdering done during "Analysis" (aka symbolic
124 // factorization) phase
125 return(0);
126 }
127
128
129 template <class Matrix, class Vector>
130 int
131 CssMKL<Matrix,Vector>::symbolicFactorization_impl()
132 {
133 if (msglvl_ > 0 && this->matrixA_->getComm()->getRank() == 0) {
134 std::cout << " CssMKL::symbolicFactorization:\n" << std::endl;
135 }
136 int_t error = 0;
137 {
138#ifdef HAVE_AMESOS2_TIMERS
139 Teuchos::TimeMonitor symbFactTimer( this->timers_.symFactTime_ );
140#endif
141
142 if (css_initialized_)
143 {
144 int_t phase = -1; // release all internal solver memory
145 void *bdummy, *xdummy;
146 const MPI_Fint CssComm = CssComm_;
147 function_map::cluster_sparse_solver( pt_, const_cast<int_t*>(&maxfct_),
148 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
149 nzvals_view_.data(), rowptr_view_.data(),
150 colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
151 const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &CssComm, &error );
152 css_initialized_ = false;
153 if (msglvl_ > 0 && error != 0 && this->matrixA_->getComm()->getRank() == 0) {
154 std::cout << " CssMKL::symbolicFactorization: clean-up failed with " << error << std::endl;
155 }
156 }
157
158 error = 0;
159 int_t phase = 11; // Analysis
160 void *bdummy, *xdummy;
161 const MPI_Fint CssComm = CssComm_;
162 function_map::cluster_sparse_solver( pt_, const_cast<int_t*>(&maxfct_),
163 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
164 nzvals_view_.data(), rowptr_view_.data(),
165 colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
166 const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &CssComm, &error );
167 }
168 check_css_mkl_error(Amesos2::SYMBFACT, error);
169 if (msglvl_ > 0 && this->matrixA_->getComm()->getRank() == 0) {
170 std::cout << " CssMKL::symbolicFactorization done:" << std::endl;
171#ifdef HAVE_AMESOS2_TIMERS
172 std::cout << " * Time : " << this->timers_.symFactTime_.totalElapsedTime() << std::endl;
173#else
174 std::cout << " * Time : not enabled" << std::endl;
175#endif
176 }
177
178 // CSS only lets you retrieve the total number of factor
179 // non-zeros, not for each individually. We should document how
180 // such a situation is reported.
181 this->setNnzLU(iparm_[17]);
182 css_initialized_ = true;
183 return(0);
184 }
185
186
187 template <class Matrix, class Vector>
188 int
189 CssMKL<Matrix,Vector>::numericFactorization_impl()
190 {
191 if (msglvl_ > 0 && this->matrixA_->getComm()->getRank() == 0) {
192 std::cout << " CssMKL::numericFactorization:\n" << std::endl;
193 }
194 int_t error = 0;
195 {
196#ifdef HAVE_AMESOS2_TIMERS
197 Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
198#endif
199
200 //int_t phase = 12; // Analysis, numerical factorization
201 int_t phase = 22; // Numerical factorization
202 void *bdummy, *xdummy;
203 const MPI_Fint CssComm = CssComm_;
204 function_map::cluster_sparse_solver( pt_, const_cast<int_t*>(&maxfct_),
205 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
206 nzvals_view_.data(), rowptr_view_.data(),
207 colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
208 const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &CssComm, &error );
209 }
210 check_css_mkl_error(Amesos2::NUMFACT, error);
211 if (msglvl_ > 0 && this->matrixA_->getComm()->getRank() == 0) {
212 std::cout << " CssMKL::numericFactorization done:" << std::endl;
213#ifdef HAVE_AMESOS2_TIMERS
214 std::cout << " * Time : " << this->timers_.numFactTime_.totalElapsedTime() << std::endl;
215#else
216 std::cout << " * Time : not enabled" << std::endl;
217#endif
218 }
219
220 return( 0 );
221 }
222
223
224 template <class Matrix, class Vector>
225 int
226 CssMKL<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
227 const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
228 {
229 using Teuchos::as;
230
231 // Get B data
232 const local_ordinal_type ld_rhs = this->matrixA_->getLocalNumRows();
233 nrhs_ = as<int_t>(X->getGlobalNumVectors());
234
235 const size_t val_store_size = as<size_t>(ld_rhs * nrhs_);
236 xvals_.resize(val_store_size);
237 bvals_.resize(val_store_size);
238 {
239#ifdef HAVE_AMESOS2_TIMERS
240 Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
241 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
242#endif
243
244 Util::get_1d_copy_helper<
245 MultiVecAdapter<Vector>,
246 solver_scalar_type>::do_get(B, bvals_(),
247 as<size_t>(ld_rhs),
248 Teuchos::ptrInArg(*css_rowmap_));
249 }
250
251 int_t error = 0;
252 {
253#ifdef HAVE_AMESOS2_TIMERS
254 Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
255#endif
256
257 const int_t phase = 33; // Solve, iterative refinement
258 const MPI_Fint CssComm = CssComm_;
259 function_map::cluster_sparse_solver( pt_,
260 const_cast<int_t*>(&maxfct_),
261 const_cast<int_t*>(&mnum_),
262 const_cast<int_t*>(&mtype_),
263 const_cast<int_t*>(&phase),
264 const_cast<int_t*>(&n_),
265 const_cast<solver_scalar_type*>(nzvals_view_.data()),
266 const_cast<int_t*>(rowptr_view_.data()),
267 const_cast<int_t*>(colind_view_.data()),
268 const_cast<int_t*>(perm_.getRawPtr()),
269 &nrhs_,
270 const_cast<int_t*>(iparm_),
271 const_cast<int_t*>(&msglvl_),
272 as<void*>(bvals_.getRawPtr()),
273 as<void*>(xvals_.getRawPtr()), &CssComm, &error );
274 }
275 check_css_mkl_error(Amesos2::SOLVE, error);
276
277 /* Get values to X */
278 {
279#ifdef HAVE_AMESOS2_TIMERS
280 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
281#endif
282
283 Util::put_1d_data_helper<
284 MultiVecAdapter<Vector>,
285 solver_scalar_type>::do_put(X, xvals_(),
286 as<size_t>(ld_rhs),
287 Teuchos::ptrInArg(*css_rowmap_));
288 }
289 if (msglvl_ > 0 && this->matrixA_->getComm()->getRank() == 0) {
290 std::cout << " CssMKL::solve done:" << std::endl;
291#ifdef HAVE_AMESOS2_TIMERS
292 std::cout << " * Time : " << this->timers_.vecRedistTime_.totalElapsedTime()
293 << " + " << this->timers_.solveTime_.totalElapsedTime() << std::endl;
294#else
295 std::cout << " * Time : not enabled" << std::endl;
296#endif
297 }
298
299 return( 0 );
300 }
301
302
303 template <class Matrix, class Vector>
304 bool
305 CssMKL<Matrix,Vector>::matrixShapeOK_impl() const
306 {
307 // CssMKL supports square matrices
308 return( this->globalNumRows_ == this->globalNumCols_ );
309 }
310
311
312 template <class Matrix, class Vector>
313 void
314 CssMKL<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
315 {
316 using Teuchos::RCP;
317 using Teuchos::getIntegralValue;
318 using Teuchos::ParameterEntryValidator;
319
320 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
321
322 // 2: Fill-in reordering from METIS, 3: thread dissection, 10: MPI version of the nested dissection
323 if( parameterList->isParameter("IPARM(2)") )
324 {
325 RCP<const ParameterEntryValidator> fillin_validator = valid_params->getEntry("IPARM(2)").validator();
326 parameterList->getEntry("IPARM(2)").setValidator(fillin_validator);
327 iparm_[1] = getIntegralValue<int>(*parameterList, "IPARM(2)");
328 }
329
330 // Max numbers of iterative refinement steps
331 if( parameterList->isParameter("IPARM(8)") )
332 {
333 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry("IPARM(8)").validator();
334 parameterList->getEntry("IPARM(8)").setValidator(refine_validator);
335 iparm_[7] = getIntegralValue<int>(*parameterList, "IPARM(8)");
336 }
337
338 // Perturb the pivot elements
339 if( parameterList->isParameter("IPARM(10)") )
340 {
341 RCP<const ParameterEntryValidator> pivot_perturb_validator = valid_params->getEntry("IPARM(10)").validator();
342 parameterList->getEntry("IPARM(10)").setValidator(pivot_perturb_validator);
343 iparm_[9] = getIntegralValue<int>(*parameterList, "IPARM(10)");
344 }
345
346 // First check if the control object requests a transpose solve.
347 // Then solver specific options can override this.
348 iparm_[11] = this->control_.useTranspose_ ? 2 : 0;
349 // Normal solve (0), or a transpose solve (1)
350 if( parameterList->isParameter("IPARM(12)") )
351 {
352 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("IPARM(12)").validator();
353 parameterList->getEntry("IPARM(12)").setValidator(trans_validator);
354 iparm_[11] = getIntegralValue<int>(*parameterList, "IPARM(12)");
355 }
356
357 // (Non-)symmetric matchings : detault 1 for nonsymmetric and 0 for symmetric matrix (default is nonsymmetric)
358 if( parameterList->isParameter("IPARM(13)") )
359 {
360 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry("IPARM(13)").validator();
361 parameterList->getEntry("IPARM(13)").setValidator(trans_validator);
362 iparm_[12] = getIntegralValue<int>(*parameterList, "IPARM(13)");
363 }
364
365 // Output: Number of nonzeros in the factor LU
366 if( parameterList->isParameter("IPARM(18)") )
367 {
368 RCP<const ParameterEntryValidator> report_validator = valid_params->getEntry("IPARM(18)").validator();
369 parameterList->getEntry("IPARM(18)").setValidator(report_validator);
370 iparm_[17] = getIntegralValue<int>(*parameterList, "IPARM(18)");
371 }
372
373 // Check input matrix is sorted
374 if( parameterList->isParameter("IPARM(27)") )
375 {
376 RCP<const ParameterEntryValidator> report_validator = valid_params->getEntry("IPARM(27)").validator();
377 parameterList->getEntry("IPARM(27)").setValidator(report_validator);
378 iparm_[26] = getIntegralValue<int>(*parameterList, "IPARM(27)");
379 }
380
381 if( parameterList->isParameter("IsContiguous") ){
382 is_contiguous_ = parameterList->get<bool>("IsContiguous");
383 }
384
385 if( parameterList->isParameter("verbose") ){
386 msglvl_ = parameterList->get<int>("verbose");
387 }
388 }
389
390
391/*
392 * TODO: It would be nice if the parameters could be expressed as
393 * either all string or as all integers. I see no way of doing this
394 * at present with the standard validators. However, we could create
395 * our own validators or kindly ask the Teuchos team to add some
396 * features for use.
397 *
398 * The issue is that with the current validators we cannot specify
399 * arbitrary sets of numbers that are the only allowed parameters.
400 * For example the IPARM(2) parameter can take only the values 0, 2,
401 * and 3. The EnhancedNumberValidator can take a min value, and max
402 * value, and a step size, but with those options there is no way to
403 * specify the needed set.
404 *
405 * Another missing feature is the ability to give docstrings for such
406 * numbers. For example IPARM(25) can take on the values 0 and 1.
407 * This would be easy enough to accomplish with just a number
408 * validator, but then have no way to document the effect of each
409 * value.
410 */
411template <class Matrix, class Vector>
412Teuchos::RCP<const Teuchos::ParameterList>
413CssMKL<Matrix,Vector>::getValidParameters_impl() const
414{
415 using std::string;
416 using Teuchos::as;
417 using Teuchos::RCP;
418 using Teuchos::tuple;
419 using Teuchos::toString;
420 using Teuchos::EnhancedNumberValidator;
421 using Teuchos::setStringToIntegralParameter;
422 using Teuchos::anyNumberParameterEntryValidator;
423
424 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
425
426 if( is_null(valid_params) ){
427 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
428
429 void* pt_temp[64];
430 int_t iparm_temp[64];
431 set_css_mkl_default_parameters(pt_temp, iparm_temp);
432 setStringToIntegralParameter<int>("IPARM(2)", toString(iparm_temp[1]),
433 "Fill-in reducing ordering for the input matrix",
434 tuple<string>("2", "3", "10"),
435 tuple<string>("Nested dissection algorithm from METIS",
436 "Parallel version of the nested dissection algorithm",
437 "MPI version of the nested dissection and symbolic factorization algorithms"),
438 tuple<int>(2, 3, 10),
439 pl.getRawPtr());
440
441 setStringToIntegralParameter<int>("IPARM(12)", toString(iparm_temp[11]),
442 "Solve with transposed or conjugate transposed matrix A",
443 tuple<string>("0", "1", "2"),
444 tuple<string>("Non-transposed",
445 "Conjugate-transposed",
446 "Transposed"),
447 tuple<int>(0, 1, 2),
448 pl.getRawPtr());
449
450 setStringToIntegralParameter<int>("IPARM(13)", toString(iparm_temp[12]),
451 "Use weighted matching",
452 tuple<string>("0", "1"),
453 tuple<string>("No matching", "Use matching"),
454 tuple<int>(0, 1),
455 pl.getRawPtr());
456
457 Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
458 Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
459
460 Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int( false );
461 accept_int.allowInt( true );
462
463 pl->set("IPARM(8)" , as<int>(iparm_temp[7]) , "Iterative refinement step",
464 anyNumberParameterEntryValidator(preferred_int, accept_int));
465
466 pl->set("IPARM(10)", as<int>(iparm_temp[9]) , "Pivoting perturbation",
467 anyNumberParameterEntryValidator(preferred_int, accept_int));
468
469 pl->set("IPARM(18)", as<int>(iparm_temp[17]), "Report the number of non-zero elements in the factors",
470 anyNumberParameterEntryValidator(preferred_int, accept_int));
471
472 pl->set("IPARM(28)", as<int>(iparm_temp[27]), "Check input matrix is sorted",
473 anyNumberParameterEntryValidator(preferred_int, accept_int));
474
475 pl->set("IsContiguous", true, "Whether GIDs contiguous");
476
477 pl->set("verbose", 0, "Verbosity Message Level");
478
479 valid_params = pl;
480 }
481
482 return valid_params;
483}
484
485
486
487template <class Matrix, class Vector>
488bool
489CssMKL<Matrix,Vector>::loadA_impl(EPhase current_phase)
490{
491#ifdef HAVE_AMESOS2_TIMERS
492 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
493#endif
494
495 // CssMKL does not need matrix data in the pre-ordering phase
496 if( current_phase == PREORDERING ) return( false );
497
498 // is_contiguous : input is contiguous
499 // CONTIGUOUS_AND_ROOTED : input is not contiguous, so make output contiguous
500 EDistribution dist_option = (iparm_[39] != 0 ? DISTRIBUTED_NO_OVERLAP : ((is_contiguous_ == true) ? ROOTED : CONTIGUOUS_AND_ROOTED));
501 if (dist_option == DISTRIBUTED_NO_OVERLAP && !is_contiguous_) {
502 // Neeed to form contiguous matrix
503 #if 1
504 // Only reinex GIDs
505 css_rowmap_ = this->matrixA_->getRowMap(); // use original map to redistribute vectors in solve
506 Teuchos::RCP<const MatrixAdapter<Matrix> > contig_mat = this->matrixA_->reindex(css_contig_rowmap_, css_contig_colmap_, current_phase);
507 #else
508 // Redistribued matrixA into contiguous GIDs
509 Teuchos::RCP<const MatrixAdapter<Matrix> > contig_mat = this->matrixA_->get(ptrInArg(*css_rowmap_));
510 //css_rowmap_ = contig_mat->getRowMap(); // use new map to redistribute vectors in solve
511 #endif
512 // Copy into local views
513 if (current_phase == SYMBFACT) {
514 Kokkos::resize(nzvals_temp_, contig_mat->getLocalNNZ());
515 Kokkos::resize(nzvals_view_, contig_mat->getLocalNNZ());
516 Kokkos::resize(colind_view_, contig_mat->getLocalNNZ());
517 Kokkos::resize(rowptr_view_, contig_mat->getLocalNumRows() + 1);
518 }
519 int_t nnz_ret = 0;
520 {
521#ifdef HAVE_AMESOS2_TIMERS
522 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
523#endif
524 Util::get_crs_helper_kokkos_view<MatrixAdapter<Matrix>,
525 host_value_type_array,host_ordinal_type_array, host_size_type_array >::do_get(
526 contig_mat.ptr(),
527 nzvals_temp_, colind_view_, rowptr_view_,
528 nnz_ret,
529 ptrInArg(*(contig_mat->getRowMap())),
530 #if 1
531 DISTRIBUTED_NO_OVERLAP,
532 #else
533 ROOTED,
534 #endif
535 SORTED_INDICES);
536 Kokkos::deep_copy(nzvals_view_, nzvals_temp_);
537 }
538 } else {
539 if (current_phase == SYMBFACT) {
540 if (dist_option == DISTRIBUTED_NO_OVERLAP) {
541 Kokkos::resize(nzvals_temp_, this->matrixA_->getLocalNNZ());
542 Kokkos::resize(nzvals_view_, this->matrixA_->getLocalNNZ());
543 Kokkos::resize(colind_view_, this->matrixA_->getLocalNNZ());
544 Kokkos::resize(rowptr_view_, this->matrixA_->getLocalNumRows() + 1);
545 } else {
546 if( this->root_ ) {
547 Kokkos::resize(nzvals_temp_, this->matrixA_->getGlobalNNZ());
548 Kokkos::resize(nzvals_view_, this->matrixA_->getGlobalNNZ());
549 Kokkos::resize(colind_view_, this->matrixA_->getGlobalNNZ());
550 Kokkos::resize(rowptr_view_, this->matrixA_->getGlobalNumRows() + 1);
551 } else {
552 Kokkos::resize(nzvals_temp_, 0);
553 Kokkos::resize(nzvals_view_, 0);
554 Kokkos::resize(colind_view_, 0);
555 Kokkos::resize(rowptr_view_, 0);
556 }
557 }
558 }
559 int_t nnz_ret = 0;
560 {
561#ifdef HAVE_AMESOS2_TIMERS
562 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
563#endif
564 Util::get_crs_helper_kokkos_view<MatrixAdapter<Matrix>,
565 host_value_type_array,host_ordinal_type_array, host_size_type_array >::do_get(
566 this->matrixA_.ptr(),
567 nzvals_temp_, colind_view_, rowptr_view_,
568 nnz_ret,
569 ptrInArg(*(this->matrixA_->getRowMap())),
570 dist_option,
571 SORTED_INDICES);
572 Kokkos::deep_copy(nzvals_view_, nzvals_temp_);
573 }
574 }
575 return( true );
576}
577
578
579template <class Matrix, class Vector>
580void
581CssMKL<Matrix,Vector>::describe_impl(Teuchos::FancyOStream &out,
582 const Teuchos::EVerbosityLevel verbLevel) const
583{
584 out << " CssMKL current parameters:" << std::endl;
585 out << " > IPARM(2) = " << iparm_[1] << std::endl;
586 out << " > IPARM(8) = " << iparm_[7] << std::endl;
587 out << " > IPARM(10) = " << iparm_[9] << std::endl;
588 out << " > IPARM(12) = " << iparm_[11] << std::endl;
589 out << " > IPARM(13) = " << iparm_[12] << std::endl;
590 out << " > IPARM(18) = " << iparm_[17] << std::endl;
591 out << " > IPARM(27) = " << iparm_[26] << std::endl;
592 out << " > IsContiguous = " << (is_contiguous_ ? "YES" : "NO") << std::endl;
593 out << " > verbose = " << msglvl_ << std::endl;
594 out << std::endl;
595}
596
597
598template <class Matrix, class Vector>
599void
600CssMKL<Matrix,Vector>::check_css_mkl_error(EPhase phase,
601 int_t error) const
602{
603 int error_i = error;
604 Teuchos::broadcast(*(this->getComm()), 0, &error_i); // We only care about root's value
605
606 if( error == 0 ) return; // No error
607
608 std::string errmsg = "Other error";
609 switch( error ){
610 case -1:
611 errmsg = "CssMKL reported error: 'Input inconsistent'";
612 break;
613 case -2:
614 errmsg = "CssMKL reported error: 'Not enough memory'";
615 break;
616 case -3:
617 errmsg = "CssMKL reported error: 'Reordering problem'";
618 break;
619 case -4:
620 errmsg =
621 "CssMKL reported error: 'Zero pivot, numerical "
622 "factorization or iterative refinement problem'";
623 break;
624 case -5:
625 errmsg = "CssMKL reported error: 'Unclassified (internal) error'";
626 break;
627 case -6:
628 errmsg = "CssMKL reported error: 'Reordering failed'";
629 break;
630 case -7:
631 errmsg = "CssMKL reported error: 'Diagonal matrix is singular'";
632 break;
633 case -8:
634 errmsg = "CssMKL reported error: '32-bit integer overflow problem'";
635 break;
636 case -9:
637 errmsg = "CssMKL reported error: 'Not enough memory for OOC'";
638 break;
639 case -10:
640 errmsg = "CssMKL reported error: 'Problems with opening OOC temporary files'";
641 break;
642 case -11:
643 errmsg = "CssMKL reported error: 'Read/write problem with OOC data file'";
644 break;
645 }
646 errmsg += (" at phase = "+std::to_string(phase));
647
648 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error, errmsg );
649}
650
651
652template <class Matrix, class Vector>
653void
654CssMKL<Matrix,Vector>::set_css_mkl_matrix_type(int_t mtype)
655{
656 if( mtype == 0 ){
657 if( complex_ ){
658 mtype_ = 13; // complex, unsymmetric
659 } else {
660 mtype_ = 11; // real, unsymmetric
661 }
662 } else {
663 switch( mtype ){
664 case 11:
665 TEUCHOS_TEST_FOR_EXCEPTION( complex_,
666 std::invalid_argument,
667 "Cannot set a real CSS matrix type with scalar type complex" );
668 mtype_ = 11; break;
669 case 13:
670 TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
671 std::invalid_argument,
672 "Cannot set a complex CSS matrix type with non-complex scalars" );
673 mtype_ = 13; break;
674 default:
675 TEUCHOS_TEST_FOR_EXCEPTION( true,
676 std::invalid_argument,
677 "Symmetric matrices are not yet supported by the Amesos2 interface" );
678 }
679 }
680}
681
682template <class Matrix, class Vector>
683void
684CssMKL<Matrix,Vector>::set_css_mkl_default_parameters(void* pt[], int_t iparm[]) const
685{
686 for( int i = 0; i < 64; ++i ){
687 pt[i] = nullptr;
688 iparm[i] = 0;
689 }
690 iparm[0] = 1; /* No solver default */
691 // Reset some of the default parameters
692 iparm[1] = 10; /* 2: Fill-in reordering from METIS, 3: thread dissection, 10: MPI version of the nested dissection and symbolic factorization*/
693 iparm[7] = 0; /* Max numbers of iterative refinement steps */
694 iparm[10] = 0; /* Disable nonsymmetric permutation and scaling MPS */
695 iparm[11] = 0; /* Normal solve (0), or a transpose solve (1) */
696 iparm[12] = 0; /* Do not use (non-)symmetric matchings */
697 iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
698 iparm[20] = 1; /* Pivoting for symmetric indefinite matrices */
699 iparm[26] = 1; /* Check input matrix is sorted */
700
701 // diagonal pertubation
702 if (mtype_ == -2 || mtype_ == -4) {
703 // symmetric indefinite
704 iparm[9] = 8; /* Perturb the pivot elements with 1E-8 */
705 } else {
706 // non-symmetric
707 iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
708 }
709
710 // set single or double precision
711 if constexpr ( std::is_same_v<solver_magnitude_type, PMKL::_REAL_t> ) {
712 iparm[27] = 1; // single-precision
713 } else {
714 iparm[27] = 0; // double-precision
715 }
716 iparm[34] = 1; /* Use zero-based indexing */
717}
718
719
720template <class Matrix, class Vector>
721const char* CssMKL<Matrix,Vector>::name = "CSSMKL";
722
723template <class Matrix, class Vector>
724const typename CssMKL<Matrix,Vector>::int_t
725CssMKL<Matrix,Vector>::maxfct_ = 1;
726
727template <class Matrix, class Vector>
728const typename CssMKL<Matrix,Vector>::int_t
729CssMKL<Matrix,Vector>::mnum_ = 1;
730
731
732} // end namespace Amesos
733
734#endif // AMESOS2_CSSMKL_DEF_HPP
Amesos2_CssMKL_decl.hpp
A template class that does nothing useful besides show developers what, in general,...
Amesos2_SolverCore_def.hpp
Amesos2::EDistribution
EDistribution
Definition Amesos2_TypeDecl.hpp:89
Amesos2::DISTRIBUTED_NO_OVERLAP
@ DISTRIBUTED_NO_OVERLAP
Definition Amesos2_TypeDecl.hpp:91
Amesos2::ROOTED
@ ROOTED
Definition Amesos2_TypeDecl.hpp:93
Amesos2::CONTIGUOUS_AND_ROOTED
@ CONTIGUOUS_AND_ROOTED
Definition Amesos2_TypeDecl.hpp:94
Amesos2::SORTED_INDICES
@ SORTED_INDICES
Definition Amesos2_TypeDecl.hpp:108
Amesos2::CssMKL
Amesos2 interface to the CssMKL package.
Definition Amesos2_CssMKL_decl.hpp:50
Amesos2::CssMKL::pt_
void * pt_[64]
CssMKL internal data address pointer.
Definition Amesos2_CssMKL_decl.hpp:282
Amesos2::CssMKL::loadA_impl
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition Amesos2_CssMKL_def.hpp:489
Amesos2::CssMKL::~CssMKL
~CssMKL()
Destructor.
Definition Amesos2_CssMKL_def.hpp:97
Amesos2::CssMKL::check_css_mkl_error
void check_css_mkl_error(EPhase phase, int_t error) const
Throws an appropriate runtime error in the event that error < 0 .
Definition Amesos2_CssMKL_def.hpp:600
Amesos2::CssMKL::name
static const char * name
The name of this solver interface.
Definition Amesos2_CssMKL_decl.hpp:57
Amesos2::CssMKL::preOrdering_impl
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition Amesos2_CssMKL_def.hpp:121
Amesos2::CssMKL::numericFactorization_impl
int numericFactorization_impl()
CssMKL specific numeric factorization.
Definition Amesos2_CssMKL_def.hpp:189
Amesos2::CssMKL::symbolicFactorization_impl
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using CssMKL.
Definition Amesos2_CssMKL_def.hpp:131
Amesos2::CssMKL::matrixShapeOK_impl
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition Amesos2_CssMKL_def.hpp:305
Amesos2::CssMKL::solve_impl
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
CssMKL specific solve.
Definition Amesos2_CssMKL_def.hpp:226
Amesos2::CssMKL::getValidParameters_impl
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition Amesos2_CssMKL_def.hpp:413
Amesos2::CssMKL::iparm_
int_t iparm_[64]
Definition Amesos2_CssMKL_decl.hpp:300
Amesos2::CssMKL::set_css_mkl_matrix_type
void set_css_mkl_matrix_type(int_t mtype=0)
Definition Amesos2_CssMKL_def.hpp:654
Amesos2::CssMKL::setParameters_impl
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition Amesos2_CssMKL_def.hpp:314
Amesos2::CssMKL::CssMKL
CssMKL(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition Amesos2_CssMKL_def.hpp:40
Amesos2::CssMKL::describe_impl
void describe_impl(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the status information about the current solver with some level of verbosity.
Definition Amesos2_CssMKL_def.hpp:581
Amesos2::SolverCore
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition Amesos2_SolverCore_decl.hpp:72
Amesos2::SolverCore< Amesos2::CssMKL, Matrix, Vector >::matrixA_
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition Amesos2_SolverCore_decl.hpp:421
Amesos2::SolverCore::getComm
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const override
Returns a pointer to the Teuchos::Comm communicator with this operator.
Definition Amesos2_SolverCore_decl.hpp:329
Amesos2::SolverCore< Amesos2::CssMKL, Matrix, Vector >::globalNumRows_
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition Amesos2_SolverCore_decl.hpp:442
Amesos2::EPhase
EPhase
Used to indicate a phase in the direct solution.
Definition Amesos2_TypeDecl.hpp:31
Amesos2::MultiVecAdapter
A templated MultiVector class adapter for Amesos2.
Definition Amesos2_MultiVecAdapter_decl.hpp:142
Amesos2::Util::get_1d_copy_helper
Helper class for getting 1-D copies of multivectors.
Definition Amesos2_MultiVecAdapter_decl.hpp:233
Amesos2::Util::get_crs_helper_kokkos_view
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition Amesos2_Util.hpp:600
Amesos2::Util::put_1d_data_helper
Helper class for putting 1-D data arrays into multivectors.
Definition Amesos2_MultiVecAdapter_decl.hpp:339
Generated by doxygen 1.9.8