Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_SolverCore_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 
19 #ifndef AMESOS2_SOLVERCORE_DEF_HPP
20 #define AMESOS2_SOLVERCORE_DEF_HPP
21 
22 #if KOKKOS_VERSION >= 40799
23 #include "KokkosKernels_ArithTraits.hpp"
24 #else
25 #include "Kokkos_ArithTraits.hpp"
26 #endif
27 
28 #include "Amesos2_MatrixAdapter_def.hpp"
29 #include "Amesos2_MultiVecAdapter_def.hpp"
30 
31 #include "Amesos2_Util.hpp"
32 
33 #include "KokkosSparse_spmv.hpp"
34 #include "KokkosBlas.hpp"
35 
36 namespace Amesos2 {
37 
38 
39 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
40 SolverCore<ConcreteSolver,Matrix,Vector>::SolverCore(
41  Teuchos::RCP<const Matrix> A,
42  Teuchos::RCP<Vector> X,
43  Teuchos::RCP<const Vector> B )
44  : matrixA_(createConstMatrixAdapter<Matrix>(A))
45  , multiVecX_(X) // may be null
46  , multiVecB_(B) // may be null
47  , globalNumRows_(matrixA_->getGlobalNumRows())
48  , globalNumCols_(matrixA_->getGlobalNumCols())
49  , globalNumNonZeros_(matrixA_->getGlobalNNZ())
50  , rowIndexBase_(matrixA_->getRowIndexBase())
51  , columnIndexBase_(matrixA_->getColumnIndexBase())
52  , rank_(Teuchos::rank(*this->getComm()))
53  , root_(rank_ == 0)
54  , nprocs_(Teuchos::size(*this->getComm()))
55 {
56  TEUCHOS_TEST_FOR_EXCEPTION(
57  !matrixShapeOK(),
58  std::invalid_argument,
59  "Matrix shape inappropriate for this solver");
60 }
61 
62 
64 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
65 SolverCore<ConcreteSolver,Matrix,Vector>::~SolverCore( )
66 {
67  // Nothing
68 }
69 
70 
71 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
72 Solver<Matrix,Vector>&
73 SolverCore<ConcreteSolver,Matrix,Vector>::preOrdering()
74 {
75 #ifdef HAVE_AMESOS2_TIMERS
76  Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
77 #endif
78 
79  loadA(PREORDERING);
80 
81  int error_code = static_cast<solver_type*>(this)->preOrdering_impl();
82  if (error_code == EXIT_SUCCESS){
83  ++status_.numPreOrder_;
84  status_.last_phase_ = PREORDERING;
85  }
86 
87  return *this;
88 }
89 
90 
91 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
92 Solver<Matrix,Vector>&
93 SolverCore<ConcreteSolver,Matrix,Vector>::symbolicFactorization()
94 {
95 #ifdef HAVE_AMESOS2_TIMERS
96  Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
97 #endif
98 
99  {
100 #ifdef HAVE_AMESOS2_TIMERS
101  Teuchos::TimeMonitor LocalTimer2(timers_.coreSymFactTime_);
102 #endif
103  if( !status_.preOrderingDone() ){
104  preOrdering();
105  if( !matrix_loaded_ ) loadA(SYMBFACT);
106  } else {
107  loadA(SYMBFACT);
108  }
109 
110  int error_code = static_cast<solver_type*>(this)->symbolicFactorization_impl();
111  if (error_code == EXIT_SUCCESS){
112  ++status_.numSymbolicFact_;
113  status_.last_phase_ = SYMBFACT;
114  }
115  }
116  return *this;
117 }
118 
119 
120 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
121 Solver<Matrix,Vector>&
122 SolverCore<ConcreteSolver,Matrix,Vector>::numericFactorization()
123 {
124 #ifdef HAVE_AMESOS2_TIMERS
125  Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
126 #endif
127  {
128 #ifdef HAVE_AMESOS2_TIMERS
129  Teuchos::TimeMonitor LocalTimer2(timers_.coreNumFactTime_);
130 #endif
131  if( !status_.symbolicFactorizationDone() ){
132  symbolicFactorization();
133  if( !matrix_loaded_ ) loadA(NUMFACT);
134  } else {
135  loadA(NUMFACT);
136  }
137 
138  int error_code = static_cast<solver_type*>(this)->numericFactorization_impl();
139  if (error_code == EXIT_SUCCESS){
140  ++status_.numNumericFact_;
141  status_.last_phase_ = NUMFACT;
142  }
143  }
144 
145  return *this;
146 }
147 
148 
149 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
150 void
151 SolverCore<ConcreteSolver,Matrix,Vector>::solve()
152 {
153  solve(multiVecX_.ptr(), multiVecB_.ptr());
154 }
155 
156 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
157 void
158 SolverCore<ConcreteSolver,Matrix,Vector>::solve(const Teuchos::Ptr<Vector> X,
159  const Teuchos::Ptr<const Vector> B) const
160 {
161 #ifdef HAVE_AMESOS2_TIMERS
162  Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
163 #endif
164 
165  X.assert_not_null();
166  B.assert_not_null();
167 
168  if (control_.useIterRefine_) {
169  solve_ir(X, B, control_.maxNumIterRefines_, control_.verboseIterRefine_);
170  return;
171  }
172 
173  const Teuchos::RCP<MultiVecAdapter<Vector> > x =
174  createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(X));
175  const Teuchos::RCP<const MultiVecAdapter<Vector> > b =
176  createConstMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(B));
177 
178 #ifdef HAVE_AMESOS2_DEBUG
179  // Check some required properties of X and B
180  TEUCHOS_TEST_FOR_EXCEPTION
181  (x->getGlobalLength() != matrixA_->getGlobalNumCols(),
182  std::invalid_argument,
183  "MultiVector X must have length equal to the number of "
184  "global columns in A. X->getGlobalLength() = "
185  << x->getGlobalLength() << " != A->getGlobalNumCols() = "
186  << matrixA_->getGlobalNumCols() << ".");
187 
188  TEUCHOS_TEST_FOR_EXCEPTION(b->getGlobalLength() != matrixA_->getGlobalNumRows(),
189  std::invalid_argument,
190  "MultiVector B must have length equal to the number of "
191  "global rows in A");
192 
193  TEUCHOS_TEST_FOR_EXCEPTION(x->getGlobalNumVectors() != b->getGlobalNumVectors(),
194  std::invalid_argument,
195  "X and B MultiVectors must have the same number of vectors");
196 #endif // HAVE_AMESOS2_DEBUG
197 
198  if( !status_.numericFactorizationDone() ){
199  // This casting-away of constness is probably OK because this
200  // function is meant to be "logically const"
201  const_cast<type&>(*this).numericFactorization();
202  }
203 
204  {
205 #ifdef HAVE_AMESOS2_TIMERS
206  Teuchos::TimeMonitor LocalTimer2(timers_.coreSolveTime_);
207 #endif
208  int error_code = static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*x), Teuchos::ptrInArg(*b));
209  if (error_code == EXIT_SUCCESS){
210  ++status_.numSolve_;
211  status_.last_phase_ = SOLVE;
212  }
213  }
214 }
215 
216 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
217 void
218 SolverCore<ConcreteSolver,Matrix,Vector>::solve(Vector* X, const Vector* B) const
219 {
220  solve(Teuchos::ptr(X), Teuchos::ptr(B));
221 }
222 
223 
224 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
225 int
226 SolverCore<ConcreteSolver,Matrix,Vector>::solve_ir(const int maxNumIters, const bool verbose)
227 {
228  return solve_ir(multiVecX_.ptr(), multiVecB_.ptr(), maxNumIters, verbose);
229 }
230 
231 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
232 int
233 SolverCore<ConcreteSolver,Matrix,Vector>::solve_ir(Vector* X, const Vector* B, const int maxNumIters, const bool verbose) const
234 {
235  return solve_ir(Teuchos::ptr(X), Teuchos::ptr(B), maxNumIters, verbose);
236 }
237 
238 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
239 int
240 SolverCore<ConcreteSolver,Matrix,Vector>::solve_ir(const Teuchos::Ptr< Vector> x,
241  const Teuchos::Ptr<const Vector> b,
242  const int maxNumIters,
243  const bool verbose) const
244 {
245 #if KOKKOS_VERSION >= 40799
246  using KAT = KokkosKernels::ArithTraits<scalar_type>;
247 #else
248  using KAT = Kokkos::ArithTraits<scalar_type>;
249 #endif
250  using impl_scalar_type = typename KAT::val_type;
251  using magni_type = typename KAT::mag_type;
252  using host_execution_space = Kokkos::DefaultHostExecutionSpace;
253  using host_crsmat_t = KokkosSparse::CrsMatrix<impl_scalar_type, int, host_execution_space, void, int>;
254  using host_graph_t = typename host_crsmat_t::StaticCrsGraphType;
255  using host_values_t = typename host_crsmat_t::values_type::non_const_type;
256  using host_row_map_t = typename host_graph_t::row_map_type::non_const_type;
257  using host_colinds_t = typename host_graph_t::entries_type::non_const_type;
258  using host_mvector_t = Kokkos::View<impl_scalar_type **, Kokkos::LayoutLeft, host_execution_space>;
259  using host_vector_t = Kokkos::View<impl_scalar_type *, Kokkos::LayoutLeft, host_execution_space>;
260  using host_magni_view = Kokkos::View<magni_type *, Kokkos::LayoutLeft, host_execution_space>;
261 
262  const impl_scalar_type one(1.0);
263  const impl_scalar_type mone = impl_scalar_type(-one);
264  const magni_type eps = KAT::eps ();
265 
266  // get data needed for IR
267  using MVAdapter = MultiVecAdapter<Vector>;
268  Teuchos::RCP< MVAdapter> X = createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(x));
269  Teuchos::RCP<const MVAdapter> B = createConstMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(b));
270 
271  auto r_ = B->clone();
272  auto e_ = X->clone();
273  auto r = r_.ptr();
274  auto e = e_.ptr();
275  Teuchos::RCP< MVAdapter> R = createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(r));
276  Teuchos::RCP< MVAdapter> E = createMultiVecAdapter<Vector>(Teuchos::rcpFromPtr(e));
277 
278  const size_t nrhs = X->getGlobalNumVectors();
279  const int nnz = this->matrixA_->getGlobalNNZ();
280  const int nrows = this->matrixA_->getGlobalNumRows();
281 
282  // get local matriix
283  host_crsmat_t crsmat;
284  host_graph_t static_graph;
285  host_row_map_t rowmap_view;
286  host_colinds_t colind_view;
287  host_values_t values_view;
288  if (this->root_) {
289  Kokkos::resize(rowmap_view, 1+nrows);
290  Kokkos::resize(colind_view, nnz);
291  Kokkos::resize(values_view, nnz);
292  } else {
293  Kokkos::resize(rowmap_view, 1);
294  Kokkos::resize(colind_view, 0);
295  Kokkos::resize(values_view, 0);
296  }
297 
298  int nnz_ret = 0;
299  Util::get_crs_helper_kokkos_view<
300  MatrixAdapter<Matrix>, host_values_t, host_colinds_t, host_row_map_t>::do_get(
301  this->matrixA_.ptr(),
302  values_view, colind_view, rowmap_view,
303  nnz_ret, ROOTED, ARBITRARY, this->rowIndexBase_);
304 
305  if (this->root_) {
306  static_graph = host_graph_t(colind_view, rowmap_view);
307  crsmat = host_crsmat_t("CrsMatrix", nrows, values_view, static_graph);
308  }
309 
310  //
311  // ** First Solve **
312  static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*X), Teuchos::ptrInArg(*B));
313 
314 
315  // auxiliary scalar Kokkos views
316  const int ldx = (this->root_ ? X->getGlobalLength() : 0);
317  const int ldb = (this->root_ ? B->getGlobalLength() : 0);
318  const int ldr = (this->root_ ? R->getGlobalLength() : 0);
319  const int lde = (this->root_ ? E->getGlobalLength() : 0);
320  const bool initialize_data = true;
321  const bool not_initialize_data = true;
322  host_mvector_t X_view;
323  host_mvector_t B_view;
324  host_mvector_t R_view;
325  host_mvector_t E_view;
326 
327  global_size_type rowIndexBase = this->rowIndexBase_;
328  auto Xptr = Teuchos::Ptr< MVAdapter>(X.ptr());
329  auto Bptr = Teuchos::Ptr<const MVAdapter>(B.ptr());
330  auto Rptr = Teuchos::Ptr< MVAdapter>(R.ptr());
331  auto Eptr = Teuchos::Ptr< MVAdapter>(E.ptr());
332  Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
333  do_get( initialize_data, Xptr, X_view, ldx, CONTIGUOUS_AND_ROOTED, rowIndexBase);
334  Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
335  do_get( initialize_data, Bptr, B_view, ldb, CONTIGUOUS_AND_ROOTED, rowIndexBase);
336  Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
337  do_get(not_initialize_data, Rptr, R_view, ldr, CONTIGUOUS_AND_ROOTED, rowIndexBase);
338  Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
339  do_get(not_initialize_data, Eptr, E_view, lde, CONTIGUOUS_AND_ROOTED, rowIndexBase);
340 
341 
342  host_magni_view x0norms("x0norms", nrhs);
343  host_magni_view bnorms("bnorms", nrhs);
344  host_magni_view enorms("enorms", nrhs);
345  if (this->root_) {
346  // compute initial solution norms (used for stopping criteria)
347  for (size_t j = 0; j < nrhs; j++) {
348  auto x_subview = Kokkos::subview(X_view, Kokkos::ALL(), j);
349  host_vector_t x_1d (const_cast<impl_scalar_type*>(x_subview.data()), x_subview.extent(0));
350  x0norms(j) = KokkosBlas::nrm2(x_1d);
351  }
352  if (verbose) {
353  std::cout << std::endl
354  << " SolverCore :: solve_ir (maxNumIters = " << maxNumIters
355  << ", tol = " << x0norms(0) << " * " << eps << " = " << x0norms(0)*eps
356  << ")" << std::endl;
357  }
358 
359  // compute residual norm
360  if (verbose) {
361  std::cout << " bnorm = ";
362  for (size_t j = 0; j < nrhs; j++) {
363  auto b_subview = Kokkos::subview(B_view, Kokkos::ALL(), j);
364  host_vector_t b_1d (const_cast<impl_scalar_type*>(b_subview.data()), b_subview.extent(0));
365  bnorms(j) = KokkosBlas::nrm2(b_1d);
366  std::cout << bnorms(j) << ", ";
367  }
368  std::cout << std::endl;
369  }
370  }
371 
372 
373  //
374  // ** Iterative Refinement **
375  int numIters = 0;
376  int converged = 0; // 0 = has not converged, 1 = converged
377  for (numIters = 0; numIters < maxNumIters && converged == 0; ++numIters) {
378  // r = b - Ax (on rank-0)
379  if (this->root_) {
380  Kokkos::deep_copy(R_view, B_view);
381  KokkosSparse::spmv("N", mone, crsmat, X_view, one, R_view);
382  Kokkos::fence();
383 
384  if (verbose) {
385  // compute residual norm
386  std::cout << " > " << numIters << " : norm(r,x,e) = ";
387  for (size_t j = 0; j < nrhs; j++) {
388  auto r_subview = Kokkos::subview(R_view, Kokkos::ALL(), j);
389  auto x_subview = Kokkos::subview(X_view, Kokkos::ALL(), j);
390  host_vector_t r_1d (const_cast<impl_scalar_type*>(r_subview.data()), r_subview.extent(0));
391  host_vector_t x_1d (const_cast<impl_scalar_type*>(x_subview.data()), x_subview.extent(0));
392  impl_scalar_type rnorm = KokkosBlas::nrm2(r_1d);
393  impl_scalar_type xnorm = KokkosBlas::nrm2(x_1d);
394  std::cout << rnorm << " -> " << rnorm/bnorms(j) << " " << xnorm << " " << enorms(j) << ", ";
395  }
396  std::cout << std::endl;
397  }
398  }
399 
400  // e = A^{-1} r
401  Util::put_1d_data_helper_kokkos_view<MVAdapter, host_mvector_t>::
402  do_put(Rptr, R_view, ldr, CONTIGUOUS_AND_ROOTED, rowIndexBase);
403  static_cast<const solver_type*>(this)->solve_impl(Teuchos::outArg(*E), Teuchos::ptrInArg(*R));
404  Util::get_1d_copy_helper_kokkos_view<MVAdapter, host_mvector_t>::
405  do_get(initialize_data, Eptr, E_view, lde, CONTIGUOUS_AND_ROOTED, rowIndexBase);
406 
407  // x = x + e (on rank-0)
408  if (this->root_) {
409  KokkosBlas::axpy(one, E_view, X_view);
410 
411  if (numIters < maxNumIters-1) {
412  // compute norm of corrections for "convergence" check
413  converged = 1;
414  for (size_t j = 0; j < nrhs; j++) {
415  auto e_subview = Kokkos::subview(E_view, Kokkos::ALL(), j);
416  host_vector_t e_1d (const_cast<impl_scalar_type*>(e_subview.data()), e_subview.extent(0));
417  enorms(j) = KokkosBlas::nrm2(e_1d);
418  if (enorms(j) > eps * x0norms(j)) {
419  converged = 0;
420  }
421  }
422  if (verbose && converged) {
423  std::cout << " converged " << std::endl;
424  }
425  }
426  }
427 
428  // broadcast "converged"
429  Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &converged);
430  } // end of for-loop for IR iteration
431 
432  if (verbose && this->root_) {
433  // r = b - Ax
434  Kokkos::deep_copy(R_view, B_view);
435  KokkosSparse::spmv("N", mone, crsmat, X_view, one, R_view);
436  Kokkos::fence();
437  std::cout << " > final residual norm = ";
438  for (size_t j = 0; j < nrhs; j++) {
439  auto r_subview = Kokkos::subview(R_view, Kokkos::ALL(), j);
440  host_vector_t r_1d (const_cast<impl_scalar_type*>(r_subview.data()), r_subview.extent(0));
441  scalar_type rnorm = KokkosBlas::nrm2(r_1d);
442  std::cout << rnorm << " -> " << rnorm/bnorms(j) << ", ";
443  }
444  std::cout << std::endl << std::endl;
445  }
446 
447  // put X for output
448  Util::put_1d_data_helper_kokkos_view<MVAdapter, host_mvector_t>::
449  do_put(Xptr, X_view, ldx, CONTIGUOUS_AND_ROOTED, rowIndexBase);
450 
451  return numIters;
452 }
453 
454 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
455 bool
456 SolverCore<ConcreteSolver,Matrix,Vector>::matrixShapeOK()
457 {
458 #ifdef HAVE_AMESOS2_TIMERS
459  Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
460 #endif
461 
462  return( static_cast<solver_type*>(this)->matrixShapeOK_impl() );
463 }
464 
465  // The RCP should probably be to a const Matrix, but that throws a
466  // wrench in some of the traits mechanisms that aren't expecting
467  // const types
468 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
469 void
470 SolverCore<ConcreteSolver,Matrix,Vector>::setA( const Teuchos::RCP<const Matrix> a,
471  EPhase keep_phase )
472 {
473  matrixA_ = createConstMatrixAdapter(a);
474 
475 #ifdef HAVE_AMESOS2_DEBUG
476  TEUCHOS_TEST_FOR_EXCEPTION( (keep_phase != CLEAN) &&
477  (globalNumRows_ != matrixA_->getGlobalNumRows() ||
478  globalNumCols_ != matrixA_->getGlobalNumCols()),
479  std::invalid_argument,
480  "Dimensions of new matrix be the same as the old matrix if "
481  "keeping any solver phase" );
482 #endif
483 
484  status_.last_phase_ = keep_phase;
485 
486  // Reset phase counters
487  switch( status_.last_phase_ ){
488  case CLEAN:
489  status_.numPreOrder_ = 0;
490  // Intentional fallthrough.
491  case PREORDERING:
492  status_.numSymbolicFact_ = 0;
493  // Intentional fallthrough.
494  case SYMBFACT:
495  status_.numNumericFact_ = 0;
496  // Intentional fallthrough.
497  case NUMFACT: // probably won't ever happen by itself
498  status_.numSolve_ = 0;
499  // Intentional fallthrough.
500  case SOLVE: // probably won't ever happen
501  break;
502  }
503 
504  // Re-get the matrix dimensions in case they have changed
505  globalNumNonZeros_ = matrixA_->getGlobalNNZ();
506  globalNumCols_ = matrixA_->getGlobalNumCols();
507  globalNumRows_ = matrixA_->getGlobalNumRows();
508 }
509 
510 
511 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
512 Solver<Matrix,Vector>&
513 SolverCore<ConcreteSolver,Matrix,Vector>::setParameters(
514  const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
515 {
516 #ifdef HAVE_AMESOS2_TIMERS
517  Teuchos::TimeMonitor LocalTimer1(timers_.totalTime_);
518 #endif
519 
520  if( parameterList->name() == "Amesos2" ){
521  // First, validate the parameterList
522  Teuchos::RCP<const Teuchos::ParameterList> valid_params = getValidParameters();
523  parameterList->validateParameters(*valid_params);
524 
525  // Do everything here that is for generic status and control parameters
526  control_.setControlParameters(parameterList);
527 
528  // Finally, hook to the implementation's parameter list parser
529  // First check if there is a dedicated sublist for this solver and use that if there is
530  if( parameterList->isSublist(name()) ){
531  // Have control look through the solver's parameter list to see if
532  // there is anything it recognizes (mostly the "Transpose" parameter)
533  control_.setControlParameters(Teuchos::sublist(parameterList, name()));
534 
535  static_cast<solver_type*>(this)->setParameters_impl(Teuchos::sublist(parameterList, name()));
536  }
537  }
538 
539  return *this;
540 }
541 
542 
543 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
544 Teuchos::RCP<const Teuchos::ParameterList>
545 SolverCore<ConcreteSolver,Matrix,Vector>::getValidParameters() const
546 {
547 #ifdef HAVE_AMESOS2_TIMERS
548  Teuchos::TimeMonitor LocalTimer1( timers_.totalTime_ );
549 #endif
550 
551  using Teuchos::ParameterList;
552  using Teuchos::RCP;
553  using Teuchos::rcp;
554 
555  //RCP<ParameterList> control_params = rcp(new ParameterList("Amesos2 Control"));
556  RCP<ParameterList> control_params = rcp(new ParameterList("Amesos2"));
557  control_params->set("Transpose", false, "Whether to solve with the matrix transpose");
558  control_params->set("Iterative refinement", false, "Whether to solve with iterative refinement");
559  control_params->set("Number of iterative refinements", 2, "Number of iterative refinements");
560  control_params->set("Verboes for iterative refinement", false, "Verbosity for iterative refinements");
561  // control_params->set("AddToDiag", "");
562  // control_params->set("AddZeroToDiag", false);
563  // control_params->set("MatrixProperty", "general");
564  // control_params->set("Reindex", false);
565 
566  RCP<const ParameterList>
567  solver_params = static_cast<const solver_type*>(this)->getValidParameters_impl();
568  // inject the "Transpose" parameter into the solver's valid parameters
569  Teuchos::rcp_const_cast<ParameterList>(solver_params)->set("Transpose", false,
570  "Whether to solve with the "
571  "matrix transpose");
572 
573  RCP<ParameterList> amesos2_params = rcp(new ParameterList("Amesos2"));
574  amesos2_params->setParameters(*control_params);
575  amesos2_params->set(name(), *solver_params);
576 
577  return amesos2_params;
578 }
579 
580 
581 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
582 std::string
583 SolverCore<ConcreteSolver,Matrix,Vector>::description() const
584 {
585  std::ostringstream oss;
586  oss << name() << " solver interface";
587  return oss.str();
588 }
589 
590 
591 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
592 void
593 SolverCore<ConcreteSolver,Matrix,Vector>::describe(
594  Teuchos::FancyOStream &out,
595  const Teuchos::EVerbosityLevel verbLevel) const
596 {
597  if( matrixA_.is_null() || (rank_ != 0) ){ return; }
598  using std::endl;
599  using std::setw;
600  using Teuchos::VERB_DEFAULT;
601  using Teuchos::VERB_NONE;
602  using Teuchos::VERB_LOW;
603  using Teuchos::VERB_MEDIUM;
604  using Teuchos::VERB_HIGH;
605  using Teuchos::VERB_EXTREME;
606  Teuchos::EVerbosityLevel vl = verbLevel;
607  if (vl == VERB_DEFAULT) vl = VERB_LOW;
608  Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getComm();
609  size_t width = 1;
610  for( size_t dec = 10; dec < globalNumRows_; dec *= 10 ) {
611  ++width;
612  }
613  width = std::max<size_t>(width,size_t(11)) + 2;
614  Teuchos::OSTab tab(out);
615  // none: print nothing
616  // low: print O(1) info from node 0
617  // medium: print O(P) info, num entries per node
618  // high: print O(N) info, num entries per row
619  // extreme: print O(NNZ) info: print indices and values
620  //
621  // for medium and higher, print constituent objects at specified verbLevel
622  if( vl != VERB_NONE ) {
623  std::string p = name();
624  Util::printLine(out);
625  out << this->description() << std::endl << std::endl;
626 
627  out << p << "Matrix has " << globalNumRows_ << " rows"
628  << " and " << globalNumNonZeros_ << " nonzeros"
629  << std::endl;
630  if( vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME ){
631  out << p << "Nonzero elements per row = "
632  << globalNumNonZeros_ / globalNumRows_
633  << std::endl;
634  out << p << "Percentage of nonzero elements = "
635  << 100.0 * globalNumNonZeros_ / (globalNumRows_ * globalNumCols_)
636  << std::endl;
637  }
638  if( vl == VERB_HIGH || vl == VERB_EXTREME ){
639  out << p << "Use transpose = " << control_.useTranspose_
640  << std::endl;
641  out << p << "Use iterative refinement = " << control_.useIterRefine_
642  << std::endl;
643  }
644  if ( vl == VERB_EXTREME ){
645  printTiming(out,vl);
646  }
647  Util::printLine(out);
648  }
649 }
650 
651 
652 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
653 void
654 SolverCore<ConcreteSolver,Matrix,Vector>::printTiming(
655  Teuchos::FancyOStream &out,
656  const Teuchos::EVerbosityLevel verbLevel) const
657 {
658  if( matrixA_.is_null() || (rank_ != 0) ){ return; }
659 
660  double preTime = timers_.preOrderTime_.totalElapsedTime();
661  double symTime = timers_.symFactTime_.totalElapsedTime();
662  double numTime = timers_.numFactTime_.totalElapsedTime();
663  double solTime = timers_.solveTime_.totalElapsedTime();
664  double totTime = timers_.totalTime_.totalElapsedTime();
665  double overhead = totTime - (preTime + symTime + numTime + solTime);
666 
667  std::string p = name() + " : ";
668  Util::printLine(out);
669 
670  if(verbLevel != Teuchos::VERB_NONE)
671  {
672  out << p << "Time to convert matrix to implementation format = "
673  << timers_.mtxConvTime_.totalElapsedTime() << " (s)"
674  << std::endl;
675  out << p << "Time to redistribute matrix = "
676  << timers_.mtxRedistTime_.totalElapsedTime() << " (s)"
677  << std::endl;
678 
679  out << p << "Time to convert vectors to implementation format = "
680  << timers_.vecConvTime_.totalElapsedTime() << " (s)"
681  << std::endl;
682  out << p << "Time to redistribute vectors = "
683  << timers_.vecRedistTime_.totalElapsedTime() << " (s)"
684  << std::endl;
685 
686  out << p << "Number of pre-orderings = "
687  << status_.getNumPreOrder()
688  << std::endl;
689  out << p << "Time for pre-ordering = "
690  << preTime << " (s), avg = "
691  << preTime / status_.getNumPreOrder() << " (s)"
692  << std::endl;
693 
694  out << p << "Number of symbolic factorizations = "
695  << status_.getNumSymbolicFact()
696  << std::endl;
697  out << p << "Time for sym fact = "
698  << symTime << " (s), avg = "
699  << symTime / status_.getNumSymbolicFact() << " (s)"
700  << std::endl;
701 
702  out << p << "Number of numeric factorizations = "
703  << status_.getNumNumericFact()
704  << std::endl;
705  out << p << "Time for num fact = "
706  << numTime << " (s), avg = "
707  << numTime / status_.getNumNumericFact() << " (s)"
708  << std::endl;
709 
710  out << p << "Number of solve phases = "
711  << status_.getNumSolve()
712  << std::endl;
713  out << p << "Time for solve = "
714  << solTime << " (s), avg = "
715  << solTime / status_.getNumSolve() << " (s)"
716  << std::endl;
717 
718  out << p << "Total time spent in Amesos2 = "
719  << totTime << " (s)"
720  << std::endl;
721  out << p << "Total time spent in the Amesos2 interface = "
722  << overhead << " (s)"
723  << std::endl;
724  out << p << " (the above time does not include solver time)"
725  << std::endl;
726  out << p << "Amesos2 interface time / total time = "
727  << overhead / totTime
728  << std::endl;
729  Util::printLine(out);
730  }
731 }
732 
733 
734 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
735 void
736 SolverCore<ConcreteSolver,Matrix,Vector>::getTiming(
737  Teuchos::ParameterList& timingParameterList) const
738 {
739  Teuchos::ParameterList temp;
740  timingParameterList = temp.setName("NULL");
741 }
742 
743 
744 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
745 std::string
746 SolverCore<ConcreteSolver,Matrix,Vector>::name() const
747 {
748  std::string solverName = solver_type::name;
749  return solverName;
750 }
751 
752 template <template <class,class> class ConcreteSolver, class Matrix, class Vector >
753 void
754 SolverCore<ConcreteSolver,Matrix,Vector>::loadA(EPhase current_phase)
755 {
756  matrix_loaded_ = static_cast<solver_type*>(this)->loadA_impl(current_phase);
757 }
758 
759 
760 } // end namespace Amesos2
761 
762 #endif // AMESOS2_SOLVERCORE_DEF_HPP
Amesos2::SolverCore::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Return a const parameter list of all of the valid parameters that this-&gt;setParameterList(...) will accept.
Definition: Amesos2_SolverCore_def.hpp:545
Amesos2::SolverCore::solve
void solve() override
Solves (or )
Definition: Amesos2_SolverCore_def.hpp:151
Amesos2::SolverCore
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:71
Amesos2::SolverCore::~SolverCore
~SolverCore()
Destructor.
Definition: Amesos2_SolverCore_def.hpp:65
size
const int size
Definition: klu2_simple.cpp:50
Amesos2::EPhase
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:31
Amesos2::SolverCore::setParameters
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList) override
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:513
Amesos2::SolverCore::description
std::string description() const override
Returns a short description of this Solver.
Definition: Amesos2_SolverCore_def.hpp:583
Amesos2::SolverCore::setA
void setA(const Teuchos::RCP< const Matrix > a, EPhase keep_phase=CLEAN) override
Sets the matrix A of this solver.
Definition: Amesos2_SolverCore_def.hpp:470
Amesos2_Util.hpp
Utility functions for Amesos2.
Amesos2::SolverCore::preOrdering
super_type & preOrdering() override
Pre-orders the matrix A for minimal fill-in.
Definition: Amesos2_SolverCore_def.hpp:73
Amesos2::SolverCore::describe
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const override
Definition: Amesos2_SolverCore_def.hpp:593
Amesos2::SolverCore::printTiming
void printTiming(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Prints timing information about the current solver.
Definition: Amesos2_SolverCore_def.hpp:654
Amesos2::SolverCore::SolverCore
SolverCore(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize a Solver instance.
Definition: Amesos2_SolverCore_def.hpp:40
Amesos2::Solver::numericFactorization
virtual type & numericFactorization(void)=0
Performs numeric factorization on the matrix.
Amesos2::SolverCore::matrixShapeOK
bool matrixShapeOK() override
Returns true if the solver can handle this matrix shape.
Definition: Amesos2_SolverCore_def.hpp:456
Amesos2::Util::printLine
void printLine(Teuchos::FancyOStream &out)
Prints a line of 70 &quot;-&quot;s on std::cout.
Definition: Amesos2_Util.cpp:85
Amesos2::SolverCore::symbolicFactorization
super_type & symbolicFactorization() override
Performs symbolic factorization on the matrix A.
Definition: Amesos2_SolverCore_def.hpp:93
Amesos2::SolverCore::loadA
void loadA(EPhase current_phase)
Refresh this solver&#39;s internal data about A.
Definition: Amesos2_SolverCore_def.hpp:754
Amesos2::SolverCore::name
std::string name() const override
Return the name of this solver.
Definition: Amesos2_SolverCore_def.hpp:746
Amesos2::Solver
Interface to Amesos2 solver objects.
Definition: Amesos2_Solver_decl.hpp:44
Amesos2::SolverCore::getTiming
void getTiming(Teuchos::ParameterList &timingParameterList) const override
Extracts timing information from the current solver.
Definition: Amesos2_SolverCore_def.hpp:736
Amesos2::SolverCore::numericFactorization
super_type & numericFactorization() override
Performs numeric factorization on the matrix A.
Definition: Amesos2_SolverCore_def.hpp:122
Generated on Mon Oct 13 2025 07:10:06 for Amesos2 - Direct Sparse Solver Interfaces by   doxygen 1.8.5