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