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-2022 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
46namespace 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 {
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>
284 {
285 public:
286 /// Constructor
288 : Iterations(0),
289 Matrix_pt(0),
290 Resolving(false),
292 {
293 }
294
295
296 /// Destructor (cleanup storage)
297 virtual ~CG()
298 {
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 {
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
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
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 {
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 {
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
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
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
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
603 : Matrix_pt(0),
604 Iterations(0),
605 Resolving(false),
607 {
608 }
609
610 /// Destructor (cleanup storage)
611 virtual ~GS()
612 {
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 {
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
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
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
789 : Matrix_pt(0),
790 Iterations(0),
791 Resolving(false),
793 {
794 }
795
796 /// Destructor (cleanup storage)
797 virtual ~GS()
798 {
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 {
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
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
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
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)
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 {
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
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),
1236 {
1237 Preconditioner_LHS = true;
1238 Iteration_restart = false;
1239 }
1240
1241 /// Destructor (cleanup storage)
1242 virtual ~GMRES()
1243 {
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 {
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
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
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
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),
1565 {
1566 Preconditioner_LHS = true;
1567 Iteration_restart = false;
1568 }
1569
1570 /// Destructor: Clean up storage
1572 {
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 {
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
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
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 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 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 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...
unsigned nrow() const
access function to the number of global rows.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
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 dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
double * values_pt()
access function to the underlying values
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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.
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 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 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.
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.
unsigned & max_iter()
Access to max. number of iterations.
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.
double & tolerance()
Access to convergence tolerance.
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 Solution_time
linear solver solution time
void operator=(const IterativeLinearSolver &)=delete
Broken assignment operator.
IterativeLinearSolver(const IterativeLinearSolver &)=delete
Broken copy constructor.
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
std::ofstream Output_file_stream
Output file stream for convergence history.
Preconditioner *const & preconditioner_pt() const
Access function to preconditioner (const version)
void disable_setup_preconditioner_before_solve()
Don't set up the preconditioner before the solve.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...