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// Version 1.0; svn revision $LastChangedRevision: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //This header defines a class for linear iterative solvers
31 
32 //Include guards
33 #ifndef OOMPH_ITERATIVE_LINEAR_SOLVER_HEADER
34 #define OOMPH_ITERATIVE_LINEAR_SOLVER_HEADER
35 
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39 #include <oomph-lib-config.h>
40 #endif
41 
42 
43 //oomph-lib headers
44 #include "matrices.h"
45 #include "problem.h"
46 #include "linear_solver.h"
47 #include "preconditioner.h"
48 
49 
50 
51 
52 namespace oomph
53 {
54 
55 //=============================================================================
56 /// \short Base class for all linear iterative solvers.
57 /// This merely defines standard interfaces for linear iterative solvers,
58 /// so that different solvers can be used in a clean and transparent manner.
59 //=============================================================================
61  {
62 
63  public:
64 
65  /// \short Constructor: Set (default) trivial preconditioner and set
66  /// defaults for tolerance and max. number of iterations
68  {
69  //Set pointer to default preconditioner
71 
72  //Set default convergence tolerance
73  Tolerance=1.0e-8;
74 
75  //Set maximum number of iterations
76  Max_iter=100;
77 
78  // set default for document convergence history
80 
81  // set default
83 
85 
86  // By default the iterative solver is not used as preconditioner
88 
89  // Indicates whether this is the first time we call the solve
90  // method
92  }
93 
94  /// Broken copy constructor
96  {
97  BrokenCopy::broken_copy("IterativeLinearSolver");
98  }
99 
100  /// Broken assignment operator
102  {
103  BrokenCopy::broken_assign("IterativeLinearSolver");
104  }
105 
106  /// Destructor (empty)
108 
109  /// Access function to preconditioner
111 
112  /// Access function to preconditioner (const version)
114 
115  /// Access to convergence tolerance
116  double& tolerance() {return Tolerance;}
117 
118  /// Access to max. number of iterations
119  unsigned& max_iter() {return Max_iter;}
120 
121  /// Number of iterations taken
122  virtual unsigned iterations() const = 0;
123 
124  /// Enable documentation of the convergence history
126 
127  /// Disable documentation of the convergence history
129 
130  /// \short Write convergence history into file with specified filename
131  /// (automatically switches on doc). Optional second argument is a string
132  /// that can be used (as a zone title) to identify what case
133  /// we're running (e.g. what combination of linear solver and
134  /// preconditioner or parameter values are used).
136  const std::string& zone_title="")
137  {
138  // start docing
140 
141  // Close if it's open
142  if (Output_file_stream.is_open())
143  {
144  Output_file_stream.close();
145  }
146 
147  // Open new one
148  Output_file_stream.open(file_name.c_str());
149 
150  // Write tecplot zone header
151  Output_file_stream << "VARIABLES=\"iterations\", \"scaled residual\""
152  << std::endl;
153  Output_file_stream << "ZONE T=\"" << zone_title << "\""
154  << std::endl;
155  Output_file_stream << 0 << " " << 1.0 << std::endl;
156  }
157 
158  /// Close convergence history output stream
160  {
161  if (Output_file_stream.is_open()) Output_file_stream.close();
162  }
163 
164  /// \short returns the time taken to assemble the jacobian matrix and
165  /// residual vector
166  double jacobian_setup_time() const
167  {
168  return Jacobian_setup_time;
169  }
170 
171  /// \short return the time taken to solve the linear system
173  {
174  return Solution_time;
175  }
176 
177  /// \short returns the the time taken to setup the preconditioner
178  virtual double preconditioner_setup_time() const
179  {
181  }
182 
183  /// Setup the preconditioner before the solve
186 
187  /// Don't set up the preconditioner before the solve
190 
191  /// Throw an error if we don't converge within max_iter
194 
195  /// Don't throw an error if we don't converge within max_iter (default).
197  {Throw_error_after_max_iter = false;}
198 
199  /// Enables the iterative solver be used as preconditioner (when
200  /// calling the solve method it bypass the setup solver method --
201  /// currently only used by Trilinos solver ---)
204 
205  /// Disables the iterative solver be used as preconditioner (when
206  /// calling the solve method it bypass the setup solver method --
207  /// currently only used by Trilinos solver ---)
210 
211  protected:
212 
213  /// \short Flag indicating if the convergence history is to be
214  /// documented
216 
217  /// Output file stream for convergence history
218  std::ofstream Output_file_stream;
219 
220  /// \short Default preconditioner: The base
221  /// class for preconditioners is a fully functional (if trivial!)
222  /// preconditioner.
224 
225  ///Convergence tolerance
226  double Tolerance;
227 
228  ///Maximum number of iterations
229  unsigned Max_iter;
230 
231  /// Pointer to the preconditioner
233 
234  /// Jacobian setup time
236 
237  /// linear solver solution time
239 
240  /// Preconditioner setup time
242 
243  /// \short indicates whether the preconditioner should be setup before solve.
244  /// Default = true;
246 
247  /// \short Should we throw an error instead of just returning when we hit
248  /// the max iterations?
250 
251  /// \short Use the iterative solver as preconditioner
253 
254  /// When the iterative solver is used a preconditioner then we call
255  /// the setup of solver method only once (the first time the solve
256  /// method is called)
258  };
259 
260 
261 ///////////////////////////////////////////////////////////////////////////////
262 ///////////////////////////////////////////////////////////////////////////////
263 ///////////////////////////////////////////////////////////////////////////////
264 
265 
266 //======================================================================
267 /// \short The conjugate gradient method.
268 //======================================================================
269  template<typename MATRIX>
270  class CG : public IterativeLinearSolver
271  {
272 
273  public:
274 
275  ///Constructor
276  CG() : Iterations(0), Matrix_pt(0), Resolving(false),
278  {}
279 
280 
281  /// Destructor (cleanup storage)
282  virtual ~CG()
283  {
284  clean_up_memory();
285  }
286 
287  /// Broken copy constructor
288  CG(const CG&)
289  {
291  }
292 
293  /// Broken assignment operator
294  void operator=(const CG&)
295  {
297  }
298 
299  /// Overload disable resolve so that it cleans up memory too
301  {
303  clean_up_memory();
304  }
305 
306 
307  /// \short Solver: Takes pointer to problem and returns the results vector
308  /// which contains the solution of the linear system defined by
309  /// the problem's fully assembled Jacobian and residual vector.
310  void solve(Problem* const &problem_pt, DoubleVector &result);
311 
312  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
313  /// vector and returns the solution of the linear system.
314  void solve(DoubleMatrixBase* const &matrix_pt,
315  const DoubleVector &rhs,
316  DoubleVector &solution)
317  {
318  // Store the matrix if required
319  if ((Enable_resolve)&&(!Resolving))
320  {
321  Matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
322 
323  // Matrix has been passed in from the outside so we must not
324  // delete it
325  Matrix_can_be_deleted=false;
326  }
327 
328  // set the distribution
329  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt))
330  {
331  // the solver has the same distribution as the matrix if possible
332  this->build_distribution(dynamic_cast<DistributableLinearAlgebraObject*>
333  (matrix_pt)->distribution_pt());
334  }
335  else
336  {
337  // the solver has the same distribution as the RHS
338  this->build_distribution(rhs.distribution_pt());
339  }
340 
341  // Call the helper function
342  this->solve_helper(matrix_pt,rhs,solution);
343  }
344 
345  /// \short Re-solve the system defined by the last assembled Jacobian
346  /// and the rhs vector specified here. Solution is returned in the
347  /// vector result.
348  void resolve(const DoubleVector &rhs, DoubleVector &result);
349 
350  /// Number of iterations taken
351  unsigned iterations() const
352  {
353  return Iterations;
354  }
355 
356 
357  private:
358 
359  /// General interface to solve function
360  void solve_helper(DoubleMatrixBase* const &matrix_pt,
361  const DoubleVector &rhs,
362  DoubleVector &solution);
363 
364 
365  /// Cleanup data that's stored for resolve (if any has been stored)
367  {
368  if ((Matrix_pt!=0)&&(Matrix_can_be_deleted))
369  {
370  delete Matrix_pt;
371  Matrix_pt=0;
372  }
373  }
374 
375  /// Number of iterations taken
376  unsigned Iterations;
377 
378  /// Pointer to matrix
379  MATRIX* Matrix_pt;
380 
381  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
382  /// bypassing setup of matrix and preconditioner
383  bool Resolving;
384 
385  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
386  /// can be deleted.
388  };
389 
390 
391 ///////////////////////////////////////////////////////////////////////////////
392 ///////////////////////////////////////////////////////////////////////////////
393 ///////////////////////////////////////////////////////////////////////////////
394 
395 
396 //======================================================================
397 /// \short The conjugate gradient method.
398 //======================================================================
399  template<typename MATRIX>
401  {
402 
403  public:
404 
405  ///Constructor
408  {}
409 
410 
411  /// Destructor (cleanup storage)
412  virtual ~BiCGStab()
413  {
414  clean_up_memory();
415  }
416 
417  /// Broken copy constructor
419  {
420  BrokenCopy::broken_copy("BiCGStab");
421  }
422 
423  /// Broken assignment operator
424  void operator=(const BiCGStab&)
425  {
426  BrokenCopy::broken_assign("BiCGStab");
427  }
428 
429 
430 
431  /// Overload disable resolve so that it cleans up memory too
433  {
435  clean_up_memory();
436  }
437 
438  /// \short Solver: Takes pointer to problem and returns the results vector
439  /// which contains the solution of the linear system defined by
440  /// the problem's fully assembled Jacobian and residual vector.
441  void solve(Problem* const &problem_pt, DoubleVector &result);
442 
443  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
444  /// vector and returns the solution of the linear system.
445  void solve(DoubleMatrixBase* const &matrix_pt,
446  const DoubleVector& rhs,
447  DoubleVector &solution)
448  {
449  // Store the matrix if required
450  if ((Enable_resolve)&&(!Resolving))
451  {
452  Matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
453 
454  // Matrix has been passed in from the outside so we must not
455  // delete it
456  Matrix_can_be_deleted=false;
457  }
458 
459  // set the distribution
460  if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt))
461  {
462  // the solver has the same distribution as the matrix if possible
463  this->build_distribution(dynamic_cast<DistributableLinearAlgebraObject*>
464  (matrix_pt)->distribution_pt());
465  }
466  else
467  {
468  // the solver has the same distribution as the RHS
469  this->build_distribution(rhs.distribution_pt());
470  }
471 
472  //Call the helper function
473  this->solve_helper(matrix_pt,rhs,solution);
474  }
475 
476  /// \short Linear-algebra-type solver: Takes pointer to a matrix
477  /// and rhs vector and returns the solution of the linear system
478  /// Call the broken base-class version. If you want this, please
479  /// implement it
480  void solve(DoubleMatrixBase* const &matrix_pt,
481  const Vector<double> &rhs,
482  Vector<double> &result)
483  {LinearSolver::solve(matrix_pt,rhs,result);}
484 
485 
486 
487  /// \short Re-solve the system defined by the last assembled Jacobian
488  /// and the rhs vector specified here. Solution is returned in the
489  /// vector result.
490  void resolve(const DoubleVector &rhs,
491  DoubleVector &result);
492 
493 
494  /// Number of iterations taken
495  unsigned iterations() const
496  {
497  return Iterations;
498  }
499 
500 
501  private:
502 
503  /// General interface to solve function
504  void solve_helper(DoubleMatrixBase* const &matrix_pt,
505  const DoubleVector &rhs,
506  DoubleVector &solution);
507 
508  /// Cleanup data that's stored for resolve (if any has been stored)
510  {
511  if ((Matrix_pt!=0)&&(Matrix_can_be_deleted))
512  {
513  delete Matrix_pt;
514  Matrix_pt=0;
515  }
516  }
517 
518  /// Number of iterations taken
519  unsigned Iterations;
520 
521  /// Pointer to matrix
522  MATRIX* Matrix_pt;
523 
524  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
525  /// bypassing setup of matrix and preconditioner
526  bool Resolving;
527 
528  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
529  /// can be deleted.
531 
532  };
533 
534 
535 ///////////////////////////////////////////////////////////////////////
536 ///////////////////////////////////////////////////////////////////////
537 ///////////////////////////////////////////////////////////////////////
538 
539 
540 //====================================================================
541 /// Smoother class:
542 /// The smoother class is designed for to be used in conjunction with
543 /// multigrid. The action of the smoother should reduce the high
544 /// frequency errors. These methods are inefficient as stand-alone
545 /// solvers.
546 //====================================================================
548  {
549 
550  public:
551 
552  /// Empty constructor
554  {}
555 
556  /// Virtual empty destructor
557  virtual ~Smoother(){};
558 
559  /// \short The smoother_solve function performs fixed number of iterations
560  /// on the system A*result=rhs. The number of (smoothing) iterations is
561  /// the same as the max. number of iterations in the underlying
562  /// IterativeLinearSolver class. Note that the result vector MUST NOT
563  /// re-initialised to zero (as it would typically be when the Smoother is called as a
564  /// iterative linear solver).
565  virtual void smoother_solve(const DoubleVector& rhs,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  /// \short 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 
580  /// \short When a derived class object is being used as a smoother in
581  /// the MG solver (or elsewhere) the residual norm does not need to be calculated
582  /// because we're simply performing a fixed number of (smoothing) iterations.
583  /// This boolean is used as a flag to indicate that the IterativeLinearSolver
584  /// (which this class is by inheritance) is supposed to act in this way.
586  }; // End of Smoother
587 
588 
589 ///////////////////////////////////////////////////////////////////////////
590 ///////////////////////////////////////////////////////////////////////////
591 ///////////////////////////////////////////////////////////////////////////
592 
593 
594 //=========================================================================
595 /// \short The Gauss Seidel method
596 //=========================================================================
597  template<typename MATRIX>
598  class GS : public virtual Smoother
599  {
600 
601  public:
602 
603  /// Constructor
604  GS() : Matrix_pt(0), Iterations(0), Resolving(false),
606  {}
607 
608  /// Destructor (cleanup storage)
609  virtual ~GS()
610  {
611  clean_up_memory();
612  }
613 
614  /// Broken copy constructor
615  GS(const GS&)
616  {
618  }
619 
620  /// Broken assignment operator
621  void operator=(const GS&)
622  {
624  }
625 
626  /// Overload disable resolve so that it cleans up memory too
628  {
630  clean_up_memory();
631  } // End of disable_resolve
632 
633  /// Set up the smoother for the matrix specified by the pointer
635  {
636  // Assume the matrix has been passed in from the outside so we must
637  // not delete it. This is needed to avoid pre- and post-smoothers
638  // deleting the same matrix in the MG solver. If this was originally
639  // set to TRUE then this will be sorted out in the other functions
640  // from which this was called
641  Matrix_can_be_deleted=false;
642 
643  // Upcast the input matrix to system matrix to the type MATRIX
644  Matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
645  } // End of smoother_setup
646 
647  /// \short The smoother_solve function performs fixed number of iterations
648  /// on the system A*result=rhs. The number of (smoothing) iterations is
649  /// the same as the max. number of iterations in the underlying
650  /// IterativeLinearSolver class.
651  void smoother_solve(const DoubleVector& rhs, DoubleVector& result)
652  {
653  // If you use a smoother but you don't want to calculate the residual
654  Use_as_smoother=true;
655 
656  // Call the helper function
657  solve_helper(Matrix_pt,rhs,result);
658  } // End of smoother_setup
659 
660  /// \short Solver: Takes pointer to problem and returns the results
661  /// vector which contains the solution of the linear system defined
662  /// by the problem's fully assembled Jacobian and residual vector.
663  void solve(Problem* const &problem_pt, DoubleVector &result);
664 
665  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
666  /// vector and returns the solution of the linear system.
667  void solve(DoubleMatrixBase* const &matrix_pt,
668  const DoubleVector &rhs,
669  DoubleVector &solution)
670  {
671  // Reset the Use_as_smoother_flag as the solver is not being used
672  // as a smoother
673  Use_as_smoother=false;
674 
675  // Set up the distribution
676  this->build_distribution(rhs.distribution_pt());
677 
678  // Store the matrix if required
679  if ((Enable_resolve)&&(!Resolving))
680  {
681  // Upcast to the appropriate matrix type
682  Matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
683  }
684 
685  // Matrix has been passed in from the outside so we must not delete it
686  Matrix_can_be_deleted=false;
687 
688  // Call the helper function
689  this->solve_helper(matrix_pt,rhs,solution);
690  } // End of solve
691 
692  /// \short Linear-algebra-type solver: Takes pointer to a matrix
693  /// and rhs vector and returns the solution of the linear system
694  /// Call the broken base-class version. If you want this, please
695  /// implement it
696  void solve(DoubleMatrixBase* const &matrix_pt,
697  const Vector<double> &rhs,
698  Vector<double> &result)
699  {
700  LinearSolver::solve(matrix_pt,rhs,result);
701  } // End of solve
702 
703  /// \short Re-solve the system defined by the last assembled Jacobian
704  /// and the rhs vector specified here. Solution is returned in the
705  /// vector result.
706  void resolve(const DoubleVector &rhs, DoubleVector &result)
707  {
708  // We are re-solving
709  Resolving=true;
710 
711 #ifdef PARANOID
712  if (Matrix_pt==0)
713  {
714  throw OomphLibError("No matrix was stored -- cannot re-solve",
715  OOMPH_CURRENT_FUNCTION,
716  OOMPH_EXCEPTION_LOCATION);
717  }
718 #endif
719 
720  // Call linear algebra-style solver
721  solve(Matrix_pt,rhs,result);
722 
723  // Reset re-solving flag
724  Resolving=false;
725  } // End of resolve
726 
727  /// Returns the time taken to set up the preconditioner
729  {
730  throw OomphLibError(
731  "Gauss Seidel is not a preconditionable iterative solver",
732  OOMPH_CURRENT_FUNCTION,
733  OOMPH_EXCEPTION_LOCATION);
734  return 0;
735  } // End of preconditioner_setup_time
736 
737  /// Number of iterations taken
738  unsigned iterations() const
739  {
740  return Iterations;
741  } // End of iterations
742 
743  private:
744 
745  /// General interface to solve function
746  void solve_helper(DoubleMatrixBase* const &matrix_pt,
747  const DoubleVector& rhs,
748  DoubleVector& solution);
749 
750  /// Cleanup data that's stored for resolve (if any has been stored)
752  {
753  // If the matrix pointer isn't null and we're allowed to delete it
754  // delete the matrix and assign the pointer the value NULL
755  if ((Matrix_pt!=0)&&(Matrix_can_be_deleted))
756  {
757  // Destroy the matrix
758  delete Matrix_pt;
759 
760  // Make it a null pointer
761  Matrix_pt=0;
762  }
763  } // End of clean_up_memory
764 
765  /// System matrix pointer in the format specified by the template argument
766  MATRIX* Matrix_pt;
767 
768  /// Number of iterations taken
769  unsigned Iterations;
770 
771  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
772  /// bypassing setup of matrix and preconditioner
773  bool Resolving;
774 
775  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
776  /// can be deleted.
778 
779  };
780 
781 ///////////////////////////////////////////////////////////////////////////
782 ///////////////////////////////////////////////////////////////////////////
783 ///////////////////////////////////////////////////////////////////////////
784 
785 //=========================================================================
786 /// \short Explicit template specialisation of the Gauss Seidel method for
787 /// compressed row format matrices
788 //=========================================================================
789  template<>
790  class GS<CRDoubleMatrix> : public virtual Smoother
791  {
792  public:
793 
794  /// Constructor
795  GS() : Matrix_pt(0), Iterations(0), Resolving(false),
797  {}
798 
799  /// Destructor (cleanup storage)
800  virtual ~GS()
801  {
802  clean_up_memory();
803  }
804 
805  /// Broken copy constructor
806  GS(const GS&)
807  {
809  }
810 
811  /// Broken assignment operator
812  void operator=(const GS&)
813  {
815  }
816 
817  /// \short The smoother_solve function performs fixed number of iterations
818  /// on the system A*result=rhs. The number of (smoothing) iterations is
819  /// the same as the max. number of iterations in the underlying
820  /// IterativeLinearSolver class.
821  void smoother_solve(const DoubleVector& rhs, DoubleVector& result)
822  {
823  // If you use a smoother but you don't want to calculate the residual
824  Use_as_smoother=true;
825 
826  // Call the helper function
827  solve_helper(Matrix_pt,rhs,result);
828  } // End of smoother_solve
829 
830  /// Overload disable resolve so that it cleans up memory too
832  {
834  clean_up_memory();
835  } // End of disable_resolve
836 
837  /// Set up the smoother for the matrix specified by the pointer
839  {
840  // Assume the matrix has been passed in from the outside so we must
841  // not delete it. This is needed to avoid pre- and post-smoothers
842  // deleting the same matrix in the MG solver. If this was originally
843  // set to TRUE then this will be sorted out in the other functions
844  // from which this was called
845  Matrix_can_be_deleted=false;
846 
847  // Call the generic setup helper function
848  setup_helper(matrix_pt);
849  } // End of smoother_setup
850 
851  /// \short Generic setup function to sort out everything that needs to be
852  /// set up with regards to the input matrix
853  void setup_helper(DoubleMatrixBase* matrix_pt);
854 
855  /// \short Solver: Takes pointer to problem and returns the results vector
856  /// which contains the solution of the linear system defined by
857  /// the problem's fully assembled Jacobian and residual vector.
858  void solve(Problem* const &problem_pt, DoubleVector &result);
859 
860  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
861  /// vector and returns the solution of the linear system.
862  void solve(DoubleMatrixBase* const &matrix_pt,
863  const DoubleVector& rhs,
864  DoubleVector& solution)
865  {
866  // Reset the Use_as_smoother_flag as the solver is not being used
867  // as a smoother
868  Use_as_smoother=false;
869 
870  // Set up the distribution
871  this->build_distribution(rhs.distribution_pt());
872 
873  // Store the matrix if required
874  if ((Enable_resolve)&&(!Resolving))
875  {
876  // Upcast to the appropriate matrix type and sort the matrix entries
877  // out so that the CRDoubleMatrix implementation of the Gauss-Seidel
878  // solver can be used
879  Matrix_pt=dynamic_cast<CRDoubleMatrix*>(matrix_pt);
880  }
881  // We still need to sort the entries
882  else
883  {
884  // The system matrix here is a CRDoubleMatrix. To make use of the
885  // specific implementation of the solver for this type of matrix we
886  // need to make sure the entries are arranged correctly
887  dynamic_cast<CRDoubleMatrix*>(matrix_pt)->sort_entries();
888 
889  // Now get access to the vector Index_of_diagonal_entries
890  Index_of_diagonal_entries=dynamic_cast<CRDoubleMatrix*>
891  (matrix_pt)->get_index_of_diagonal_entries();
892  }
893 
894  // Matrix has been passed in from the outside so we must not delete it
895  Matrix_can_be_deleted=false;
896 
897  // Call the helper function
898  solve_helper(matrix_pt,rhs,solution);
899  } // End of solve
900 
901  /// \short Linear-algebra-type solver: Takes pointer to a matrix
902  /// and rhs vector and returns the solution of the linear system
903  /// Call the broken base-class version. If you want this, please
904  /// implement it
905  void solve(DoubleMatrixBase* const &matrix_pt,
906  const Vector<double> &rhs,
907  Vector<double> &result)
908  {
909  LinearSolver::solve(matrix_pt,rhs,result);
910  } // End of solve
911 
912  /// \short Re-solve the system defined by the last assembled Jacobian
913  /// and the rhs vector specified here. Solution is returned in the
914  /// vector result.
915  void resolve(const DoubleVector& rhs,DoubleVector& result)
916  {
917  // We are re-solving
918  Resolving=true;
919 
920 #ifdef PARANOID
921  // If the matrix pointer is null
922  if (this->Matrix_pt==0)
923  {
924  throw OomphLibError("No matrix was stored -- cannot re-solve",
925  OOMPH_CURRENT_FUNCTION,
926  OOMPH_EXCEPTION_LOCATION);
927  }
928 #endif
929 
930  // Call linear algebra-style solver
931  solve(Matrix_pt,rhs,result);
932 
933  // Reset re-solving flag
934  Resolving=false;
935  } // End of resolve
936 
937  /// Returns the time taken to set up the preconditioner
939  {
940  // Create the error message
941  std::string error_output_string="Gauss Seidel is not a ";
942  error_output_string+="preconditionable iterative solver";
943 
944  // Throw an error
945  throw OomphLibError(error_output_string,
946  OOMPH_CURRENT_FUNCTION,
947  OOMPH_EXCEPTION_LOCATION);
948 
949  // Return a value so the compiler doesn't throw up an error about
950  // no input being returned
951  return 0;
952  } // End of preconditioner_setup_time
953 
954  /// Number of iterations taken
955  unsigned iterations() const
956  {
957  // Return the number of iterations
958  return Iterations;
959  } // End of iterations
960 
961  private:
962 
963  /// General interface to solve function
964  void solve_helper(DoubleMatrixBase* const &matrix_pt,
965  const DoubleVector& rhs,
966  DoubleVector& solution);
967 
968  /// Clean up data that's stored for resolve (if any has been stored)
970  {
971  // If the matrix pointer isn't null AND we're allowed to delete the
972  // matrix which is only when we create the matrix ourselves
973  if ((Matrix_pt!=0)&&(Matrix_can_be_deleted))
974  {
975  // Delete the matrix
976  delete Matrix_pt;
977 
978  // Assign the associated pointer the value NULL
979  Matrix_pt=0;
980  }
981  } // End of clean_up_memory
982 
983  /// System matrix pointer in the format specified by the template argument
985 
986  /// Number of iterations taken
987  unsigned Iterations;
988 
989  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
990  /// bypassing setup of matrix and preconditioner
991  bool Resolving;
992 
993  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
994  /// can be deleted.
996 
997  /// \short Vector whose i'th entry contains the index of the last entry
998  /// below or on the diagonal of the i'th row of the matrix
1000  };
1001 
1002 ///////////////////////////////////////////////////////////////////////
1003 ///////////////////////////////////////////////////////////////////////
1004 ///////////////////////////////////////////////////////////////////////
1005 
1006 //=========================================================================
1007 /// Damped Jacobi "solver" templated by matrix type. The "solver"
1008 /// exists in many different incarnations: It's an IterativeLinearSolver,
1009 /// and a Smoother, all of which use the same basic iteration.
1010 //=========================================================================
1011  template<typename MATRIX>
1012  class DampedJacobi : public virtual Smoother
1013  {
1014 
1015  public:
1016 
1017  /// Empty constructor
1018  DampedJacobi(const double& omega=2.0/3.0) : Matrix_can_be_deleted(true)
1019  {
1020  // Damping factor
1021  Omega=omega;
1022  }
1023 
1024  /// Empty destructor
1026  {
1027  // Run the generic clean up function
1028  clean_up_memory();
1029  }
1030 
1031  /// Broken copy constructor
1033  {
1034  BrokenCopy::broken_copy("DampedJacobi");
1035  }
1036 
1037  /// Broken assignment operator
1039  {
1040  BrokenCopy::broken_assign("DampedJacobi");
1041  }
1042 
1043  /// Cleanup data that's stored for resolve (if any has been stored)
1045  {
1046  // If the matrix pointer isn't null AND we're allowed to delete the
1047  // matrix which is only when we create the matrix ourselves
1048  if ((Matrix_pt!=0) && (Matrix_can_be_deleted))
1049  {
1050  // Delete the matrix
1051  delete Matrix_pt;
1052 
1053  // Assign the associated pointer the value NULL
1054  Matrix_pt=0;
1055  }
1056  } // End of clean_up_memory
1057 
1058  /// Setup: Pass pointer to the matrix and store in cast form
1060  {
1061  // Assume the matrix has been passed in from the outside so we must not
1062  // delete it
1063  Matrix_can_be_deleted=false;
1064 
1065  // Upcast to the appropriate matrix type
1066  Matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
1067 
1068  // Extract the diagonal entries of the matrix and store them
1069  extract_diagonal_entries(matrix_pt);
1070  } // End of smoother_setup
1071 
1072  /// Function to extract the diagonal entries from the matrix
1074  {
1075  // If we're using a CRDoubleMatrix object
1076  if (dynamic_cast<CRDoubleMatrix*>(matrix_pt))
1077  {
1078  // The matrix diagonal (we need this when we need to calculate inv(D)
1079  // where D is the diagonal of A and it remains the same for all uses
1080  // of the iterative scheme so we can store it and call it in each
1081  // iteration)
1082  Matrix_diagonal=dynamic_cast<CRDoubleMatrix*>
1083  (Matrix_pt)->diagonal_entries();
1084  }
1085  // If we're using a complex matrix then diagonal entries has to be a
1086  // complex vector rather than a vector of doubles.
1087  else if (dynamic_cast<CCDoubleMatrix*>(matrix_pt))
1088  {
1089  // Make an ostringstream object to create an error message
1090  std::ostringstream error_message_stream;
1091 
1092  // Create the error message
1093  error_message_stream << "Damped Jacobi can only cater to real-valued "
1094  << "matrices. If you require a complex-valued "
1095  << "version, please write this yourself. "
1096  << "It is likely that the only difference will be "
1097  << "the use of complex vectors.";
1098 
1099  // Throw an error
1100  throw OomphLibError(error_message_stream.str(),
1101  OOMPH_CURRENT_FUNCTION,
1102  OOMPH_EXCEPTION_LOCATION);
1103  }
1104  // Just extract the diagonal entries normally
1105  else
1106  {
1107  // Calculate the number of rows in the matrix
1108  unsigned n_row=Matrix_pt->nrow();
1109 
1110  // Loop over the rows of the matrix
1111  for (unsigned i=0;i<n_row;i++)
1112  {
1113  // Assign the i-th value of Matrix_diagonal
1114  Matrix_diagonal[i]=(*Matrix_pt)(i,i);
1115  }
1116  } // if (dynamic_cast<CRDoubleMatrix*>(matrix_pt))
1117 
1118  // Calculate the n.d.o.f.
1119  unsigned n_dof=Matrix_diagonal.size();
1120 
1121  // Find the reciprocal of the entries of Matrix_diagonal
1122  for (unsigned i=0;i<n_dof;i++)
1123  {
1125  }
1126  } // End of extract_diagonal_entries
1127 
1128  /// \short The smoother_solve function performs fixed number of iterations
1129  /// on the system A*result=rhs. The number of (smoothing) iterations is
1130  /// the same as the max. number of iterations in the underlying
1131  /// IterativeLinearSolver class.
1132  void smoother_solve(const DoubleVector& rhs, DoubleVector& solution)
1133  {
1134  // If you use a smoother but you don't want to calculate the residual
1135  Use_as_smoother=true;
1136 
1137  // Call the helper function
1138  solve_helper(Matrix_pt,rhs,solution);
1139  } // End of smoother_solve
1140 
1141  /// \short Use damped Jacobi iteration as an IterativeLinearSolver:
1142  /// This obtains the Jacobian matrix J and the residual vector r
1143  /// (needed for the Newton method) from the problem's get_jacobian
1144  /// function and returns the result of Jx=r.
1145  void solve(Problem* const& problem_pt, DoubleVector& result);
1146 
1147  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
1148  /// vector and returns the solution of the linear system.
1149  void solve(DoubleMatrixBase* const &matrix_pt,
1150  const DoubleVector& rhs,
1151  DoubleVector& solution)
1152  {
1153  // Matrix has been passed in from the outside so we must not delete it
1154  Matrix_can_be_deleted=false;
1155 
1156  // Indicate that the solver is not being used as a smoother
1157  Use_as_smoother=false;
1158 
1159  // Set up the distribution
1160  this->build_distribution(rhs.distribution_pt());
1161 
1162  // Store the matrix if required
1163  if ((Enable_resolve)&&(!Resolving))
1164  {
1165  // Upcast to the appropriate matrix type
1166  Matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
1167  }
1168 
1169  // Extract the diagonal entries of the matrix and store them
1170  extract_diagonal_entries(matrix_pt);
1171 
1172  // Call the helper function
1173  solve_helper(matrix_pt,rhs,solution);
1174  } // End of solve
1175 
1176  /// \short Re-solve the system defined by the last assembled Jacobian
1177  /// and the rhs vector specified here. Solution is returned in the
1178  /// vector result.
1179  void resolve(const DoubleVector &rhs, DoubleVector &result)
1180  {
1181  // We are re-solving
1182  Resolving=true;
1183 
1184 #ifdef PARANOID
1185  if (Matrix_pt==0)
1186  {
1187  throw OomphLibError("No matrix was stored -- cannot re-solve",
1188  OOMPH_CURRENT_FUNCTION,
1189  OOMPH_EXCEPTION_LOCATION);
1190  }
1191 #endif
1192 
1193  // Call linear algebra-style solver
1194  solve(Matrix_pt,rhs,result);
1195 
1196  // Reset re-solving flag
1197  Resolving=false;
1198  } // End of resolve
1199 
1200  /// Number of iterations taken
1201  unsigned iterations() const
1202  {
1203  // Return the value of Iterations
1204  return Iterations;
1205  } // End of iterations
1206 
1207  private:
1208 
1209  /// \short This is where the actual work is done -- different
1210  /// implementations for different matrix types.
1211  void solve_helper(DoubleMatrixBase* const &matrix_pt,
1212  const DoubleVector &rhs,
1213  DoubleVector &solution);
1214 
1215  /// Pointer to the matrix
1216  MATRIX* Matrix_pt;
1217 
1218  /// Vector containing the diagonal entries of A
1220 
1221  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
1222  /// bypassing setup of matrix and preconditioner
1224 
1225  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
1226  /// can be deleted.
1228 
1229  /// Number of iterations taken
1230  unsigned Iterations;
1231 
1232  /// Damping factor
1233  double Omega;
1234  };
1235 
1236 
1237 ///////////////////////////////////////////////////////////////////////
1238 ///////////////////////////////////////////////////////////////////////
1239 ///////////////////////////////////////////////////////////////////////
1240 
1241 
1242 //======================================================================
1243 /// \short The GMRES method.
1244 //======================================================================
1245  template<typename MATRIX>
1247  {
1248 
1249  public:
1250 
1251  /// Constructor
1253  Matrix_pt(0),
1254  Resolving(false),
1255  Matrix_can_be_deleted(true)
1256  {
1257  Preconditioner_LHS=true;
1258  Iteration_restart=false;
1259  }
1260 
1261  /// Destructor (cleanup storage)
1262  virtual ~GMRES()
1263  {
1264  clean_up_memory();
1265  }
1266 
1267  /// Broken copy constructor
1268  GMRES(const GMRES&)
1269  {
1270  BrokenCopy::broken_copy("GMRES");
1271  }
1272 
1273  /// Broken assignment operator
1274  void operator=(const GMRES&)
1275  {
1276  BrokenCopy::broken_assign("GMRES");
1277  }
1278 
1279  /// Overload disable resolve so that it cleans up memory too
1281  {
1283  clean_up_memory();
1284  }
1285 
1286  /// \short Solver: Takes pointer to problem and returns the results vector
1287  /// which contains the solution of the linear system defined by
1288  /// the problem's fully assembled Jacobian and residual vector.
1289  void solve(Problem* const &problem_pt, DoubleVector &result);
1290 
1291  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
1292  /// vector and returns the solution of the linear system.
1293  void solve(DoubleMatrixBase* const &matrix_pt,
1294  const DoubleVector &rhs,
1295  DoubleVector &solution)
1296  {
1297  // setup the distribution
1298  this->build_distribution(rhs.distribution_pt());
1299 
1300  // Store the matrix if required
1301  if ((Enable_resolve)&&(!Resolving))
1302  {
1303  Matrix_pt=dynamic_cast<MATRIX*>(matrix_pt);
1304 
1305  // Matrix has been passed in from the outside so we must not
1306  // delete it
1307  Matrix_can_be_deleted=false;
1308  }
1309 
1310  // Call the helper function
1311  this->solve_helper(matrix_pt,rhs,solution);
1312  }
1313 
1314 
1315  /// \short Linear-algebra-type solver: Takes pointer to a matrix
1316  /// and rhs vector and returns the solution of the linear system
1317  /// Call the broken base-class version. If you want this, please
1318  /// implement it
1319  void solve(DoubleMatrixBase* const &matrix_pt,
1320  const Vector<double> &rhs,
1321  Vector<double> &result)
1322  {LinearSolver::solve(matrix_pt,rhs,result);}
1323 
1324  /// \short Re-solve the system defined by the last assembled Jacobian
1325  /// and the rhs vector specified here. Solution is returned in the
1326  /// vector result.
1327  void resolve(const DoubleVector &rhs,
1328  DoubleVector &result);
1329 
1330  /// Number of iterations taken
1331  unsigned iterations() const
1332  {
1333  return Iterations;
1334  }
1335 
1336  /// \short access function indicating whether restarted GMRES is used
1337  bool iteration_restart() const
1338  {
1339  return Iteration_restart;
1340  }
1341 
1342  /// \short switches on iteration restarting and takes as an argument the
1343  /// number of iterations after which the construction of the orthogonalisation
1344  /// basis vectors should be restarted
1345  void enable_iteration_restart(const unsigned& restart)
1346  {
1347  Restart = restart;
1348  Iteration_restart = true;
1349  }
1350 
1351  /// switches off iteration restart
1353  {
1354  Iteration_restart = false;
1355  }
1356 
1357  /// \short Set left preconditioning (the default)
1359 
1360  /// \short Enable right preconditioning
1362 
1363  private:
1364 
1365  /// General interface to solve function
1366  void solve_helper(DoubleMatrixBase* const &matrix_pt,
1367  const DoubleVector &rhs,
1368  DoubleVector &solution);
1369 
1370  /// Cleanup data that's stored for resolve (if any has been stored)
1372  {
1373  if ((Matrix_pt!=0)&&(Matrix_can_be_deleted))
1374  {
1375  delete Matrix_pt;
1376  Matrix_pt=0;
1377  }
1378  }
1379 
1380  /// Helper function to update the result vector using the result, x=x_0+V_m*y
1381  void update(const unsigned& k,const Vector<Vector<double> >& H,
1382  const Vector<double>& s,const Vector<DoubleVector>& v,
1383  DoubleVector& x)
1384  {
1385  // Make a local copy of s
1386  Vector<double> y(s);
1387 
1388  // Backsolve:
1389  for (int i=int(k);i>=0;i--)
1390  {
1391  // Divide the i-th entry of y by the i-th diagonal entry of H
1392  y[i]/=H[i][i];
1393 
1394  // Loop over the previous values of y and update
1395  for (int j=i-1;j>=0;j--)
1396  {
1397  // Update the j-th entry of y
1398  y[j] -= H[i][j]*y[i];
1399  }
1400  } // for (int i=int(k);i>=0;i--)
1401 
1402  // Store the number of rows in the result vector
1403  unsigned n_x=x.nrow();
1404 
1405  // Build a temporary vector with entries initialised to 0.0
1406  DoubleVector temp(x.distribution_pt(),0.0);
1407 
1408  // Build a temporary vector with entries initialised to 0.0
1409  DoubleVector z(x.distribution_pt(),0.0);
1410 
1411  // Get access to the underlying values
1412  double* temp_pt=temp.values_pt();
1413 
1414  // Calculate x=Vy
1415  for (unsigned j=0;j<=k;j++)
1416  {
1417  // Get access to j-th column of v
1418  const double* vj_pt=v[j].values_pt();
1419 
1420  // Loop over the entries of the vector, temp
1421  for (unsigned i=0;i<n_x;i++)
1422  {
1423  temp_pt[i]+=vj_pt[i]*y[j];
1424  }
1425  } // for (unsigned j=0;j<=k;j++)
1426 
1427  // If we're using LHS preconditioning
1428  if(Preconditioner_LHS)
1429  {
1430  // Since we're using LHS preconditioning the preconditioner is applied
1431  // to the matrix and RHS vector so we simply update the value of x
1432  x+=temp;
1433  }
1434  // If we're using RHS preconditioning
1435  else
1436  {
1437  // Since we're using RHS preconditioning the preconditioner is applied
1438  // to the solution vector
1440 
1441  // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
1442  // sparse linear systems", p.284]
1443  x+=z;
1444  }
1445  } // End of update
1446 
1447  /// \short Helper function: Generate a plane rotation. This is done by
1448  /// finding the values of \f$ \cos(\theta) \f$ (i.e. cs) and \sin(\theta)
1449  /// (i.e. sn) such that:
1450  /// \f[
1451  /// \begin{bmatrix}
1452  /// \cos\theta & \sin\theta \newline
1453  /// -\sin\theta & \cos\theta
1454  /// \end{bmatrix}
1455  /// \begin{bmatrix}
1456  /// dx \newline
1457  /// dy
1458  /// \end{bmatrix}
1459  /// =
1460  /// \begin{bmatrix}
1461  /// r \newline
1462  /// 0
1463  /// \end{bmatrix},
1464  /// \f]
1465  /// where \f$ r=\sqrt{pow(dx,2)+pow(dy,2)} \f$. The values of a and b are
1466  /// given by:
1467  /// \f[
1468  /// \cos\theta&=\dfrac{dx}{\sqrt{pow(dx,2)+pow(dy,2)}},
1469  /// \f]
1470  /// and
1471  /// \f[
1472  /// \sin\theta&=\dfrac{dy}{\sqrt{pow(dx,2)+pow(dy,2)}}.
1473  /// \f]
1474  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.192
1475  void generate_plane_rotation(double &dx,double &dy,double &cs,double &sn)
1476  {
1477  // If dy=0 then we do not need to apply a rotation
1478  if (dy==0.0)
1479  {
1480  // Using theta=0 gives cos(theta)=1
1481  cs=1.0;
1482 
1483  // Using theta=0 gives sin(theta)=0
1484  sn=0.0;
1485  }
1486  // If dx or dy is large using the normal form of calculting cs and sn
1487  // is naive since this may overflow or underflow so instead we calculate
1488  // r=sqrt(pow(dx,2)+pow(dy,2)) by using r=|dy|sqrt(1+pow(dx/dy,2)) if
1489  // |dy|>|dx| [see <A HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
1490  else if(fabs(dy)>fabs(dx))
1491  {
1492  // Since |dy|>|dx| calculate the ratio dx/dy
1493  double temp=dx/dy;
1494 
1495  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))
1496  sn=1.0/sqrt(1.0+temp*temp);
1497 
1498  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))=(dx/dy)*sin(theta)
1499  cs=temp*sn;
1500  }
1501  // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
1502  // calculate the values of cs and sn using the method above
1503  else
1504  {
1505  // Since |dx|>=|dy| calculate the ratio dy/dx
1506  double temp=dy/dx;
1507 
1508  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))
1509  cs=1.0/sqrt(1.0+temp*temp);
1510 
1511  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))=(dy/dx)*cos(theta)
1512  sn=temp*cs;
1513  }
1514  } // End of generate_plane_rotation
1515 
1516  /// \short Helper function: Apply plane rotation. This is done using the
1517  /// update:
1518  /// \f[
1519  ///\begin{bmatrix}
1520  /// dx \newline
1521  /// dy
1522  /// \end{bmatrix}
1523  /// \leftarrow
1524  /// \begin{bmatrix}
1525  /// \cos\theta & \sin\theta \newline
1526  /// -\sin\theta & \cos\theta
1527  /// \end{bmatrix}
1528  /// \begin{bmatrix}
1529  /// dx \newline
1530  /// dy
1531  /// \end{bmatrix}.
1532  /// \f]
1533  void apply_plane_rotation(double &dx,double &dy,double &cs,double &sn)
1534  {
1535  // Calculate the value of dx but don't update it yet
1536  double temp=cs*dx+sn*dy;
1537 
1538  // Set the value of dy
1539  dy=-sn*dx+cs*dy;
1540 
1541  // Set the value of dx using the correct values of dx and dy
1542  dx=temp;
1543  }
1544 
1545  /// Number of iterations taken
1546  unsigned Iterations;
1547 
1548  /// \short The number of iterations before the iteration proceedure is
1549  /// restarted if iteration restart is used
1550  unsigned Restart;
1551 
1552  /// boolean indicating if iteration restarting is used
1554 
1555  /// Pointer to matrix
1556  MATRIX* Matrix_pt;
1557 
1558  /// \short Boolean flag to indicate if the solve is done in re-solve mode,
1559  /// bypassing setup of matrix and preconditioner
1561 
1562  /// \short Boolean flag to indicate if the matrix pointed to be Matrix_pt
1563  /// can be deleted.
1565 
1566  /// \short boolean indicating use of left hand preconditioning (if true)
1567  /// or right hand preconditioning (if false)
1569  };
1570 } // End of namespace oomph
1571 
1572 #endif
virtual double preconditioner_setup_time() const
returns the the time taken to setup the preconditioner
bool Throw_error_after_max_iter
Should we throw an error instead of just returning when we hit the max iterations?
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
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...
GS(const GS &)
Broken copy constructor.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
void smoother_solve(const DoubleVector &rhs, DoubleVector &result)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs. The number of (smoothing) iterations is the same as the max. number of iterations in the underlying IterativeLinearSolver class.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
The Identity Preconditioner.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
static IdentityPreconditioner Default_preconditioner
Default preconditioner: The base class for preconditioners is a fully functional (if trivial!) precon...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
virtual unsigned iterations() const =0
Number of iterations taken.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
bool Iteration_restart
boolean indicating if iteration restarting is used
unsigned iterations() const
Number of iterations taken.
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
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.
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...
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
double & tolerance()
Access to convergence tolerance.
virtual ~BiCGStab()
Destructor (cleanup storage)
void operator=(const CG &)
Broken assignment operator.
void smoother_setup(DoubleMatrixBase *matrix_pt)
Setup: Pass pointer to the matrix and store in cast form.
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 enable_error_after_max_iter()
Throw an error if we don't converge within max_iter.
double Tolerance
Convergence tolerance.
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 ...
double * values_pt()
access function to the underlying values
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
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 clean_up_memory()
Clean up data that's stored for resolve (if any has been stored)
cstr elem_len * i
Definition: cfortran.h:607
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 apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update: .
unsigned iterations() const
Number of iterations taken.
virtual ~GS()
Destructor (cleanup storage)
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. The number of (smoothing) iterations is the same as the max. number of iterations in the underlying IterativeLinearSolver class. Note that the result vector MUST NOT re-initialised to zero (as it would typically be when the Smoother is called as a iterative linear solver).
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
unsigned Max_iter
Maximum number of iterations.
unsigned nrow() const
access function to the number of global rows.
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:80
The Problem class.
Definition: problem.h:152
Preconditioner * Preconditioner_pt
Pointer to the preconditioner.
~DampedJacobi()
Empty destructor.
virtual ~IterativeLinearSolver()
Destructor (empty)
MATRIX * Matrix_pt
Pointer to the matrix.
void operator=(const GS &)
Broken assignment operator.
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)
This is where the actual work is done – different implementations for different matrix types...
Preconditioner *const & preconditioner_pt() const
Access function to preconditioner (const version)
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 ...
IterativeLinearSolver()
Constructor: Set (default) trivial preconditioner and set defaults for tolerance and max...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
std::ofstream Output_file_stream
Output file stream for convergence history.
The GMRES method.
virtual void smoother_setup(DoubleMatrixBase *matrix_pt)=0
Set up the smoother for the matrix specified by the pointer.
void disable_iteration_restart()
switches off iteration restart
MATRIX * Matrix_pt
Pointer to matrix.
void enable_setup_preconditioner_before_solve()
Setup the preconditioner before the solve.
The conjugate gradient method.
void disable_doc_convergence_history()
Disable documentation of the convergence history.
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 solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
GS()
Constructor.
double preconditioner_setup_time() const
Returns the time taken to set up the preconditioner.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void operator=(const BiCGStab &)
Broken assignment operator.
bool Use_as_smoother
When a derived class object is being used as a smoother in the MG solver (or elsewhere) the residual ...
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 enable_doc_convergence_history()
Enable documentation of the convergence history.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
void solve(Problem *const &problem_pt, DoubleVector &result)
Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the r...
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...
double linear_solver_solution_time() const
return the time taken to solve the linear system
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
unsigned Restart
The number of iterations before the iteration proceedure is restarted if iteration restart is used...
bool Matrix_can_be_deleted
Boolean flag to indicate if the matrix pointed to be Matrix_pt can be deleted.
void smoother_solve(const DoubleVector &rhs, DoubleVector &solution)
The smoother_solve function performs fixed number of iterations on the system A*result=rhs. The number of (smoothing) iterations is the same as the max. number of iterations in the underlying IterativeLinearSolver class.
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 ...
GS(const GS &)
Broken copy constructor.
void extract_diagonal_entries(DoubleMatrixBase *matrix_pt)
Function to extract the diagonal entries from the matrix.
double Jacobian_setup_time
Jacobian setup time.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
double jacobian_setup_time() const
returns the time taken to assemble the jacobian matrix and residual vector
unsigned iterations() const
Number of iterations taken.
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
static char t char * s
Definition: cfortran.h:572
void operator=(const DampedJacobi &)
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 ...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Re-solve the system defined by the last assembled Jacobian and the rhs vector specified here...
Smoother()
Empty constructor.
MATRIX * Matrix_pt
Pointer to matrix.
unsigned Iterations
Number of iterations taken.
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 disable_error_after_max_iter()
Don't throw an error if we don't converge within max_iter (default).
MATRIX * Matrix_pt
System matrix pointer in the format specified by the template argument.
void solve_helper(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &solution)
General interface to solve function.
double Omega
Damping factor.
CG()
Constructor.
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...
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
void enable_iteration_restart(const unsigned &restart)
switches on iteration restarting and takes as an argument the number of iterations after which the co...
unsigned Iterations
Number of iterations taken.
CG(const CG &)
Broken copy constructor.
double Solution_time
linear solver solution time
virtual ~GMRES()
Destructor (cleanup storage)
IterativeLinearSolver(const IterativeLinearSolver &)
Broken copy constructor.
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
void clean_up_memory()
Cleanup data that's stored for resolve (if any has been stored)
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
DampedJacobi(const double &omega=2.0/3.0)
Empty constructor.
DampedJacobi(const DampedJacobi &)
Broken copy constructor.
void operator=(const IterativeLinearSolver &)
Broken assignment operator.
The conjugate gradient method.
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 disable_setup_preconditioner_before_solve()
Don't set up the preconditioner before the solve.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void set_preconditioner_LHS()
Set left preconditioning (the default)
The Gauss Seidel method.
virtual ~Smoother()
Virtual empty destructor.
unsigned iterations() const
Number of iterations taken.
bool Doc_convergence_history
Flag indicating if the convergence history is to be documented.
unsigned Iterations
Number of iterations taken.
void close_convergence_history_file_stream()
Close convergence history output stream.
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
void smoother_setup(DoubleMatrixBase *matrix_pt)
Set up the smoother for the matrix specified by the pointer.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned & max_iter()
Access to max. number of iterations.
bool Use_iterative_solver_as_preconditioner
Use the iterative solver as preconditioner.
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...
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...
double preconditioner_setup_time() const
Returns the time taken to set up the preconditioner.
void operator=(const GS &)
Broken assignment operator.
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...
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 ...
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...
Vector< double > Matrix_diagonal
Vector containing the diagonal entries of A.
bool iteration_restart() const
access function indicating whether restarted GMRES 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. The number of (smoothing) iterations is the same as the max. number of iterations in the underlying IterativeLinearSolver class.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
virtual ~GS()
Destructor (cleanup storage)
GMRES(const GMRES &)
Broken copy constructor.
virtual ~CG()
Destructor (cleanup storage)
BiCGStab(const BiCGStab &)
Broken copy constructor.
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.
bool Setup_preconditioner_before_solve
indicates whether the preconditioner should be setup before solve. Default = true; ...
void set_preconditioner_RHS()
Enable right preconditioning.
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 operator=(const GMRES &)
Broken assignment operator.
Abstract base class for matrices of doubles – adds abstract interfaces for solving, LU decomposition and multiplication by vectors.
Definition: matrices.h:275
unsigned Iterations
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...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
double Preconditioner_setup_time
Preconditioner setup time.
MATRIX * Matrix_pt
Pointer to matrix.
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
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 ...
bool Resolving
Boolean flag to indicate if the solve is done in re-solve mode, bypassing setup of matrix and precond...