HermesCommon  3.0
linear_matrix_solver.cpp
1 // This file is part of HermesCommon
2 //
3 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
4 // Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
5 //
6 // Hermes2D is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published
8 // by the Free Software Foundation; either version 2 of the License,
9 // or (at your option) any later version.
10 //
11 // Hermes2D is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Hermes2D; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #include "linear_matrix_solver.h"
30 #include "api.h"
31 #include "exceptions.h"
32 #include "util/memory_handling.h"
33 
34 using namespace Hermes::Algebra;
35 
36 namespace Hermes
37 {
38  namespace Solvers
39  {
40  template<typename Scalar>
41  LinearMatrixSolver<Scalar>::LinearMatrixSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs) : reuse_scheme(HERMES_CREATE_STRUCTURE_FROM_SCRATCH), general_matrix(matrix), general_rhs(rhs)
42  {
43  sln = nullptr;
44  time = -1.0;
45  n_eq = 1;
46  node_wise_ordering = false;
47  }
48 
49  template<typename Scalar>
51  {
52  free_with_check(sln);
53  }
54 
55  template<typename Scalar>
57  {
58  return sln;
59  }
60 
61  template<typename Scalar>
63  {
64  return this->general_matrix;
65  }
66 
67  template<typename Scalar>
69  {
70  return this->general_rhs;
71  }
72 
73  template<typename Scalar>
75  {
76  return time;
77  }
78 
79  template<typename Scalar>
81  {
82  set_reuse_scheme(HERMES_CREATE_STRUCTURE_FROM_SCRATCH);
83  }
84 
85  template<typename Scalar>
87  {
88  return reuse_scheme;
89  }
90 
91  template<typename Scalar>
93  {
94  this->reuse_scheme = reuse_scheme;
95  }
96 
97  template<typename Scalar>
99  {
100  this->reuse_scheme = HERMES_CREATE_STRUCTURE_FROM_SCRATCH;
101  this->n_eq = num_pdes;
102  this->node_wise_ordering = true;
103  }
104 
105  template<typename Scalar>
107  {
108  this->reuse_scheme = HERMES_CREATE_STRUCTURE_FROM_SCRATCH;
109  this->node_wise_ordering = false;
110  }
111 
112  template<typename Scalar>
114  {
115  throw Exceptions::MethodNotOverridenException("LinearMatrixSolver<Scalar>::get_residual_norm");
116  return 0.;
117  }
118 
119  template<typename Scalar>
121  {
122  DirectSolver<Scalar>* solver = dynamic_cast<DirectSolver<Scalar>*>(const_cast<LinearMatrixSolver<Scalar>*>(this));
123  if (solver)
124  return solver;
125  else
126  {
127  throw Hermes::Exceptions::LinearMatrixSolverException("Can not cast to DirectSolver.");
128  }
129  }
130 
131  template<typename Scalar>
133  {
134  LoopSolver<Scalar>* solver = dynamic_cast<LoopSolver<Scalar>*>(const_cast<LinearMatrixSolver<Scalar>*>(this));
135  if (solver)
136  return solver;
137  else
138  {
139  throw Hermes::Exceptions::LinearMatrixSolverException("Can not cast to LoopSolver.");
140  }
141  }
142 
143  template<typename Scalar>
145  {
146  IterSolver<Scalar>* solver = dynamic_cast<IterSolver<Scalar>*>(const_cast<LinearMatrixSolver<Scalar>*>(this));
147  if (solver)
148  return solver;
149  else
150  {
151  throw Hermes::Exceptions::LinearMatrixSolverException("Can not cast to IterSolver.");
152  }
153  }
154 
155  template<typename Scalar>
157  {
158  AMGSolver<Scalar>* solver = dynamic_cast<AMGSolver<Scalar>*>(const_cast<LinearMatrixSolver<Scalar>*>(this));
159  if (solver)
160  return solver;
161  else
162  {
163  throw Hermes::Exceptions::LinearMatrixSolverException("Can not cast to AMGSolver.");
164  }
165  }
166 
167  template<typename Scalar>
168  ExternalSolver<Scalar>* static_create_external_solver(CSCMatrix<Scalar> *m, SimpleVector<Scalar> *rhs)
169  {
170  throw Exceptions::MethodNotOverridenException("ExternalSolver<Scalar>::create_external_solver");
171  return nullptr;
172  }
173 
174  template<>
175  ExternalSolver<double>::creation ExternalSolver<double>::create_external_solver = static_create_external_solver < double > ;
176  template<>
177  ExternalSolver<std::complex<double> >::creation ExternalSolver<std::complex<double> >::create_external_solver = static_create_external_solver < std::complex<double> > ;
178 
179  template<>
180  HERMES_API LinearMatrixSolver<double>* create_linear_solver(Matrix<double>* matrix, Vector<double>* rhs, bool use_direct_solver)
181  {
182  Vector<double>* rhs_dummy = nullptr;
183  switch (use_direct_solver ? Hermes::HermesCommonApi.get_integral_param_value(Hermes::directMatrixSolverType) : Hermes::HermesCommonApi.get_integral_param_value(Hermes::matrixSolverType))
184  {
185  case Hermes::SOLVER_EXTERNAL:
186  {
187  if (rhs != nullptr)
188  return ExternalSolver<double>::create_external_solver(static_cast<CSCMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs));
189  else
190  return ExternalSolver<double>::create_external_solver(static_cast<CSCMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs_dummy));
191  }
192  case Hermes::SOLVER_AZTECOO:
193  {
194  if (use_direct_solver)
195  throw Hermes::Exceptions::Exception("The iterative solver AztecOO selected as a direct solver.");
196 #if defined HAVE_AZTECOO && defined HAVE_EPETRA
197  if (rhs != nullptr) return new AztecOOSolver<double>(static_cast<EpetraMatrix<double>*>(matrix), static_cast<EpetraVector<double>*>(rhs));
198  else return new AztecOOSolver<double>(static_cast<EpetraMatrix<double>*>(matrix), static_cast<EpetraVector<double>*>(rhs_dummy));
199 #else
200  throw Hermes::Exceptions::Exception("AztecOO not installed.");
201 #endif
202  break;
203  }
204  case Hermes::SOLVER_AMESOS:
205  {
206 #if defined HAVE_AMESOS && defined HAVE_EPETRA
207  if (rhs != nullptr) return new AmesosSolver<double>("Amesos_Klu", static_cast<EpetraMatrix<double>*>(matrix), static_cast<EpetraVector<double>*>(rhs));
208  else return new AmesosSolver<double>("Amesos_Klu", static_cast<EpetraMatrix<double>*>(matrix), static_cast<EpetraVector<double>*>(rhs_dummy));
209 #else
210  throw Hermes::Exceptions::Exception("Amesos not installed.");
211 #endif
212  break;
213  }
214  case Hermes::SOLVER_MUMPS:
215  {
216 #ifdef WITH_MUMPS
217  if (rhs != nullptr) return new MumpsSolver<double>(static_cast<MumpsMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs));
218  else return new MumpsSolver<double>(static_cast<MumpsMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs_dummy));
219 #else
220  throw Hermes::Exceptions::Exception("MUMPS was not installed.");
221 #endif
222  break;
223  }
224  case Hermes::SOLVER_PETSC:
225  {
226  if (use_direct_solver)
227  throw Hermes::Exceptions::Exception("The iterative solver PETSc selected as a direct solver.");
228 #ifdef WITH_PETSC
229  if (rhs != nullptr) return new PetscLinearMatrixSolver<double>(static_cast<PetscMatrix<double>*>(matrix), static_cast<PetscVector<double>*>(rhs));
230  else return new PetscLinearMatrixSolver<double>(static_cast<PetscMatrix<double>*>(matrix), static_cast<PetscVector<double>*>(rhs_dummy));
231 #else
232  throw Hermes::Exceptions::Exception("PETSc not installed.");
233 #endif
234  break;
235  }
236  case Hermes::SOLVER_UMFPACK:
237  {
238 #ifdef WITH_UMFPACK
239  if (rhs != nullptr) return new UMFPackLinearMatrixSolver<double>(static_cast<CSCMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs));
240  else return new UMFPackLinearMatrixSolver<double>(static_cast<CSCMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs_dummy));
241 #else
242  throw Hermes::Exceptions::Exception("UMFPACK was not installed.");
243 #endif
244  break;
245  }
246  case Hermes::SOLVER_PARALUTION_ITERATIVE:
247  {
248  if (use_direct_solver)
249  throw Hermes::Exceptions::Exception("The iterative solver PARALUTION selected as a direct solver.");
250 #ifdef WITH_PARALUTION
251  return new IterativeParalutionLinearMatrixSolver<double>(static_cast<ParalutionMatrix<double>*>(matrix), static_cast<ParalutionVector<double>*>(rhs));
252 #else
253  throw Hermes::Exceptions::Exception("PARALUTION was not installed.");
254 #endif
255  break;
256  }
257  case Hermes::SOLVER_PARALUTION_AMG:
258  {
259  if (use_direct_solver)
260  throw Hermes::Exceptions::Exception("The AMG solver PARALUTION selected as a direct solver.");
261 #ifdef WITH_PARALUTION
262  return new AMGParalutionLinearMatrixSolver<double>(static_cast<ParalutionMatrix<double>*>(matrix), static_cast<ParalutionVector<double>*>(rhs));
263 #else
264  throw Hermes::Exceptions::Exception("PARALUTION was not installed.");
265 #endif
266  break;
267  }
268  case Hermes::SOLVER_SUPERLU:
269  {
270 #ifdef WITH_SUPERLU
271  if (rhs != nullptr) return new SuperLUSolver<double>(static_cast<CSCMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs));
272  else return new SuperLUSolver<double>(static_cast<CSCMatrix<double>*>(matrix), static_cast<SimpleVector<double>*>(rhs_dummy));
273 #else
274  throw Hermes::Exceptions::Exception("SuperLU was not installed.");
275 #endif
276  break;
277  }
278  default:
279  throw Hermes::Exceptions::Exception("Unknown matrix solver requested in create_linear_solver().");
280  }
281  return nullptr;
282  }
283 
284  template<>
285  HERMES_API LinearMatrixSolver<std::complex<double> >* create_linear_solver(Matrix<std::complex<double> >* matrix, Vector<std::complex<double> >* rhs, bool use_direct_solver)
286  {
287  Vector<std::complex<double> >* rhs_dummy = nullptr;
288  switch (use_direct_solver ? Hermes::HermesCommonApi.get_integral_param_value(Hermes::directMatrixSolverType) : Hermes::HermesCommonApi.get_integral_param_value(Hermes::matrixSolverType))
289  {
290  case Hermes::SOLVER_EXTERNAL:
291  {
292  if (rhs != nullptr)
293  return ExternalSolver<std::complex<double> >::create_external_solver(static_cast<CSCMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs));
294  else
295  return ExternalSolver<std::complex<double> >::create_external_solver(static_cast<CSCMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs_dummy));
296  }
297  case Hermes::SOLVER_AZTECOO:
298  {
299  if (use_direct_solver)
300  throw Hermes::Exceptions::Exception("The iterative solver AztecOO selected as a direct solver.");
301 #if defined HAVE_AZTECOO && defined HAVE_EPETRA
302  if (rhs != nullptr) return new AztecOOSolver<std::complex<double> >(static_cast<EpetraMatrix<std::complex<double> >*>(matrix), static_cast<EpetraVector<std::complex<double> >*>(rhs));
303  else return new AztecOOSolver<std::complex<double> >(static_cast<EpetraMatrix<std::complex<double> >*>(matrix), static_cast<EpetraVector<std::complex<double> >*>(rhs_dummy));
304 #else
305  throw Hermes::Exceptions::Exception("AztecOO not installed.");
306 #endif
307  break;
308  }
309  case Hermes::SOLVER_AMESOS:
310  {
311 #if defined HAVE_AMESOS && defined HAVE_EPETRA
312  if (rhs != nullptr) return new AmesosSolver<std::complex<double> >("Amesos_Klu", static_cast<EpetraMatrix<std::complex<double> >*>(matrix), static_cast<EpetraVector<std::complex<double> >*>(rhs));
313  else return new AmesosSolver<std::complex<double> >("Amesos_Klu", static_cast<EpetraMatrix<std::complex<double> >*>(matrix), static_cast<EpetraVector<std::complex<double> >*>(rhs_dummy));
314 #else
315  throw Hermes::Exceptions::Exception("Amesos not installed.");
316 #endif
317  break;
318  }
319  case Hermes::SOLVER_MUMPS:
320  {
321 #ifdef WITH_MUMPS
322  if (rhs != nullptr) return new MumpsSolver<std::complex<double> >(static_cast<MumpsMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs));
323  else return new MumpsSolver<std::complex<double> >(static_cast<MumpsMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs_dummy));
324 #else
325  throw Hermes::Exceptions::Exception("MUMPS was not installed.");
326 #endif
327  break;
328  }
329  case Hermes::SOLVER_PETSC:
330  {
331  if (use_direct_solver)
332  throw Hermes::Exceptions::Exception("The iterative solver PETSc selected as a direct solver.");
333 #ifdef WITH_PETSC
334  if (rhs != nullptr) return new PetscLinearMatrixSolver<std::complex<double> >(static_cast<PetscMatrix<std::complex<double> >*>(matrix), static_cast<PetscVector<std::complex<double> >*>(rhs));
335  else return new PetscLinearMatrixSolver<std::complex<double> >(static_cast<PetscMatrix<std::complex<double> >*>(matrix), static_cast<PetscVector<std::complex<double> >*>(rhs_dummy));
336 #else
337  throw Hermes::Exceptions::Exception("PETSc not installed.");
338 #endif
339  break;
340  }
341  case Hermes::SOLVER_UMFPACK:
342  {
343 #ifdef WITH_UMFPACK
344  if (rhs != nullptr) return new UMFPackLinearMatrixSolver<std::complex<double> >(static_cast<CSCMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs));
345  else return new UMFPackLinearMatrixSolver<std::complex<double> >(static_cast<CSCMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs_dummy));
346 #else
347  throw Hermes::Exceptions::Exception("UMFPACK was not installed.");
348 #endif
349  break;
350  }
351  case Hermes::SOLVER_PARALUTION_ITERATIVE:
352  case Hermes::SOLVER_PARALUTION_AMG:
353  {
354  if (use_direct_solver)
355  throw Hermes::Exceptions::Exception("The iterative solver PARALUTION selected as a direct solver.");
356 #ifdef WITH_PARALUTION
357  throw Hermes::Exceptions::Exception("PARALUTION only works for real problems.");
358 #else
359  throw Hermes::Exceptions::Exception("PARALUTION was not installed.");
360 #endif
361  break;
362  }
363  case Hermes::SOLVER_SUPERLU:
364  {
365 #ifdef WITH_SUPERLU
366  if (rhs != nullptr) return new SuperLUSolver<std::complex<double> >(static_cast<CSCMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs));
367  else return new SuperLUSolver<std::complex<double> >(static_cast<CSCMatrix<std::complex<double> >*>(matrix), static_cast<SimpleVector<std::complex<double> >*>(rhs_dummy));
368 #else
369  throw Hermes::Exceptions::Exception("SuperLU was not installed.");
370 #endif
371  break;
372  }
373  default:
374  throw Hermes::Exceptions::Exception("Unknown matrix solver requested in create_linear_solver().");
375  }
376  return nullptr;
377  }
378 
379  template <typename Scalar>
381  {
382  }
383 
384  template <typename Scalar>
386  {
387  }
388 
389  template <typename Scalar>
391  {
392  solve(nullptr);
393  }
394 
395  template <typename Scalar>
396  void SimpleExternalSolver<Scalar>::solve(Scalar* initial_guess)
397  {
398  // Output.
399  this->process_matrix_output(this->m);
400  this->process_vector_output(this->rhs);
401 
402  // External process.
403  std::string resultFileName = this->command();
404 
405  // Handling of the result file.
406  this->sln = malloc_with_check<Scalar>(this->m->get_size());
408  temp.alloc(this->m->get_size());
409  temp.import_from_file((char*)resultFileName.c_str(), "x", EXPORT_FORMAT_PLAIN_ASCII);
410  memcpy(this->sln, temp.v, this->m->get_size() * sizeof(Scalar));
411  }
412 
413  template <typename Scalar>
415  {
416  }
417 
418  template <typename Scalar>
419  void DirectSolver<Scalar>::solve(Scalar* initial_guess)
420  {
421  this->solve();
422  }
423 
424  template <typename Scalar>
425  LoopSolver<Scalar>::LoopSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs) : LinearMatrixSolver<Scalar>(matrix, rhs), max_iters(10000), tolerance(1e-8)
426  {
427  }
428 
429  template<typename Scalar>
431  {
432  this->tolerance = tol;
433  this->toleranceType = AbsoluteTolerance;
434  }
435 
436  template<typename Scalar>
438  {
439  this->tolerance = tol;
440  this->toleranceType = toleranceType;
441  }
442 
443  template<typename Scalar>
445  {
446  this->max_iters = iters;
447  }
448 
449  template <typename Scalar>
450  IterSolver<Scalar>::IterSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs) : LoopSolver<Scalar>(matrix, rhs), precond_yes(false), iterSolverType(CG)
451  {
452  }
453 
454  template<typename Scalar>
456  {
457  this->iterSolverType = iterSolverType;
458  }
459 
460  template <typename Scalar>
461  AMGSolver<Scalar>::AMGSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs) : LoopSolver<Scalar>(matrix, rhs), smootherSolverType(CG), smootherPreconditionerType(MultiColoredSGS)
462  {
463  }
464 
465  template<typename Scalar>
467  {
468  this->smootherPreconditionerType = preconditionerType_;
469  this->smootherSolverType = solverType_;
470  }
471 
472  template class HERMES_API LinearMatrixSolver < double > ;
473  template class HERMES_API LinearMatrixSolver < std::complex<double> > ;
474  template class HERMES_API SimpleExternalSolver < double > ;
475  template class HERMES_API SimpleExternalSolver < std::complex<double> > ;
476  template class HERMES_API ExternalSolver < double > ;
477  template class HERMES_API ExternalSolver < std::complex<double> > ;
478  template class HERMES_API DirectSolver < double > ;
479  template class HERMES_API DirectSolver < std::complex<double> > ;
480  template class HERMES_API LoopSolver < double > ;
481  template class HERMES_API LoopSolver < std::complex<double> > ;
482  template class HERMES_API IterSolver < double > ;
483  template class HERMES_API IterSolver < std::complex<double> > ;
484  template class HERMES_API AMGSolver < double > ;
485  template class HERMES_API AMGSolver < std::complex<double> > ;
486  }
487 }
General (abstract) matrix representation in Hermes.
Definition: matrix.h:36
AmesosSolver class as an interface to Amesos.
General namespace for the Hermes library.
Linear matrix solver functionality.
Exception interface Basically a std::exception, but with a constructor with string and with print_msg...
Definition: exceptions.h:49
Main Hermes API.
Encapsulation of AztecOO linear solver.
General (abstract) vector representation in Hermes.
PARALUTION solver interface.
Encapsulation of PARALUTION AMG linear solver.
virtual void alloc(unsigned int ndofs)
Definition: vector.cpp:417
Namespace containing classes for vector / matrix operations.
File containing common definitions, and basic global enums etc. for HermesCommon. ...
PreconditionerType
The preconditioner type.
Definition: precond.h:42
General (abstract) sparse matrix representation in Hermes.
Abstract class for defining interface for Algebraic Multigrid solvers. Internal, though utilizable fo...
An example class for using external solvers that run a command and store the result in a file...
Plain ascii file lines contains row column and value.
Vector used with MUMPS solver.
Definition: vector.h:116
Special-purpose abstract class for using external solvers. For examples implementation, see the class SimpleExternalSolver.
Encapsulation of Amesos linear solver.
Definition: amesos_solver.h:39
AztecOOSolver class as an interface to AztecOO.
General CSC Matrix class. (can be used in umfpack, in that case use the CSCMatrix subclass...
Definition: cs_matrix.h:131
HERMES_COMMON_API Hermes::Api HermesCommonApi
Global instance used inside Hermes which is also accessible to users.
Definition: api.cpp:158
virtual void set_tolerance(double tol)
virtual void set_max_iters(int iters)
Abstract middle-class for solvers that work in a loop of a kind (iterative, multigrid, ...)
HERMES_API LinearMatrixSolver< Scalar > * create_linear_solver(Matrix< Scalar > *matrix, Vector< Scalar > *rhs, bool use_direct_solver=false)
Function returning a solver according to the users's choice.
General Paralution matrix.
Matrix used with MUMPS solver.
Definition: mumps_solver.h:77
File containing definition of exceptions classes.
UMFPACK solver interface.
Abstract class for defining solver interface.
MUMPS solver interface.
Abstract class for defining interface for iterative solvers. Internal, though utilizable for defining...
Encapsulation of UMFPACK linear solver.
Method is not overriden and should be.
Definition: exceptions.h:201
PETSc solver interface.
virtual void set_smoother(IterSolverType solverType, PreconditionerType preconditionerType)
Set smoother (an iterative linear matrix solver).
SuperLU solver interface.
Class representing the vector for UMFPACK.
Encapsulation of PARALUTION iterative linear solver.
void set_solver_type(IterSolverType iterSolverType)