iterative_linear_solver.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // This header defines a class for linear iterative solvers
27 
28 // Include guards
29 #ifndef OOMPH_ITERATIVE_LINEAR_SOLVER_HEADER
30 #define OOMPH_ITERATIVE_LINEAR_SOLVER_HEADER
31 
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35 #include <oomph-lib-config.h>
36 #endif
37 
38 
39 // oomph-lib headers
40 #include "matrices.h"
41 #include "problem.h"
42 #include "linear_solver.h"
43 #include "preconditioner.h"
44 
45 
46 namespace oomph
47 {
48  //===========================================================================
49  /// Base class for all linear iterative solvers.
50  /// This merely defines standard interfaces for linear iterative solvers,
51  /// so that different solvers can be used in a clean and transparent manner.
52  //===========================================================================
54  {
55  public:
56  /// Constructor: Set (default) trivial preconditioner and set
57  /// defaults for tolerance and max. number of iterations
59  {
60  // Set pointer to default preconditioner
62 
63  // Set default convergence tolerance
64  Tolerance = 1.0e-8;
65 
66  // Set maximum number of iterations
67  Max_iter = 100;
68 
69  // set default for document convergence history
71 
72  // set default
74 
76 
77  // By default the iterative solver is not used as preconditioner
79 
80  // Indicates whether this is the first time we call the solve
81  // method
83  }
84 
85  /// Broken copy constructor
87 
88  /// Broken assignment operator
89  void operator=(const IterativeLinearSolver&) = delete;
90 
91  /// Destructor (empty)
93 
94  /// Access function to preconditioner
96  {
97  return Preconditioner_pt;
98  }
99 
100  /// Access function to preconditioner (const version)
102  {
103  return Preconditioner_pt;
104  }
105 
106  /// Access to convergence tolerance
107  double& tolerance()
108  {
109  return Tolerance;
110  }
111 
112  /// Access to max. number of iterations
113  unsigned& max_iter()
114  {
115  return Max_iter;
116  }
117 
118  /// Number of iterations taken
119  virtual unsigned iterations() const = 0;
120 
121  /// Enable documentation of the convergence history
123  {
125  }
126 
127  /// Disable documentation of the convergence history
129  {
130  Doc_convergence_history = false;
131  }
132 
133  /// Write convergence history into file with specified filename
134  /// (automatically switches on doc). Optional second argument is a string
135  /// that can be used (as a zone title) to identify what case
136  /// we're running (e.g. what combination of linear solver and
137  /// preconditioner or parameter values are used).
139  const std::string& file_name, const std::string& zone_title = "")
140  {
141  // start docing
143 
144  // Close if it's open
145  if (Output_file_stream.is_open())
146  {
147  Output_file_stream.close();
148  }
149 
150  // Open new one
151  Output_file_stream.open(file_name.c_str());
152 
153  // Write tecplot zone header
154  Output_file_stream << "VARIABLES=\"iterations\", \"scaled residual\""
155  << std::endl;
156  Output_file_stream << "ZONE T=\"" << zone_title << "\"" << std::endl;
157  Output_file_stream << 0 << " " << 1.0 << std::endl;
158  }
159 
160  /// Close convergence history output stream
162  {
163  if (Output_file_stream.is_open()) Output_file_stream.close();
164  }
165 
166  /// returns the time taken to assemble the jacobian matrix and
167  /// residual vector
168  double jacobian_setup_time() const
169  {
170  return Jacobian_setup_time;
171  }
172 
173  /// return the time taken to solve the linear system
175  {
176  return Solution_time;
177  }
178 
179  /// returns the the time taken to setup the preconditioner
180  virtual double preconditioner_setup_time() const
181  {
183  }
184 
185  /// Setup the preconditioner before the solve
187  {
189  }
190 
191  /// Don't set up the preconditioner before the solve
193  {
195  }
196 
197  /// Throw an error if we don't converge within max_iter
199  {
201  }
202 
203  /// Don't throw an error if we don't converge within max_iter (default).
205  {
207  }
208 
209  /// Enables the iterative solver be used as preconditioner (when
210  /// calling the solve method it bypass the setup solver method --
211  /// currently only used by Trilinos solver ---)
213  {
215  }
216 
217  /// Disables the iterative solver be used as preconditioner (when
218  /// calling the solve method it bypass the setup solver method --
219  /// currently only used by Trilinos solver ---)
221  {
223  }
224 
225  protected:
226  /// Flag indicating if the convergence history is to be
227  /// documented
229 
230  /// Output file stream for convergence history
231  std::ofstream Output_file_stream;
232 
233  /// Default preconditioner: The base
234  /// class for preconditioners is a fully functional (if trivial!)
235  /// preconditioner.
237 
238  /// Convergence tolerance
239  double Tolerance;
240 
241  /// Maximum number of iterations
242  unsigned Max_iter;
243 
244  /// Pointer to the preconditioner
246 
247  /// Jacobian setup time
249 
250  /// linear solver solution time
252 
253  /// Preconditioner setup time
255 
256  /// indicates whether the preconditioner should be setup before
257  /// solve. Default = true;
259 
260  /// Should we throw an error instead of just returning when we hit
261  /// the max iterations?
263 
264  /// Use the iterative solver as preconditioner
266 
267  /// When the iterative solver is used a preconditioner then we call
268  /// the setup of solver method only once (the first time the solve
269  /// method is called)
271  };
272 
273 
274  /// ////////////////////////////////////////////////////////////////////////////
275  /// ////////////////////////////////////////////////////////////////////////////
276  /// ////////////////////////////////////////////////////////////////////////////
277 
278 
279  //======================================================================
280  /// The conjugate gradient method.
281  //======================================================================
282  template<typename MATRIX>
283  class CG : public IterativeLinearSolver
284  {
285  public:
286  /// Constructor
287  CG()
288  : Iterations(0),
289  Matrix_pt(0),
290  Resolving(false),
292  {
293  }
294 
295 
296  /// Destructor (cleanup storage)
297  virtual ~CG()
298  {
299  clean_up_memory();
300  }
301 
302  /// Broken copy constructor
303  CG(const CG&) = delete;
304 
305  /// Broken assignment operator
306  void operator=(const CG&) = delete;
307 
308  /// Overload disable resolve so that it cleans up memory too
310  {
312  clean_up_memory();
313  }
314 
315 
316  /// Solver: Takes pointer to problem and returns the results vector
317  /// which contains the solution of the linear system defined by
318  /// the problem's fully assembled Jacobian and residual vector.
319  void solve(Problem* const& problem_pt, DoubleVector& result);
320 
321  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
322  /// vector and returns the solution of the linear system.
323  void solve(DoubleMatrixBase* const& matrix_pt,
324  const DoubleVector& rhs,
325  DoubleVector& solution)
326  {
327  // Store the matrix if required
328  if ((Enable_resolve) && (!Resolving))
329  {
330  Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
331 
332  // Matrix has been passed in from the outside so we must not
333  // delete it
334  Matrix_can_be_deleted = false;
335  }
336 
337  // set the distribution
338  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt))
339  {
340  // the solver has the same distribution as the matrix if possible
341  this->build_distribution(
342  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
343  ->distribution_pt());
344  }
345  else
346  {
347  // the solver has the same distribution as the RHS
348  this->build_distribution(rhs.distribution_pt());
349  }
350 
351  // Call the helper function
352  this->solve_helper(matrix_pt, rhs, solution);
353  }
354 
355  /// Re-solve the system defined by the last assembled Jacobian
356  /// and the rhs vector specified here. Solution is returned in the
357  /// vector result.
358  void resolve(const DoubleVector& rhs, DoubleVector& result);
359 
360  /// Number of iterations taken
361  unsigned iterations() const
362  {
363  return Iterations;
364  }
365 
366 
367  private:
368  /// General interface to solve function
369  void solve_helper(DoubleMatrixBase* const& matrix_pt,
370  const DoubleVector& rhs,
371  DoubleVector& solution);
372 
373 
374  /// Cleanup data that's stored for resolve (if any has been stored)
376  {
377  if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
378  {
379  delete Matrix_pt;
380  Matrix_pt = 0;
381  }
382  }
383 
384  /// Number of iterations taken
385  unsigned Iterations;
386 
387  /// Pointer to matrix
388  MATRIX* Matrix_pt;
389 
390  /// Boolean flag to indicate if the solve is done in re-solve mode,
391  /// bypassing setup of matrix and preconditioner
392  bool Resolving;
393 
394  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
395  /// can be deleted.
397  };
398 
399 
400  /// ////////////////////////////////////////////////////////////////////////////
401  /// ////////////////////////////////////////////////////////////////////////////
402  /// ////////////////////////////////////////////////////////////////////////////
403 
404 
405  //======================================================================
406  /// The conjugate gradient method.
407  //======================================================================
408  template<typename MATRIX>
410  {
411  public:
412  /// Constructor
414  : Iterations(0),
415  Matrix_pt(0),
416  Resolving(false),
418  {
419  }
420 
421 
422  /// Destructor (cleanup storage)
423  virtual ~BiCGStab()
424  {
425  clean_up_memory();
426  }
427 
428  /// Broken copy constructor
429  BiCGStab(const BiCGStab&) = delete;
430 
431  /// Broken assignment operator
432  void operator=(const BiCGStab&) = delete;
433 
434  /// Overload disable resolve so that it cleans up memory too
436  {
438  clean_up_memory();
439  }
440 
441  /// Solver: Takes pointer to problem and returns the results vector
442  /// which contains the solution of the linear system defined by
443  /// the problem's fully assembled Jacobian and residual vector.
444  void solve(Problem* const& problem_pt, DoubleVector& result);
445 
446  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
447  /// vector and returns the solution of the linear system.
448  void solve(DoubleMatrixBase* const& matrix_pt,
449  const DoubleVector& rhs,
450  DoubleVector& solution)
451  {
452  // Store the matrix if required
453  if ((Enable_resolve) && (!Resolving))
454  {
455  Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
456 
457  // Matrix has been passed in from the outside so we must not
458  // delete it
459  Matrix_can_be_deleted = false;
460  }
461 
462  // set the distribution
463  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt))
464  {
465  // the solver has the same distribution as the matrix if possible
466  this->build_distribution(
467  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
468  ->distribution_pt());
469  }
470  else
471  {
472  // the solver has the same distribution as the RHS
473  this->build_distribution(rhs.distribution_pt());
474  }
475 
476  // Call the helper function
477  this->solve_helper(matrix_pt, rhs, solution);
478  }
479 
480  /// Linear-algebra-type solver: Takes pointer to a matrix
481  /// and rhs vector and returns the solution of the linear system
482  /// Call the broken base-class version. If you want this, please
483  /// implement it
484  void solve(DoubleMatrixBase* const& matrix_pt,
485  const Vector<double>& rhs,
486  Vector<double>& result)
487  {
488  LinearSolver::solve(matrix_pt, rhs, result);
489  }
490 
491 
492  /// Re-solve the system defined by the last assembled Jacobian
493  /// and the rhs vector specified here. Solution is returned in the
494  /// vector result.
495  void resolve(const DoubleVector& rhs, DoubleVector& result);
496 
497 
498  /// Number of iterations taken
499  unsigned iterations() const
500  {
501  return Iterations;
502  }
503 
504 
505  private:
506  /// General interface to solve function
507  void solve_helper(DoubleMatrixBase* const& matrix_pt,
508  const DoubleVector& rhs,
509  DoubleVector& solution);
510 
511  /// Cleanup data that's stored for resolve (if any has been stored)
513  {
514  if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
515  {
516  delete Matrix_pt;
517  Matrix_pt = 0;
518  }
519  }
520 
521  /// Number of iterations taken
522  unsigned Iterations;
523 
524  /// Pointer to matrix
525  MATRIX* Matrix_pt;
526 
527  /// Boolean flag to indicate if the solve is done in re-solve mode,
528  /// bypassing setup of matrix and preconditioner
529  bool Resolving;
530 
531  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
532  /// can be deleted.
534  };
535 
536 
537  /// ////////////////////////////////////////////////////////////////////
538  /// ////////////////////////////////////////////////////////////////////
539  /// ////////////////////////////////////////////////////////////////////
540 
541 
542  //====================================================================
543  /// Smoother class:
544  /// The smoother class is designed for to be used in conjunction with
545  /// multigrid. The action of the smoother should reduce the high
546  /// frequency errors. These methods are inefficient as stand-alone
547  /// solvers.
548  //====================================================================
550  {
551  public:
552  /// Empty constructor
553  Smoother() : Use_as_smoother(false) {}
554 
555  /// Virtual empty destructor
556  virtual ~Smoother(){};
557 
558  /// The smoother_solve function performs fixed number of iterations
559  /// on the system A*result=rhs. The number of (smoothing) iterations is
560  /// the same as the max. number of iterations in the underlying
561  /// IterativeLinearSolver class. Note that the result vector MUST NOT
562  /// re-initialised to zero (as it would typically be when the Smoother is
563  /// called as a iterative linear solver).
564  virtual void smoother_solve(const DoubleVector& rhs,
565  DoubleVector& result) = 0;
566 
567  /// Set up the smoother for the matrix specified by the pointer
568  virtual void smoother_setup(DoubleMatrixBase* matrix_pt) = 0;
569 
570  /// Self-test to check that all the dimensions of the inputs to
571  /// solve helper are consistent and everything that needs to be built, is.
572  template<typename MATRIX>
573  void check_validity_of_solve_helper_inputs(MATRIX* const& matrix_pt,
574  const DoubleVector& rhs,
575  DoubleVector& solution,
576  const double& n_dof);
577 
578  protected:
579  /// When a derived class object is being used as a smoother in
580  /// the MG solver (or elsewhere) the residual norm does not need to be
581  /// calculated because we're simply performing a fixed number of (smoothing)
582  /// iterations. This boolean is used as a flag to indicate that the
583  /// IterativeLinearSolver (which this class is by inheritance) is supposed
584  /// to act in this way.
586  }; // End of Smoother
587 
588 
589  /// ////////////////////////////////////////////////////////////////////////
590  /// ////////////////////////////////////////////////////////////////////////
591  /// ////////////////////////////////////////////////////////////////////////
592 
593 
594  //=========================================================================
595  /// The Gauss Seidel method
596  //=========================================================================
597  template<typename MATRIX>
598  class GS : public virtual Smoother
599  {
600  public:
601  /// Constructor
602  GS()
603  : Matrix_pt(0),
604  Iterations(0),
605  Resolving(false),
607  {
608  }
609 
610  /// Destructor (cleanup storage)
611  virtual ~GS()
612  {
613  clean_up_memory();
614  }
615 
616  /// Broken copy constructor
617  GS(const GS&) = delete;
618 
619  /// Broken assignment operator
620  void operator=(const GS&) = delete;
621 
622  /// Overload disable resolve so that it cleans up memory too
624  {
626  clean_up_memory();
627  } // End of disable_resolve
628 
629  /// Set up the smoother for the matrix specified by the pointer
631  {
632  // Assume the matrix has been passed in from the outside so we must
633  // not delete it. This is needed to avoid pre- and post-smoothers
634  // deleting the same matrix in the MG solver. If this was originally
635  // set to TRUE then this will be sorted out in the other functions
636  // from which this was called
637  Matrix_can_be_deleted = false;
638 
639  // Upcast the input matrix to system matrix to the type MATRIX
640  Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
641  } // End of smoother_setup
642 
643  /// The smoother_solve function performs fixed number of iterations
644  /// on the system A*result=rhs. The number of (smoothing) iterations is
645  /// the same as the max. number of iterations in the underlying
646  /// IterativeLinearSolver class.
647  void smoother_solve(const DoubleVector& rhs, DoubleVector& result)
648  {
649  // If you use a smoother but you don't want to calculate the residual
650  Use_as_smoother = true;
651 
652  // Call the helper function
653  solve_helper(Matrix_pt, rhs, result);
654  } // End of smoother_setup
655 
656  /// Solver: Takes pointer to problem and returns the results
657  /// vector which contains the solution of the linear system defined
658  /// by the problem's fully assembled Jacobian and residual vector.
659  void solve(Problem* const& problem_pt, DoubleVector& result);
660 
661  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
662  /// vector and returns the solution of the linear system.
663  void solve(DoubleMatrixBase* const& matrix_pt,
664  const DoubleVector& rhs,
665  DoubleVector& solution)
666  {
667  // Reset the Use_as_smoother_flag as the solver is not being used
668  // as a smoother
669  Use_as_smoother = false;
670 
671  // Set up the distribution
672  this->build_distribution(rhs.distribution_pt());
673 
674  // Store the matrix if required
675  if ((Enable_resolve) && (!Resolving))
676  {
677  // Upcast to the appropriate matrix type
678  Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
679  }
680 
681  // Matrix has been passed in from the outside so we must not delete it
682  Matrix_can_be_deleted = false;
683 
684  // Call the helper function
685  this->solve_helper(matrix_pt, rhs, solution);
686  } // End of solve
687 
688  /// Linear-algebra-type solver: Takes pointer to a matrix
689  /// and rhs vector and returns the solution of the linear system
690  /// Call the broken base-class version. If you want this, please
691  /// implement it
692  void solve(DoubleMatrixBase* const& matrix_pt,
693  const Vector<double>& rhs,
694  Vector<double>& result)
695  {
696  LinearSolver::solve(matrix_pt, rhs, result);
697  } // End of solve
698 
699  /// Re-solve the system defined by the last assembled Jacobian
700  /// and the rhs vector specified here. Solution is returned in the
701  /// vector result.
702  void resolve(const DoubleVector& rhs, DoubleVector& result)
703  {
704  // We are re-solving
705  Resolving = true;
706 
707 #ifdef PARANOID
708  if (Matrix_pt == 0)
709  {
710  throw OomphLibError("No matrix was stored -- cannot re-solve",
711  OOMPH_CURRENT_FUNCTION,
712  OOMPH_EXCEPTION_LOCATION);
713  }
714 #endif
715 
716  // Call linear algebra-style solver
717  solve(Matrix_pt, rhs, result);
718 
719  // Reset re-solving flag
720  Resolving = false;
721  } // End of resolve
722 
723  /// Returns the time taken to set up the preconditioner
725  {
726  throw OomphLibError(
727  "Gauss Seidel is not a preconditionable iterative solver",
728  OOMPH_CURRENT_FUNCTION,
729  OOMPH_EXCEPTION_LOCATION);
730  return 0;
731  } // End of preconditioner_setup_time
732 
733  /// Number of iterations taken
734  unsigned iterations() const
735  {
736  return Iterations;
737  } // End of iterations
738 
739  private:
740  /// General interface to solve function
741  void solve_helper(DoubleMatrixBase* const& matrix_pt,
742  const DoubleVector& rhs,
743  DoubleVector& solution);
744 
745  /// Cleanup data that's stored for resolve (if any has been stored)
747  {
748  // If the matrix pointer isn't null and we're allowed to delete it
749  // delete the matrix and assign the pointer the value NULL
750  if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
751  {
752  // Destroy the matrix
753  delete Matrix_pt;
754 
755  // Make it a null pointer
756  Matrix_pt = 0;
757  }
758  } // End of clean_up_memory
759 
760  /// System matrix pointer in the format specified by the template argument
761  MATRIX* Matrix_pt;
762 
763  /// Number of iterations taken
764  unsigned Iterations;
765 
766  /// Boolean flag to indicate if the solve is done in re-solve mode,
767  /// bypassing setup of matrix and preconditioner
768  bool Resolving;
769 
770  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
771  /// can be deleted.
773  };
774 
775  /// ////////////////////////////////////////////////////////////////////////
776  /// ////////////////////////////////////////////////////////////////////////
777  /// ////////////////////////////////////////////////////////////////////////
778 
779  //=========================================================================
780  /// Explicit template specialisation of the Gauss Seidel method for
781  /// compressed row format matrices
782  //=========================================================================
783  template<>
784  class GS<CRDoubleMatrix> : public virtual Smoother
785  {
786  public:
787  /// Constructor
788  GS()
789  : Matrix_pt(0),
790  Iterations(0),
791  Resolving(false),
793  {
794  }
795 
796  /// Destructor (cleanup storage)
797  virtual ~GS()
798  {
799  clean_up_memory();
800  }
801 
802  /// Broken copy constructor
803  GS(const GS&) = delete;
804 
805  /// Broken assignment operator
806  void operator=(const GS&) = delete;
807 
808  /// The smoother_solve function performs fixed number of iterations
809  /// on the system A*result=rhs. The number of (smoothing) iterations is
810  /// the same as the max. number of iterations in the underlying
811  /// IterativeLinearSolver class.
812  void smoother_solve(const DoubleVector& rhs, DoubleVector& result)
813  {
814  // If you use a smoother but you don't want to calculate the residual
815  Use_as_smoother = true;
816 
817  // Call the helper function
818  solve_helper(Matrix_pt, rhs, result);
819  } // End of smoother_solve
820 
821  /// Overload disable resolve so that it cleans up memory too
823  {
825  clean_up_memory();
826  } // End of disable_resolve
827 
828  /// Set up the smoother for the matrix specified by the pointer
830  {
831  // Assume the matrix has been passed in from the outside so we must
832  // not delete it. This is needed to avoid pre- and post-smoothers
833  // deleting the same matrix in the MG solver. If this was originally
834  // set to TRUE then this will be sorted out in the other functions
835  // from which this was called
836  Matrix_can_be_deleted = false;
837 
838  // Call the generic setup helper function
839  setup_helper(matrix_pt);
840  } // End of smoother_setup
841 
842  /// Generic setup function to sort out everything that needs to be
843  /// set up with regards to the input matrix
844  void setup_helper(DoubleMatrixBase* matrix_pt);
845 
846  /// Solver: Takes pointer to problem and returns the results vector
847  /// which contains the solution of the linear system defined by
848  /// the problem's fully assembled Jacobian and residual vector.
849  void solve(Problem* const& problem_pt, DoubleVector& result);
850 
851  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
852  /// vector and returns the solution of the linear system.
853  void solve(DoubleMatrixBase* const& matrix_pt,
854  const DoubleVector& rhs,
855  DoubleVector& solution)
856  {
857  // Reset the Use_as_smoother_flag as the solver is not being used
858  // as a smoother
859  Use_as_smoother = false;
860 
861  // Set up the distribution
862  this->build_distribution(rhs.distribution_pt());
863 
864  // Store the matrix if required
865  if ((Enable_resolve) && (!Resolving))
866  {
867  // Upcast to the appropriate matrix type and sort the matrix entries
868  // out so that the CRDoubleMatrix implementation of the Gauss-Seidel
869  // solver can be used
870  Matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
871  }
872  // We still need to sort the entries
873  else
874  {
875  // The system matrix here is a CRDoubleMatrix. To make use of the
876  // specific implementation of the solver for this type of matrix we
877  // need to make sure the entries are arranged correctly
878  dynamic_cast<CRDoubleMatrix*>(matrix_pt)->sort_entries();
879 
880  // Now get access to the vector Index_of_diagonal_entries
881  Index_of_diagonal_entries = dynamic_cast<CRDoubleMatrix*>(matrix_pt)
882  ->get_index_of_diagonal_entries();
883  }
884 
885  // Matrix has been passed in from the outside so we must not delete it
886  Matrix_can_be_deleted = false;
887 
888  // Call the helper function
889  solve_helper(matrix_pt, rhs, solution);
890  } // End of solve
891 
892  /// Linear-algebra-type solver: Takes pointer to a matrix
893  /// and rhs vector and returns the solution of the linear system
894  /// Call the broken base-class version. If you want this, please
895  /// implement it
896  void solve(DoubleMatrixBase* const& matrix_pt,
897  const Vector<double>& rhs,
898  Vector<double>& result)
899  {
900  LinearSolver::solve(matrix_pt, rhs, result);
901  } // End of solve
902 
903  /// Re-solve the system defined by the last assembled Jacobian
904  /// and the rhs vector specified here. Solution is returned in the
905  /// vector result.
906  void resolve(const DoubleVector& rhs, DoubleVector& result)
907  {
908  // We are re-solving
909  Resolving = true;
910 
911 #ifdef PARANOID
912  // If the matrix pointer is null
913  if (this->Matrix_pt == 0)
914  {
915  throw OomphLibError("No matrix was stored -- cannot re-solve",
916  OOMPH_CURRENT_FUNCTION,
917  OOMPH_EXCEPTION_LOCATION);
918  }
919 #endif
920 
921  // Call linear algebra-style solver
922  solve(Matrix_pt, rhs, result);
923 
924  // Reset re-solving flag
925  Resolving = false;
926  } // End of resolve
927 
928  /// Returns the time taken to set up the preconditioner
930  {
931  // Create the error message
932  std::string error_output_string = "Gauss Seidel is not a ";
933  error_output_string += "preconditionable iterative solver";
934 
935  // Throw an error
936  throw OomphLibError(
937  error_output_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
938 
939  // Return a value so the compiler doesn't throw up an error about
940  // no input being returned
941  return 0;
942  } // End of preconditioner_setup_time
943 
944  /// Number of iterations taken
945  unsigned iterations() const
946  {
947  // Return the number of iterations
948  return Iterations;
949  } // End of iterations
950 
951  private:
952  /// General interface to solve function
953  void solve_helper(DoubleMatrixBase* const& matrix_pt,
954  const DoubleVector& rhs,
955  DoubleVector& solution);
956 
957  /// Clean up data that's stored for resolve (if any has been stored)
959  {
960  // If the matrix pointer isn't null AND we're allowed to delete the
961  // matrix which is only when we create the matrix ourselves
962  if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
963  {
964  // Delete the matrix
965  delete Matrix_pt;
966 
967  // Assign the associated pointer the value NULL
968  Matrix_pt = 0;
969  }
970  } // End of clean_up_memory
971 
972  /// System matrix pointer in the format specified by the template argument
974 
975  /// Number of iterations taken
976  unsigned Iterations;
977 
978  /// Boolean flag to indicate if the solve is done in re-solve mode,
979  /// bypassing setup of matrix and preconditioner
980  bool Resolving;
981 
982  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
983  /// can be deleted.
985 
986  /// Vector whose i'th entry contains the index of the last entry
987  /// below or on the diagonal of the i'th row of the matrix
989  };
990 
991  /// ////////////////////////////////////////////////////////////////////
992  /// ////////////////////////////////////////////////////////////////////
993  /// ////////////////////////////////////////////////////////////////////
994 
995  //=========================================================================
996  /// Damped Jacobi "solver" templated by matrix type. The "solver"
997  /// exists in many different incarnations: It's an IterativeLinearSolver,
998  /// and a Smoother, all of which use the same basic iteration.
999  //=========================================================================
1000  template<typename MATRIX>
1001  class DampedJacobi : public virtual Smoother
1002  {
1003  public:
1004  /// Empty constructor
1005  DampedJacobi(const double& omega = 2.0 / 3.0) : Matrix_can_be_deleted(true)
1006  {
1007  // Damping factor
1008  Omega = omega;
1009  }
1010 
1011  /// Empty destructor
1013  {
1014  // Run the generic clean up function
1015  clean_up_memory();
1016  }
1017 
1018  /// Broken copy constructor
1019  DampedJacobi(const DampedJacobi&) = delete;
1020 
1021  /// Broken assignment operator
1022  void operator=(const DampedJacobi&) = delete;
1023 
1024  /// Cleanup data that's stored for resolve (if any has been stored)
1026  {
1027  // If the matrix pointer isn't null AND we're allowed to delete the
1028  // matrix which is only when we create the matrix ourselves
1029  if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
1030  {
1031  // Delete the matrix
1032  delete Matrix_pt;
1033 
1034  // Assign the associated pointer the value NULL
1035  Matrix_pt = 0;
1036  }
1037  } // End of clean_up_memory
1038 
1039  /// Setup: Pass pointer to the matrix and store in cast form
1041  {
1042  // Assume the matrix has been passed in from the outside so we must not
1043  // delete it
1044  Matrix_can_be_deleted = false;
1045 
1046  // Upcast to the appropriate matrix type
1047  Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
1048 
1049  // Extract the diagonal entries of the matrix and store them
1050  extract_diagonal_entries(matrix_pt);
1051  } // End of smoother_setup
1052 
1053  /// Function to extract the diagonal entries from the matrix
1055  {
1056  // If we're using a CRDoubleMatrix object
1057  if (dynamic_cast<CRDoubleMatrix*>(matrix_pt))
1058  {
1059  // The matrix diagonal (we need this when we need to calculate inv(D)
1060  // where D is the diagonal of A and it remains the same for all uses
1061  // of the iterative scheme so we can store it and call it in each
1062  // iteration)
1063  Matrix_diagonal =
1064  dynamic_cast<CRDoubleMatrix*>(Matrix_pt)->diagonal_entries();
1065  }
1066  // If we're using a complex matrix then diagonal entries has to be a
1067  // complex vector rather than a vector of doubles.
1068  else if (dynamic_cast<CCDoubleMatrix*>(matrix_pt))
1069  {
1070  // Make an ostringstream object to create an error message
1071  std::ostringstream error_message_stream;
1072 
1073  // Create the error message
1074  error_message_stream << "Damped Jacobi can only cater to real-valued "
1075  << "matrices. If you require a complex-valued "
1076  << "version, please write this yourself. "
1077  << "It is likely that the only difference will be "
1078  << "the use of complex vectors.";
1079 
1080  // Throw an error
1081  throw OomphLibError(error_message_stream.str(),
1082  OOMPH_CURRENT_FUNCTION,
1083  OOMPH_EXCEPTION_LOCATION);
1084  }
1085  // Just extract the diagonal entries normally
1086  else
1087  {
1088  // Calculate the number of rows in the matrix
1089  unsigned n_row = Matrix_pt->nrow();
1090 
1091  // Loop over the rows of the matrix
1092  for (unsigned i = 0; i < n_row; i++)
1093  {
1094  // Assign the i-th value of Matrix_diagonal
1095  Matrix_diagonal[i] = (*Matrix_pt)(i, i);
1096  }
1097  } // if (dynamic_cast<CRDoubleMatrix*>(matrix_pt))
1098 
1099  // Calculate the n.d.o.f.
1100  unsigned n_dof = Matrix_diagonal.size();
1101 
1102  // Find the reciprocal of the entries of Matrix_diagonal
1103  for (unsigned i = 0; i < n_dof; i++)
1104  {
1105  Matrix_diagonal[i] = 1.0 / Matrix_diagonal[i];
1106  }
1107  } // End of extract_diagonal_entries
1108 
1109  /// The smoother_solve function performs fixed number of iterations
1110  /// on the system A*result=rhs. The number of (smoothing) iterations is
1111  /// the same as the max. number of iterations in the underlying
1112  /// IterativeLinearSolver class.
1113  void smoother_solve(const DoubleVector& rhs, DoubleVector& solution)
1114  {
1115  // If you use a smoother but you don't want to calculate the residual
1116  Use_as_smoother = true;
1117 
1118  // Call the helper function
1119  solve_helper(Matrix_pt, rhs, solution);
1120  } // End of smoother_solve
1121 
1122  /// Use damped Jacobi iteration as an IterativeLinearSolver:
1123  /// This obtains the Jacobian matrix J and the residual vector r
1124  /// (needed for the Newton method) from the problem's get_jacobian
1125  /// function and returns the result of Jx=r.
1126  void solve(Problem* const& problem_pt, DoubleVector& result);
1127 
1128  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
1129  /// vector and returns the solution of the linear system.
1130  void solve(DoubleMatrixBase* const& matrix_pt,
1131  const DoubleVector& rhs,
1132  DoubleVector& solution)
1133  {
1134  // Matrix has been passed in from the outside so we must not delete it
1135  Matrix_can_be_deleted = false;
1136 
1137  // Indicate that the solver is not being used as a smoother
1138  Use_as_smoother = false;
1139 
1140  // Set up the distribution
1141  this->build_distribution(rhs.distribution_pt());
1142 
1143  // Store the matrix if required
1144  if ((Enable_resolve) && (!Resolving))
1145  {
1146  // Upcast to the appropriate matrix type
1147  Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
1148  }
1149 
1150  // Extract the diagonal entries of the matrix and store them
1151  extract_diagonal_entries(matrix_pt);
1152 
1153  // Call the helper function
1154  solve_helper(matrix_pt, rhs, solution);
1155  } // End of solve
1156 
1157  /// Re-solve the system defined by the last assembled Jacobian
1158  /// and the rhs vector specified here. Solution is returned in the
1159  /// vector result.
1160  void resolve(const DoubleVector& rhs, DoubleVector& result)
1161  {
1162  // We are re-solving
1163  Resolving = true;
1164 
1165 #ifdef PARANOID
1166  if (Matrix_pt == 0)
1167  {
1168  throw OomphLibError("No matrix was stored -- cannot re-solve",
1169  OOMPH_CURRENT_FUNCTION,
1170  OOMPH_EXCEPTION_LOCATION);
1171  }
1172 #endif
1173 
1174  // Call linear algebra-style solver
1175  solve(Matrix_pt, rhs, result);
1176 
1177  // Reset re-solving flag
1178  Resolving = false;
1179  } // End of resolve
1180 
1181  /// Number of iterations taken
1182  unsigned iterations() const
1183  {
1184  // Return the value of Iterations
1185  return Iterations;
1186  } // End of iterations
1187 
1188  private:
1189  /// This is where the actual work is done -- different
1190  /// implementations for different matrix types.
1191  void solve_helper(DoubleMatrixBase* const& matrix_pt,
1192  const DoubleVector& rhs,
1193  DoubleVector& solution);
1194 
1195  /// Pointer to the matrix
1196  MATRIX* Matrix_pt;
1197 
1198  /// Vector containing the diagonal entries of A
1200 
1201  /// Boolean flag to indicate if the solve is done in re-solve mode,
1202  /// bypassing setup of matrix and preconditioner
1204 
1205  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
1206  /// can be deleted.
1208 
1209  /// Number of iterations taken
1210  unsigned Iterations;
1211 
1212  /// Damping factor
1213  double Omega;
1214  };
1215 
1216 
1217  /// ////////////////////////////////////////////////////////////////////
1218  /// ////////////////////////////////////////////////////////////////////
1219  /// ////////////////////////////////////////////////////////////////////
1220 
1221 
1222  //======================================================================
1223  /// The GMRES method.
1224  //======================================================================
1225  template<typename MATRIX>
1227  {
1228  public:
1229  /// Constructor
1231  : Iterations(0),
1232  Matrix_pt(0),
1233  Resolving(false),
1234  Matrix_can_be_deleted(true),
1236  {
1237  Preconditioner_LHS = true;
1238  Iteration_restart = false;
1239  }
1240 
1241  /// Destructor (cleanup storage)
1242  virtual ~GMRES()
1243  {
1244  clean_up_memory();
1245  }
1246 
1247  /// Broken copy constructor
1248  GMRES(const GMRES&) = delete;
1249 
1250  /// Broken assignment operator
1251  void operator=(const GMRES&) = delete;
1252 
1253  /// Overload disable resolve so that it cleans up memory too
1255  {
1257  clean_up_memory();
1258  }
1259 
1260  /// function to enable the computation of the gradient
1262  {
1263  Compute_gradient = true;
1264  }
1265 
1266  /// Solver: Takes pointer to problem and returns the results vector
1267  /// which contains the solution of the linear system defined by
1268  /// the problem's fully assembled Jacobian and residual vector.
1269  void solve(Problem* const& problem_pt, DoubleVector& result);
1270 
1271  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
1272  /// vector and returns the solution of the linear system.
1273  void solve(DoubleMatrixBase* const& matrix_pt,
1274  const DoubleVector& rhs,
1275  DoubleVector& solution)
1276  {
1277  // setup the distribution
1278  this->build_distribution(rhs.distribution_pt());
1279 
1280  // Store the matrix if required
1281  if ((Enable_resolve) && (!Resolving))
1282  {
1283  Matrix_pt = dynamic_cast<MATRIX*>(matrix_pt);
1284 
1285  // Matrix has been passed in from the outside so we must not
1286  // delete it
1287  Matrix_can_be_deleted = false;
1288  }
1289 
1290  // Call the helper function
1291  this->solve_helper(matrix_pt, rhs, solution);
1292  }
1293 
1294 
1295  /// Linear-algebra-type solver: Takes pointer to a matrix
1296  /// and rhs vector and returns the solution of the linear system
1297  /// Call the broken base-class version. If you want this, please
1298  /// implement it
1299  void solve(DoubleMatrixBase* const& matrix_pt,
1300  const Vector<double>& rhs,
1301  Vector<double>& result)
1302  {
1303  LinearSolver::solve(matrix_pt, rhs, result);
1304  }
1305 
1306  /// Re-solve the system defined by the last assembled Jacobian
1307  /// and the rhs vector specified here. Solution is returned in the
1308  /// vector result.
1309  void resolve(const DoubleVector& rhs, DoubleVector& result);
1310 
1311  /// Number of iterations taken
1312  unsigned iterations() const
1313  {
1314  return Iterations;
1315  }
1316 
1317  /// access function indicating whether restarted GMRES is used
1318  bool iteration_restart() const
1319  {
1320  return Iteration_restart;
1321  }
1322 
1323  /// switches on iteration restarting and takes as an argument the
1324  /// number of iterations after which the construction of the
1325  /// orthogonalisation basis vectors should be restarted
1326  void enable_iteration_restart(const unsigned& restart)
1327  {
1328  Restart = restart;
1329  Iteration_restart = true;
1330  }
1331 
1332  /// switches off iteration restart
1334  {
1335  Iteration_restart = false;
1336  }
1337 
1338  /// Set left preconditioning (the default)
1340  {
1341  Preconditioner_LHS = true;
1342  }
1343 
1344  /// Enable right preconditioning
1346  {
1347  Preconditioner_LHS = false;
1348  }
1349 
1350  private:
1351  /// General interface to solve function
1352  void solve_helper(DoubleMatrixBase* const& matrix_pt,
1353  const DoubleVector& rhs,
1354  DoubleVector& solution);
1355 
1356  /// Cleanup data that's stored for resolve (if any has been stored)
1358  {
1359  if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
1360  {
1361  delete Matrix_pt;
1362  Matrix_pt = 0;
1363  }
1364  }
1365 
1366  /// Helper function to update the result vector using the result,
1367  /// x=x_0+V_m*y
1368  void update(const unsigned& k,
1369  const Vector<Vector<double>>& H,
1370  const Vector<double>& s,
1371  const Vector<DoubleVector>& v,
1372  DoubleVector& x)
1373  {
1374  // Make a local copy of s
1375  Vector<double> y(s);
1376 
1377  // Backsolve:
1378  for (int i = int(k); i >= 0; i--)
1379  {
1380  // Divide the i-th entry of y by the i-th diagonal entry of H
1381  y[i] /= H[i][i];
1382 
1383  // Loop over the previous values of y and update
1384  for (int j = i - 1; j >= 0; j--)
1385  {
1386  // Update the j-th entry of y
1387  y[j] -= H[i][j] * y[i];
1388  }
1389  } // for (int i=int(k);i>=0;i--)
1390 
1391  // Store the number of rows in the result vector
1392  unsigned n_x = x.nrow();
1393 
1394  // Build a temporary vector with entries initialised to 0.0
1395  DoubleVector temp(x.distribution_pt(), 0.0);
1396 
1397  // Build a temporary vector with entries initialised to 0.0
1398  DoubleVector z(x.distribution_pt(), 0.0);
1399 
1400  // Get access to the underlying values
1401  double* temp_pt = temp.values_pt();
1402 
1403  // Calculate x=Vy
1404  for (unsigned j = 0; j <= k; j++)
1405  {
1406  // Get access to j-th column of v
1407  const double* vj_pt = v[j].values_pt();
1408 
1409  // Loop over the entries of the vector, temp
1410  for (unsigned i = 0; i < n_x; i++)
1411  {
1412  temp_pt[i] += vj_pt[i] * y[j];
1413  }
1414  } // for (unsigned j=0;j<=k;j++)
1415 
1416  // If we're using LHS preconditioning
1417  if (Preconditioner_LHS)
1418  {
1419  // Since we're using LHS preconditioning the preconditioner is applied
1420  // to the matrix and RHS vector so we simply update the value of x
1421  x += temp;
1422  }
1423  // If we're using RHS preconditioning
1424  else
1425  {
1426  // Start the timer
1427  double t_start_prec = TimingHelpers::timer();
1428 
1429  // Since we're using RHS preconditioning the preconditioner is applied
1430  // to the solution vector
1432 
1433  // Calculate the time taken for the preconditioner solve
1435  (TimingHelpers::timer() - t_start_prec);
1436 
1437  // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
1438  // sparse linear systems", p.284]
1439  x += z;
1440  }
1441  } // End of update
1442 
1443  /// Helper function: Generate a plane rotation. This is done by
1444  /// finding the values of \f$ \cos(\theta) \f$ (i.e. cs) and \sin(\theta)
1445  /// (i.e. sn) such that:
1446  /// \f[ \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix} = \begin{bmatrix} r \newline 0 \end{bmatrix}, \f]
1447  /// where \f$ r=\sqrt{pow(dx,2)+pow(dy,2)} \f$. The values of a and b are
1448  /// given by:
1449  /// \f[ \cos\theta=\dfrac{dx}{\sqrt{pow(dx,2)+pow(dy,2)}}, \f]
1450  /// and
1451  /// \f[ \sin\theta=\dfrac{dy}{\sqrt{pow(dx,2)+pow(dy,2)}}. \f]
1452  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.192
1453  void generate_plane_rotation(double& dx, double& dy, double& cs, double& sn)
1454  {
1455  // If dy=0 then we do not need to apply a rotation
1456  if (dy == 0.0)
1457  {
1458  // Using theta=0 gives cos(theta)=1
1459  cs = 1.0;
1460 
1461  // Using theta=0 gives sin(theta)=0
1462  sn = 0.0;
1463  }
1464  // If dx or dy is large using the normal form of calculting cs and sn
1465  // is naive since this may overflow or underflow so instead we calculate
1466  // r=sqrt(pow(dx,2)+pow(dy,2)) by using r=|dy|sqrt(1+pow(dx/dy,2)) if
1467  // |dy|>|dx| [see <A
1468  // HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
1469  else if (fabs(dy) > fabs(dx))
1470  {
1471  // Since |dy|>|dx| calculate the ratio dx/dy
1472  double temp = dx / dy;
1473 
1474  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))
1475  sn = 1.0 / sqrt(1.0 + temp * temp);
1476 
1477  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))=(dx/dy)*sin(theta)
1478  cs = temp * sn;
1479  }
1480  // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
1481  // calculate the values of cs and sn using the method above
1482  else
1483  {
1484  // Since |dx|>=|dy| calculate the ratio dy/dx
1485  double temp = dy / dx;
1486 
1487  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))
1488  cs = 1.0 / sqrt(1.0 + temp * temp);
1489 
1490  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))=(dy/dx)*cos(theta)
1491  sn = temp * cs;
1492  }
1493  } // End of generate_plane_rotation
1494 
1495  /// Helper function: Apply plane rotation. This is done using the
1496  /// update:
1497  /// \f[ \begin{bmatrix} dx \newline dy \end{bmatrix} \leftarrow \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix}. \f]
1498  void apply_plane_rotation(double& dx, double& dy, double& cs, double& sn)
1499  {
1500  // Calculate the value of dx but don't update it yet
1501  double temp = cs * dx + sn * dy;
1502 
1503  // Set the value of dy
1504  dy = -sn * dx + cs * dy;
1505 
1506  // Set the value of dx using the correct values of dx and dy
1507  dx = temp;
1508  }
1509 
1510  /// Number of iterations taken
1511  unsigned Iterations;
1512 
1513  /// The number of iterations before the iteration proceedure is
1514  /// restarted if iteration restart is used
1515  unsigned Restart;
1516 
1517  /// boolean indicating if iteration restarting is used
1519 
1520  /// Pointer to matrix
1521  MATRIX* Matrix_pt;
1522 
1523  /// Boolean flag to indicate if the solve is done in re-solve mode,
1524  /// bypassing setup of matrix and preconditioner
1526 
1527  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
1528  /// can be deleted.
1530 
1531  /// boolean indicating use of left hand preconditioning (if true)
1532  /// or right hand preconditioning (if false)
1534 
1535  /// Storage for the time spent applying the preconditioner
1537  };
1538 
1539 
1540  /// ////////////////////////////////////////////////////////////////////
1541  /// ////////////////////////////////////////////////////////////////////
1542  /// ////////////////////////////////////////////////////////////////////
1543 
1544 
1545  //======================================================================
1546  /// The GMRES method.
1547  //======================================================================
1549  {
1550  public:
1551  /// Constructor
1553  DoubleVector* c_pt,
1554  double* x_pt,
1555  double* rhs_pt)
1556  : Iterations(0),
1557  Matrix_pt(0),
1558  B_pt(b_pt),
1559  C_pt(c_pt),
1560  X_pt(x_pt),
1561  Rhs_pt(rhs_pt),
1563  Resolving(false),
1564  Matrix_can_be_deleted(true)
1565  {
1566  Preconditioner_LHS = true;
1567  Iteration_restart = false;
1568  }
1569 
1570  /// Destructor: Clean up storage
1572  {
1573  clean_up_memory();
1574  }
1575 
1576  /// Broken copy constructor
1578 
1579  /// Broken assignment operator
1580  void operator=(const AugmentedProblemGMRES&) = delete;
1581 
1582  /// Overload disable resolve so that it cleans up memory too
1584  {
1586  clean_up_memory();
1587  } // End of disable_resolve
1588 
1589  /// Solver: Takes pointer to problem and returns the results vector
1590  /// which contains the solution of the linear system defined by
1591  /// the problem's fully assembled Jacobian and residual vector.
1592  void solve(Problem* const& problem_pt, DoubleVector& result);
1593 
1594  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
1595  /// vector and returns the solution of the linear system.
1596  void solve(DoubleMatrixBase* const& matrix_pt,
1597  const DoubleVector& rhs,
1598  DoubleVector& solution)
1599  {
1600  // Upcast to a CRDoubleMatrix
1601  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
1602 
1603  // If we can't upcast to a CRDoubleMatrix then this won't work...
1604  if (cr_matrix_pt == 0)
1605  {
1606  // Throw an error
1607  throw OomphLibError("Can't upcast input matrix to a CRDoubleMatrix!",
1608  OOMPH_CURRENT_FUNCTION,
1609  OOMPH_EXCEPTION_LOCATION);
1610  }
1611 
1612  // Set up the distribution
1613  this->build_distribution(cr_matrix_pt->distribution_pt());
1614 
1615  // Store the matrix if required
1616  if ((Enable_resolve) && (!Resolving))
1617  {
1618  // Store a pointer to the upcasted matrix
1619  Matrix_pt = cr_matrix_pt;
1620 
1621  // Matrix has been passed in from the outside so we must not
1622  // delete it
1623  Matrix_can_be_deleted = false;
1624  }
1625 
1626  // Call the helper function
1627  this->solve_helper(matrix_pt, rhs, solution);
1628  } // End of solve
1629 
1630  /// Linear-algebra-type solver: Takes pointer to a matrix
1631  /// and rhs vector and returns the solution of the linear system
1632  /// Call the broken base-class version. If you want this, please
1633  /// implement it
1634  void solve(DoubleMatrixBase* const& matrix_pt,
1635  const Vector<double>& rhs,
1636  Vector<double>& result)
1637  {
1638  LinearSolver::solve(matrix_pt, rhs, result);
1639  }
1640 
1641  /// Re-solve the system defined by the last assembled Jacobian
1642  /// and the rhs vector specified here. Solution is returned in the
1643  /// vector result.
1644  void resolve(const DoubleVector& rhs, DoubleVector& result);
1645 
1646  /// Number of iterations taken
1647  unsigned iterations() const
1648  {
1649  return Iterations;
1650  }
1651 
1652  /// access function indicating whether restarted GMRES is used
1653  bool iteration_restart() const
1654  {
1655  return Iteration_restart;
1656  }
1657 
1658  /// switches on iteration restarting and takes as an argument the
1659  /// number of iterations after which the construction of the
1660  /// orthogonalisation basis vectors should be restarted
1661  void enable_iteration_restart(const unsigned& restart)
1662  {
1663  Restart = restart;
1664  Iteration_restart = true;
1665  }
1666 
1667  /// Switches off iteration restart
1669  {
1670  Iteration_restart = false;
1671  } // End of disable_iteration_restart
1672 
1673  /// Set left preconditioning (the default)
1675  {
1676  Preconditioner_LHS = true;
1677  }
1678 
1679  /// Enable right preconditioning
1681  {
1682  Preconditioner_LHS = false;
1683  }
1684 
1685  private:
1686  /// Cleanup data that's stored for resolve (if any has been stored)
1688  {
1689  if ((Matrix_pt != 0) && (Matrix_can_be_deleted))
1690  {
1691  delete Matrix_pt;
1692  Matrix_pt = 0;
1693  }
1694  } // End of clean_up_memory
1695 
1696  /// Multiply the vector x by the augmented system matrix
1698  const DoubleVector& x,
1699  DoubleVector& soln)
1700  {
1701  // How many dofs are there in the non-augmented system?
1702  unsigned n_dof = matrix_pt->nrow();
1703 
1704  // Allocate space for the non-augmented component of x
1705  DoubleVector x_small(this->distribution_pt(), 0.0);
1706 
1707  // Loop over the first n_dof entries of x
1708  for (unsigned i = 0; i < n_dof; i++)
1709  {
1710  // Copy the i-th entry over
1711  x_small[i] = x[i];
1712  }
1713 
1714  // Allocate space for the product of the matrix and x_small
1715  DoubleVector a_prod_xsmall(this->distribution_pt(), 0.0);
1716 
1717  // Now multiply the matrix and x_small
1718  matrix_pt->multiply(x_small, a_prod_xsmall);
1719 
1720  // Get the scalar component of x associated with the system augmentation
1721  double y = x[n_dof];
1722 
1723  // Loop over the first n_dof entries of soln
1724  for (unsigned i = 0; i < n_dof; i++)
1725  {
1726  // Compute the i-th entry
1727  soln[i] = a_prod_xsmall[i] + y * (*B_pt)[i];
1728  }
1729 
1730  // Compute the final entry of soln
1731  soln[n_dof] = C_pt->dot(x_small);
1732  } // End of augmented_matrix_multiply
1733 
1734  /// Apply the block-diagonal Schur complement preconditioner to
1735  /// compute the LHS which has size N+1 (the augmented system size)
1737  DoubleVector& soln)
1738  {
1739  // How many dofs are there in the non-augmented system?
1740  unsigned n_dof = this->distribution_pt()->nrow();
1741 
1742  // Compute the final entry of r first
1743  soln[n_dof] = rhs[n_dof] / Schur_complement_scalar;
1744 
1745  // Do we want to use a block-diagonal approximation? The default is to
1746  // use a block upper-triangular approximation for the preconditioner
1747  bool use_block_diagonal_preconditioner = false;
1748 
1749  // Allocate space for (the non-augmented part of) r satisfying b-Jx = Mr
1750  DoubleVector r_small(this->distribution_pt(), 0.0);
1751 
1752  // Allocate space for the non-augmented system part of the RHS vector
1753  DoubleVector rhs_small(this->distribution_pt(), 0.0);
1754 
1755  // Loop over the first n_dof entries of the RHS vector
1756  for (unsigned i = 0; i < n_dof; i++)
1757  {
1758  // If we want to use the block-diagonal preconditioner
1759  if (use_block_diagonal_preconditioner)
1760  {
1761  // Copy the i-th entry over
1762  rhs_small[i] = rhs[i];
1763  }
1764  // If we want to use the block upper-triangular preconditioner
1765  else
1766  {
1767  // Copy the i-th entry over
1768  rhs_small[i] = rhs[i] - soln[n_dof] * (*B_pt)[i];
1769  }
1770  } // for (unsigned i=0;i<n_dof;i++)
1771 
1772  // Apply the preconditioner
1773  preconditioner_pt()->preconditioner_solve(rhs_small, r_small);
1774 
1775  // Loop over the first n_dof entries of the solution vector
1776  for (unsigned i = 0; i < n_dof; i++)
1777  {
1778  // Copy the i-th entry over
1779  soln[i] = r_small[i];
1780  }
1781  } // End of apply_schur_complement_preconditioner
1782 
1783  /// General interface to solve function
1784  void solve_helper(DoubleMatrixBase* const& matrix_pt,
1785  const DoubleVector& rhs,
1786  DoubleVector& solution);
1787 
1788  /// Helper function to update the result vector using the result,
1789  /// x=x_0+V_m*y
1790  void update(const unsigned& k,
1791  const Vector<Vector<double>>& H,
1792  const Vector<double>& s,
1793  const Vector<DoubleVector>& v,
1794  DoubleVector& x)
1795  {
1796  // Make a local copy of s
1797  Vector<double> y(s);
1798 
1799  // Backsolve:
1800  for (int i = int(k); i >= 0; i--)
1801  {
1802  // Divide the i-th entry of y by the i-th diagonal entry of H
1803  y[i] /= H[i][i];
1804 
1805  // Loop over the previous values of y and update
1806  for (int j = i - 1; j >= 0; j--)
1807  {
1808  // Update the j-th entry of y
1809  y[j] -= H[i][j] * y[i];
1810  }
1811  } // for (int i=int(k);i>=0;i--)
1812 
1813  // Store the number of rows in the result vector
1814  unsigned n_x = x.nrow();
1815 
1816  // Build a temporary vector with entries initialised to 0.0
1817  DoubleVector temp(x.distribution_pt(), 0.0);
1818 
1819  // Build a temporary vector with entries initialised to 0.0
1820  DoubleVector z(x.distribution_pt(), 0.0);
1821 
1822  // Get access to the underlying values
1823  double* temp_pt = temp.values_pt();
1824 
1825  // Calculate x=Vy
1826  for (unsigned j = 0; j <= k; j++)
1827  {
1828  // Get access to j-th column of v
1829  const double* vj_pt = v[j].values_pt();
1830 
1831  // Loop over the entries of the vector, temp
1832  for (unsigned i = 0; i < n_x; i++)
1833  {
1834  temp_pt[i] += vj_pt[i] * y[j];
1835  }
1836  } // for (unsigned j=0;j<=k;j++)
1837 
1838  // If we're using LHS preconditioning
1839  if (Preconditioner_LHS)
1840  {
1841  // Since we're using LHS preconditioning the preconditioner is applied
1842  // to the matrix and RHS vector so we simply update the value of x
1843  x += temp;
1844  }
1845  // If we're using RHS preconditioning
1846  else
1847  {
1848  // Since we're using RHS preconditioning the preconditioner is applied
1849  // to the solution vector
1851 
1852  // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
1853  // sparse linear systems", p.284]
1854  x += z;
1855  }
1856  } // End of update
1857 
1858  /// Helper function: Generate a plane rotation. This is done by
1859  /// finding the values of \f$ \cos(\theta) \f$ (i.e. cs) and \sin(\theta)
1860  /// (i.e. sn) such that:
1861  /// \f[ \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix} = \begin{bmatrix} r \newline 0 \end{bmatrix}, \f]
1862  /// where \f$ r=\sqrt{pow(dx,2)+pow(dy,2)} \f$. The values of a and b are
1863  /// given by:
1864  /// \f[ \cos\theta=\dfrac{dx}{\sqrt{pow(dx,2)+pow(dy,2)}}, \f]
1865  /// and
1866  /// \f[ \sin\theta=\dfrac{dy}{\sqrt{pow(dx,2)+pow(dy,2)}}. \f]
1867  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.192
1868  void generate_plane_rotation(double& dx, double& dy, double& cs, double& sn)
1869  {
1870  // If dy=0 then we do not need to apply a rotation
1871  if (dy == 0.0)
1872  {
1873  // Using theta=0 gives cos(theta)=1
1874  cs = 1.0;
1875 
1876  // Using theta=0 gives sin(theta)=0
1877  sn = 0.0;
1878  }
1879  // If dx or dy is large using the normal form of calculting cs and sn
1880  // is naive since this may overflow or underflow so instead we calculate
1881  // r=sqrt(pow(dx,2)+pow(dy,2)) by using r=|dy|sqrt(1+pow(dx/dy,2)) if
1882  // |dy|>|dx| [see <A
1883  // HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
1884  else if (fabs(dy) > fabs(dx))
1885  {
1886  // Since |dy|>|dx| calculate the ratio dx/dy
1887  double temp = dx / dy;
1888 
1889  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))
1890  sn = 1.0 / sqrt(1.0 + temp * temp);
1891 
1892  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))=(dx/dy)*sin(theta)
1893  cs = temp * sn;
1894  }
1895  // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
1896  // calculate the values of cs and sn using the method above
1897  else
1898  {
1899  // Since |dx|>=|dy| calculate the ratio dy/dx
1900  double temp = dy / dx;
1901 
1902  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))
1903  cs = 1.0 / sqrt(1.0 + temp * temp);
1904 
1905  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))=(dy/dx)*cos(theta)
1906  sn = temp * cs;
1907  }
1908  } // End of generate_plane_rotation
1909 
1910  /// Helper function: Apply plane rotation. This is done using the
1911  /// update:
1912  /// \f[ \begin{bmatrix} dx \newline dy \end{bmatrix} \leftarrow \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix}. \f]
1913  void apply_plane_rotation(double& dx, double& dy, double& cs, double& sn)
1914  {
1915  // Calculate the value of dx but don't update it yet
1916  double temp = cs * dx + sn * dy;
1917 
1918  // Set the value of dy
1919  dy = -sn * dx + cs * dy;
1920 
1921  // Set the value of dx using the correct values of dx and dy
1922  dx = temp;
1923  }
1924 
1925  /// Number of iterations taken
1926  unsigned Iterations;
1927 
1928  /// The number of iterations before the iteration proceedure is
1929  /// restarted if iteration restart is used
1930  unsigned Restart;
1931 
1932  /// boolean indicating if iteration restarting is used
1934 
1935  /// Pointer to matrix
1937 
1938  /// Pointer to the column vector in the bordered system
1940 
1941  /// Pointer to the row vector in the bordered system
1943 
1944  /// Pointer to the last entry of the LHS vector in the bordered system
1945  double* X_pt;
1946 
1947  /// Pointer to the last entry of the RHS vector in the bordered system
1948  double* Rhs_pt;
1949 
1950  /// The scalar component of the Schur complement preconditioner
1952 
1953  /// Boolean flag to indicate if the solve is done in re-solve mode,
1954  /// bypassing setup of matrix and preconditioner
1956 
1957  /// Boolean flag to indicate if the matrix pointed to be Matrix_pt
1958  /// can be deleted.
1960 
1961  /// boolean indicating use of left hand preconditioning (if true)
1962  /// or right hand preconditioning (if false)
1964  };
1965 } // End of namespace oomph
1966 
1967 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
AugmentedProblemGMRES(DoubleVector *b_pt, DoubleVector *c_pt, double *x_pt, double *rhs_pt)
Constructor.
unsigned iterations() const
Number of iterations taken.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update:
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
DoubleVector * C_pt
Pointer to the row vector in the bordered system.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void update(const unsigned &k, const Vector< Vector< double >> &H, const Vector< double > &s, const Vector< DoubleVector > &v, DoubleVector &x)
Helper function to update the result vector using the result, x=x_0+V_m*y.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
CRDoubleMatrix * Matrix_pt
Pointer to matrix.
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
bool Iteration_restart
boolean indicating if iteration restarting is used
unsigned Restart
The number of iterations before the iteration proceedure is restarted if iteration restart is used.
void disable_iteration_restart()
Switches off iteration restart.
AugmentedProblemGMRES(const AugmentedProblemGMRES &)=delete
Broken copy constructor.
bool iteration_restart() const
access function indicating whether restarted GMRES is used
double * Rhs_pt
Pointer to the last entry of the RHS vector in the bordered system.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void apply_schur_complement_preconditioner(const DoubleVector &rhs, DoubleVector &soln)
Apply the block-diagonal Schur complement preconditioner to compute the LHS which has size N+1 (the a...
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void generate_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Generate a plane rotation. This is done by finding the values of (i....
void enable_iteration_restart(const unsigned &restart)
switches on iteration restarting and takes as an argument the number of iterations after which the co...
double Schur_complement_scalar
The scalar component of the Schur complement preconditioner.
double * X_pt
Pointer to the last entry of the LHS vector in the bordered system.
void augmented_matrix_multiply(CRDoubleMatrix *matrix_pt, const DoubleVector &x, DoubleVector &soln)
Multiply the vector x by the augmented system matrix.
unsigned Iterations
Number of iterations taken.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void operator=(const AugmentedProblemGMRES &)=delete
Broken assignment operator.
DoubleVector * B_pt
Pointer to the column vector in the bordered system.
virtual ~AugmentedProblemGMRES()
Destructor: Clean up storage.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
void set_preconditioner_RHS()
Enable right preconditioning.
void set_preconditioner_LHS()
Set left preconditioning (the default)
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
unsigned iterations() const
Number of iterations taken.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
MATRIX * Matrix_pt
Pointer to matrix.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
unsigned Iterations
Number of iterations taken.
virtual ~BiCGStab()
Destructor (cleanup storage)
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void operator=(const BiCGStab &)=delete
Broken assignment operator.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
BiCGStab(const BiCGStab &)=delete
Broken copy constructor.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
Definition: matrices.h:2791
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
CG()
Constructor.
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
unsigned iterations() const
Number of iterations taken.
CG(const CG &)=delete
Broken copy constructor.
virtual ~CG()
Destructor (cleanup storage)
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void operator=(const CG &)=delete
Broken assignment operator.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
unsigned Iterations
Number of iterations taken.
MATRIX * Matrix_pt
Pointer to matrix.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1782
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
DampedJacobi(const DampedJacobi &)=delete
Broken copy constructor.
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
~DampedJacobi()
Empty destructor.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
unsigned Iterations
Number of iterations taken.
void operator=(const DampedJacobi &)=delete
Broken assignment operator.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
void extract_diagonal_entries(DoubleMatrixBase *matrix_pt)
Function to extract the diagonal entries from the matrix.
double Omega
Damping factor.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
This is where the actual work is done – different implementations for different matrix types.
void smoother_setup(DoubleMatrixBase *matrix_pt)
Setup: Pass pointer to the matrix and store in cast form.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
unsigned iterations() const
Number of iterations taken.
DampedJacobi(const double &omega=2.0/3.0)
Empty constructor.
Vector< double > Matrix_diagonal
Vector containing the diagonal entries of A.
void solve(Problem *const &problem_pt, DoubleVector &result)
Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the r...
MATRIX * Matrix_pt
Pointer to the matrix.
void smoother_solve(const DoubleVector &rhs, DoubleVector &solution)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned nrow() const
access function to the number of global rows.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
Definition: matrices.h:261
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
double * values_pt()
access function to the underlying values
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
void enable_computation_of_gradient()
function to enable the computation of the gradient
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
double Preconditioner_application_time
Storage for the time spent applying the preconditioner.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
MATRIX * Matrix_pt
Pointer to matrix.
void set_preconditioner_RHS()
Enable right preconditioning.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
GMRES(const GMRES &)=delete
Broken copy constructor.
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
unsigned Iterations
Number of iterations taken.
void operator=(const GMRES &)=delete
Broken assignment operator.
void update(const unsigned &k, const Vector< Vector< double >> &H, const Vector< double > &s, const Vector< DoubleVector > &v, DoubleVector &x)
Helper function to update the result vector using the result, x=x_0+V_m*y.
unsigned iterations() const
Number of iterations taken.
void enable_iteration_restart(const unsigned &restart)
switches on iteration restarting and takes as an argument the number of iterations after which the co...
virtual ~GMRES()
Destructor (cleanup storage)
unsigned Restart
The number of iterations before the iteration proceedure is restarted if iteration restart is used.
void apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update:
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void set_preconditioner_LHS()
Set left preconditioning (the default)
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void generate_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Generate a plane rotation. This is done by finding the values of (i....
void disable_iteration_restart()
switches off iteration restart
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
bool iteration_restart() const
access function indicating whether restarted GMRES is used
bool Iteration_restart
boolean indicating if iteration restarting is used
void smoother_solve(const DoubleVector &rhs, DoubleVector &result)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
unsigned iterations() const
Number of iterations taken.
GS(const GS &)=delete
Broken copy constructor.
unsigned Iterations
Number of iterations taken.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void operator=(const GS &)=delete
Broken assignment operator.
void clean_up_memory()
Clean up data that's stored for resolve (if any has been stored)
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
CRDoubleMatrix * Matrix_pt
System matrix pointer in the format specified by the template argument.
void smoother_setup(DoubleMatrixBase *matrix_pt)
Set up the smoother for the matrix specified by the pointer.
Vector< int > Index_of_diagonal_entries
Vector whose i'th entry contains the index of the last entry below or on the diagonal of the i'th row...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
double preconditioner_setup_time() const
Returns the time taken to set up the preconditioner.
virtual ~GS()
Destructor (cleanup storage)
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
virtual ~GS()
Destructor (cleanup storage)
unsigned Iterations
Number of iterations taken.
double preconditioner_setup_time() const
Returns the time taken to set up the preconditioner.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
void smoother_solve(const DoubleVector &rhs, DoubleVector &result)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void smoother_setup(DoubleMatrixBase *matrix_pt)
Set up the smoother for the matrix specified by the pointer.
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
MATRIX * Matrix_pt
System matrix pointer in the format specified by the template argument.
void operator=(const GS &)=delete
Broken assignment operator.
GS()
Constructor.
unsigned iterations() const
Number of iterations taken.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
GS(const GS &)=delete
Broken copy constructor.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here....
The Identity Preconditioner.
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
bool Use_iterative_solver_as_preconditioner
Use the iterative solver as preconditioner.
bool First_time_solve_when_used_as_preconditioner
When the iterative solver is used a preconditioner then we call the setup of solver method only once ...
void open_convergence_history_file_stream(const std::string &file_name, const std::string &zone_title="")
Write convergence history into file with specified filename (automatically switches on doc)....
void close_convergence_history_file_stream()
Close convergence history output stream.
Preconditioner *const & preconditioner_pt() const
Access function to preconditioner (const version)
void enable_setup_preconditioner_before_solve()
Setup the preconditioner before the solve.
double Tolerance
Convergence tolerance.
void disable_error_after_max_iter()
Don't throw an error if we don't converge within max_iter (default).
Preconditioner * Preconditioner_pt
Pointer to the preconditioner.
void enable_doc_convergence_history()
Enable documentation of the convergence history.
void enable_iterative_solver_as_preconditioner()
Enables the iterative solver be used as preconditioner (when calling the solve method it bypass the s...
virtual unsigned iterations() const =0
Number of iterations taken.
double linear_solver_solution_time() const
return the time taken to solve the linear system
virtual ~IterativeLinearSolver()
Destructor (empty)
void enable_error_after_max_iter()
Throw an error if we don't converge within max_iter.
bool Throw_error_after_max_iter
Should we throw an error instead of just returning when we hit the max iterations?
IterativeLinearSolver()
Constructor: Set (default) trivial preconditioner and set defaults for tolerance and max....
void disable_doc_convergence_history()
Disable documentation of the convergence history.
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
double Preconditioner_setup_time
Preconditioner setup time.
double Jacobian_setup_time
Jacobian setup time.
static IdentityPreconditioner Default_preconditioner
Default preconditioner: The base class for preconditioners is a fully functional (if trivial!...
virtual double preconditioner_setup_time() const
returns the the time taken to setup the preconditioner
void disable_iterative_solver_as_preconditioner()
Disables the iterative solver be used as preconditioner (when calling the solve method it bypass the ...
unsigned Max_iter
Maximum number of iterations.
double & tolerance()
Access to convergence tolerance.
double Solution_time
linear solver solution time
void operator=(const IterativeLinearSolver &)=delete
Broken assignment operator.
IterativeLinearSolver(const IterativeLinearSolver &)=delete
Broken copy constructor.
std::ofstream Output_file_stream
Output file stream for convergence history.
void disable_setup_preconditioner_before_solve()
Don't set up the preconditioner before the solve.
unsigned & max_iter()
Access to max. number of iterations.
bool Doc_convergence_history
Flag indicating if the convergence history is to be documented.
double jacobian_setup_time() const
returns the time taken to assemble the jacobian matrix and residual vector
bool Setup_preconditioner_before_solve
indicates whether the preconditioner should be setup before solve. Default = true;
unsigned nrow() const
access function to the number of global rows.
Base class for all linear solvers. This merely defines standard interfaces for linear solvers,...
Definition: linear_solver.h:68
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:73
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
bool Compute_gradient
flag that indicates whether the gradient required for the globally convergent Newton method should be...
Definition: linear_solver.h:81
An OomphLibError object which should be thrown when an run-time error is encountered....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual ~Smoother()
Virtual empty destructor.
void check_validity_of_solve_helper_inputs(MATRIX *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution, const double &n_dof)
Self-test to check that all the dimensions of the inputs to solve helper are consistent and everythin...
virtual void smoother_setup(DoubleMatrixBase *matrix_pt)=0
Set up the smoother for the matrix specified by the pointer.
virtual void smoother_solve(const DoubleVector &rhs, DoubleVector &result)=0
The smoother_solve function performs fixed number of iterations on the system A*result=rhs....
bool Use_as_smoother
When a derived class object is being used as a smoother in the MG solver (or elsewhere) the residual ...
Smoother()
Empty constructor.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...