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-2024 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 
27 // This header defines a class for linear solvers
28 
29 // Include guards
30 #ifndef OOMPH_LINEAR_SOLVER_HEADER
31 #define OOMPH_LINEAR_SOLVER_HEADER
32 
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36 #include <oomph-lib-config.h>
37 #endif
38 
39 #ifdef OOMPH_HAS_MPI
40 #include "mpi.h"
41 #endif
42 
43 // oomph-lib headers
44 #include "Vector.h"
45 #include "double_vector.h"
46 #include "matrices.h"
47 
48 namespace oomph
49 {
50  // Forward declaration of problem class
51  class Problem;
52 
53  //====================================================================
54  /// Base class for all linear solvers. This merely defines standard
55  /// interfaces for linear solvers, so that different solvers can be
56  /// used in a clean and transparent manner. Note that LinearSolvers
57  /// are primarily used to solve the linear systems arising in
58  /// oomph-lib's Newton iteration. Their primary solve function
59  /// therefore takes a pointer to the associated problem, construct its
60  /// Jacobian matrix and residual vector, and return the solution
61  /// of the linear system formed by the Jacobian and the residual vector.
62  /// We also provide broken virtual interfaces
63  /// to a linear-algebra-type solve function in which the matrix
64  /// and the rhs can be specified, but this are not guaranteed to
65  /// implemented for all linear solvers (e.g. for frontal solvers).
66  //====================================================================
68  {
69  protected:
70  /// Boolean that indicates whether the matrix (or its factors, in
71  /// the case of direct solver) should be stored so that the resolve function
72  /// can be used.
74 
75  /// Boolean flag that indicates whether the time taken
76  // for the solves should be documented
77  bool Doc_time;
78 
79  /// flag that indicates whether the gradient required for the
80  /// globally convergent Newton method should be computed or not
82 
83  /// flag that indicates whether the gradient was computed or not
85 
86  /// DoubleVector storing the gradient for the globally convergent
87  /// Newton method
89 
90  public:
91  /// Empty constructor, initialise the member data
93  : Enable_resolve(false),
94  Doc_time(true),
95  Compute_gradient(false),
97  {
98  }
99 
100  /// Broken copy constructor
101  LinearSolver(const LinearSolver& dummy) = delete;
102 
103  /// Broken assignment operator
104  void operator=(const LinearSolver&) = delete;
105 
106  /// Empty virtual destructor
107  virtual ~LinearSolver() {}
108 
109  /// Enable documentation of solve times
111  {
112  Doc_time = true;
113  }
114 
115  /// Disable documentation of solve times
117  {
118  Doc_time = false;
119  }
120 
121  /// Is documentation of solve times enabled?
122  bool is_doc_time_enabled() const
123  {
124  return Doc_time;
125  }
126 
127  /// Boolean flag indicating if resolves are enabled
128  bool is_resolve_enabled() const
129  {
130  return Enable_resolve;
131  }
132 
133  /// Enable resolve (i.e. store matrix and/or LU decomposition, say)
134  /// Virtual so it can be overloaded to perform additional tasks
135  virtual void enable_resolve()
136  {
137  Enable_resolve = true;
138  }
139 
140  /// Disable resolve (i.e. store matrix and/or LU decomposition, say)
141  /// This function simply resets an internal flag. It's virtual so
142  /// it can be overloaded to perform additional tasks such as
143  /// cleaning up memory that is only required for the resolve.
144  virtual void disable_resolve()
145  {
146  Enable_resolve = false;
147  }
148 
149  /// Solver: Takes pointer to problem and returns the results vector
150  /// which contains the solution of the linear system defined by
151  /// the problem's fully assembled Jacobian and residual vector.
152  virtual void solve(Problem* const& problem_pt, DoubleVector& result) = 0;
153 
154  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
155  /// vector and returns the solution of the linear system.
156  virtual void solve(DoubleMatrixBase* const& matrix_pt,
157  const DoubleVector& rhs,
158  DoubleVector& result)
159  {
160  throw OomphLibError(
161  "DoubleVector based solve function not implemented for this solver",
162  OOMPH_CURRENT_FUNCTION,
163  OOMPH_EXCEPTION_LOCATION);
164  }
165 
166  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
167  /// vector and returns the solution of the linear system.
168  virtual void solve(DoubleMatrixBase* const& matrix_pt,
169  const Vector<double>& rhs,
170  Vector<double>& result)
171  {
172  throw OomphLibError(
173  "Vector<double> based solve function not implemented for this solver",
174  OOMPH_CURRENT_FUNCTION,
175  OOMPH_EXCEPTION_LOCATION);
176  }
177 
178  /// Solver: Takes pointer to problem and returns the results vector
179  /// which contains the solution of the linear system defined by the
180  /// problem's fully assembled Jacobian and residual vector (broken virtual).
181  virtual void solve_transpose(Problem* const& problem_pt,
182  DoubleVector& result)
183  {
184  // Create an output stream
185  std::ostringstream error_message_stream;
186 
187  // Create the error message
188  error_message_stream << "The function to solve the transposed system has "
189  << "not yet been\nimplemented for this solver."
190  << std::endl;
191 
192  // Throw the error message
193  throw OomphLibError(error_message_stream.str(),
194  OOMPH_CURRENT_FUNCTION,
195  OOMPH_EXCEPTION_LOCATION);
196  } // End of solve_transpose
197 
198  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
199  /// vector and returns the solution of the linear system.
200  virtual void solve_transpose(DoubleMatrixBase* const& matrix_pt,
201  const DoubleVector& rhs,
202  DoubleVector& result)
203  {
204  throw OomphLibError(
205  "DoubleVector based solve function not implemented for this solver",
206  OOMPH_CURRENT_FUNCTION,
207  OOMPH_EXCEPTION_LOCATION);
208  }
209 
210  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
211  /// vector and returns the solution of the linear system.
212  virtual void solve_transpose(DoubleMatrixBase* const& matrix_pt,
213  const Vector<double>& rhs,
214  Vector<double>& result)
215  {
216  throw OomphLibError(
217  "Vector<double> based solve function not implemented for this solver",
218  OOMPH_CURRENT_FUNCTION,
219  OOMPH_EXCEPTION_LOCATION);
220  }
221 
222  /// Resolve the system defined by the last assembled jacobian
223  /// and the rhs vector. Solution is returned in the vector result.
224  /// (broken virtual)
225  virtual void resolve(const DoubleVector& rhs, DoubleVector& result)
226  {
227  throw OomphLibError(
228  "Resolve function not implemented for this linear solver",
229  OOMPH_CURRENT_FUNCTION,
230  OOMPH_EXCEPTION_LOCATION);
231  }
232 
233  /// Solver: Resolve the system defined by the last assembled jacobian
234  /// and the rhs vector. Solution is returned in the vector result.
235  /// (broken virtual)
236  virtual void resolve_transpose(const DoubleVector& rhs,
237  DoubleVector& result)
238  {
239  // Create an output stream
240  std::ostringstream error_message_stream;
241 
242  // Create the error message
243  error_message_stream
244  << "The function to resolve the transposed system has "
245  << "not yet been\nimplemented for this solver." << std::endl;
246 
247  // Throw the error message
248  throw OomphLibError(error_message_stream.str(),
249  OOMPH_CURRENT_FUNCTION,
250  OOMPH_EXCEPTION_LOCATION);
251  } // End of resolve_transpose
252 
253  /// Empty virtual function that can be overloaded in specific
254  /// linear solvers to clean up any memory that may have been
255  /// allocated (e.g. when preparing for a re-solve).
256  virtual void clean_up_memory() {}
257 
258  /// returns the time taken to assemble the Jacobian matrix and
259  /// residual vector (needs to be overloaded for each solver)
260  virtual double jacobian_setup_time() const
261  {
262  throw OomphLibError(
263  "jacobian_setup_time has not been implemented for this linear solver",
264  OOMPH_CURRENT_FUNCTION,
265  OOMPH_EXCEPTION_LOCATION);
266  return 0;
267  }
268 
269  /// return the time taken to solve the linear system (needs to be
270  /// overloaded for each linear solver)
271  virtual double linear_solver_solution_time() const
272  {
273  throw OomphLibError(
274  "linear_solver_solution_time() not implemented for this linear solver",
275  OOMPH_CURRENT_FUNCTION,
276  OOMPH_EXCEPTION_LOCATION);
277  return 0;
278  }
279 
280  /// function to enable the computation of the gradient required
281  /// for the globally convergent Newton method
283  {
284  throw OomphLibError(
285  "enable_computation_of_gradient() not implemented for "
286  "this linear solver",
287  OOMPH_CURRENT_FUNCTION,
288  OOMPH_EXCEPTION_LOCATION);
289  }
290 
291  /// function to disable the computation of the gradient required
292  /// for the globally convergent Newton method
294  {
295  Compute_gradient = false;
296  }
297 
298  /// function to reset the size of the gradient before each Newton
299  /// solve
301  {
303  }
304 
305  /// function to access the gradient, provided it has been computed
306  void get_gradient(DoubleVector& gradient)
307  {
308 #ifdef PARANOID
310  {
311 #endif
313 #ifdef PARANOID
314  }
315  else
316  {
317  throw OomphLibError(
318  "The gradient has not been computed for this linear solver!",
319  OOMPH_CURRENT_FUNCTION,
320  OOMPH_EXCEPTION_LOCATION);
321  }
322 #endif
323  }
324  };
325 
326  //=============================================================================
327  /// Dense LU decomposition-based solve of full assembled linear system.
328  /// VERY inefficient but useful to illustrate the principle.
329  /// Only suitable for use with Serial matrices and vectors.
330  /// This solver will only work with non-distributed matrices and vectors
331  /// (note: DenseDoubleMatrix is not distributable)
332  //============================================================================
333  class DenseLU : public LinearSolver
334  {
335  /// The DenseDoubleMatrix class is a friend
336  friend class DenseDoubleMatrix;
337 
338  public:
339  /// Constructor, initialise storage
341  : Jacobian_setup_time(0.0),
342  Solution_time(0.0),
344  Index(0),
345  LU_factors(0)
346  {
347  // Shut up!
348  Doc_time = false;
349  }
350 
351  /// Broken copy constructor
352  DenseLU(const DenseLU& dummy) = delete;
353 
354  /// Broken assignment operator
355  void operator=(const DenseLU&) = delete;
356 
357  /// Destructor, clean up the stored LU factors
359  {
360  clean_up_memory();
361  }
362 
363  /// Solver: Takes pointer to problem and returns the results Vector
364  /// which contains the solution of the linear system defined by
365  /// the problem's fully assembled Jacobian and residual Vector.
366  void solve(Problem* const& problem_pt, DoubleVector& result);
367 
368  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
369  /// vector and returns the solution of the linear system.
370  void solve(DoubleMatrixBase* const& matrix_pt,
371  const DoubleVector& rhs,
372  DoubleVector& result);
373 
374  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
375  /// vector and returns the solution of the linear system.
376  void solve(DoubleMatrixBase* const& matrix_pt,
377  const Vector<double>& rhs,
378  Vector<double>& result);
379 
380  /// returns the time taken to assemble the jacobian matrix and
381  /// residual vector
382  double jacobian_setup_time() const
383  {
384  return Jacobian_setup_time;
385  }
386 
387  /// return the time taken to solve the linear system (needs to be
388  /// overloaded for each linear solver)
389  virtual double linear_solver_solution_time() const
390  {
391  return Solution_time;
392  }
393 
394  protected:
395  /// Perform the LU decomposition of the matrix
396  void factorise(DoubleMatrixBase* const& matrix_pt);
397 
398  /// Do the backsubstitution step to solve the system LU result = rhs
399  void backsub(const DoubleVector& rhs, DoubleVector& result);
400 
401  /// perform back substitution using Vector<double>
402  void backsub(const Vector<double>& rhs, Vector<double>& result);
403 
404  /// Clean up the stored LU factors
405  void clean_up_memory();
406 
407  /// Jacobian setup time
409 
410  /// Solution time
412 
413  /// Sign of the determinant of the matrix
414  /// (obtained during the LU decomposition)
416 
417  private:
418  /// Pointer to storage for the index of permutations in the LU solve
419  long* Index;
420 
421  /// Pointer to storage for the LU decomposition
422  double* LU_factors;
423  };
424 
425 
426  //====================================================================
427  /// Dense LU decomposition-based solve of linear system
428  /// assembled via finite differencing of the residuals Vector.
429  /// Even more inefficient than DenseLU but excellent sanity check!
430  //====================================================================
431  class FD_LU : public DenseLU
432  {
433  public:
434  /// Constructor: empty
435  FD_LU() : DenseLU() {}
436 
437  /// Broken copy constructor
438  FD_LU(const FD_LU& dummy) = delete;
439 
440  /// Broken assignment operator
441  void operator=(const FD_LU&) = delete;
442 
443  /// Solver: Takes pointer to problem and returns the results Vector
444  /// which contains the solution of the linear system defined by
445  /// the problem's residual Vector (Jacobian computed by FD approx.)
446  void solve(Problem* const& problem_pt, DoubleVector& result);
447 
448  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
449  /// vector and returns the solution of the linear system.
450  void solve(DoubleMatrixBase* const& matrix_pt,
451  const DoubleVector& rhs,
452  DoubleVector& result)
453  {
454  DenseLU::solve(matrix_pt, rhs, result);
455  }
456 
457  /// Linear-algebra-type solver: Takes pointer to a matrix
458  /// and rhs vector and returns the solution of the linear system
459  /// Call the broken base-class version. If you want this, please
460  /// implement it
461  void solve(DoubleMatrixBase* const& matrix_pt,
462  const Vector<double>& rhs,
463  Vector<double>& result)
464  {
465  LinearSolver::solve(matrix_pt, rhs, result);
466  }
467  };
468 
469 
470  /// ////////////////////////////////////////////////////////////////////////////
471  /// ////////////////////////////////////////////////////////////////////////////
472  /// ////////////////////////////////////////////////////////////////////////////
473 
474 
475  //=============================================================================
476  /// SuperLU Project Solver class. This is a combined wrapper for both
477  /// SuperLU and SuperLU Dist.
478  /// See http://crd.lbl.gov/~xiaoye/SuperLU/
479  /// Default Behaviour: If this solver is distributed over more than one
480  /// processor then SuperLU Dist is used.
481  /// Member data naming convention: member data associated with the SuperLU
482  /// Dist solver begins Dist_... and member data associated with the serial
483  /// SuperLU solver begins Serial_... .
484  //=============================================================================
486  {
487  public:
488  /// enum type to specify the solver behaviour.
489  /// Default - will employ superlu dist if more than 1 processor.
490  /// Serial - will always try use superlu (serial).
491  /// Distributed - will always try to use superlu dist.
492  enum Type
493  {
497  };
498 
499  /// Constructor. Set the defaults.
501  {
502  // Set solver wide default values and null pointers
503  Doc_stats = false;
504  Suppress_solve = false;
505  Using_dist = false;
507 
508 #ifdef OOMPH_HAS_MPI
509  // Set default values and nullify pointers for SuperLU Dist
510  Dist_use_global_solver = false;
511  Dist_delete_matrix_data = false;
516  Dist_value_pt = 0;
517  Dist_index_pt = 0;
518  Dist_start_pt = 0;
519 #endif
520 
521  // Set default values and null pointers for SuperLU (serial)
524  Serial_n_dof = 0;
525  }
526 
527  /// Broken copy constructor
528  SuperLUSolver(const SuperLUSolver& dummy) = delete;
529 
530  /// Broken assignment operator
531  void operator=(const SuperLUSolver&) = delete;
532 
533  /// Destructor, clean up the stored matrices
535  {
536  clean_up_memory();
537  }
538 
539  /// function to enable the computation of the gradient
541  {
542  Compute_gradient = true;
543  }
544 
545  // SuperLUSolver methods
546  /// /////////////////////
547 
548  /// Overload disable resolve so that it cleans up memory too
550  {
552  clean_up_memory();
553  }
554 
555  /// Solver: Takes pointer to problem and returns the results Vector
556  /// which contains the solution of the linear system defined by
557  /// the problem's fully assembled Jacobian and residual Vector.
558  void solve(Problem* const& problem_pt, DoubleVector& result);
559 
560  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
561  /// vector and returns the solution of the linear system.
562  /// The function returns the global result Vector.
563  /// Note: if Delete_matrix_data is true the function
564  /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
565  void solve(DoubleMatrixBase* const& matrix_pt,
566  const DoubleVector& rhs,
567  DoubleVector& result);
568 
569 
570  /* /// Linear-algebra-type solver: Takes pointer to a matrix */
571  /* /// and rhs vector and returns the solution of the linear system */
572  /* /// Call the broken base-class version. If you want this, please */
573  /* /// implement it */
574  /* void solve(DoubleMatrixBase* const &matrix_pt, */
575  /* const Vector<double> &rhs, */
576  /* Vector<double> &result) */
577  /* {LinearSolver::solve(matrix_pt,rhs,result);} */
578 
579  /// Solver: Takes pointer to problem and returns the results Vector
580  /// which contains the solution of the linear system defined by
581  /// the problem's fully assembled Jacobian and residual Vector.
582  void solve_transpose(Problem* const& problem_pt, DoubleVector& result);
583 
584  /// Linear-algebra-type solver: Takes pointer to a matrix and rhs
585  /// vector and returns the solution of the linear system.
586  /// The function returns the global result Vector.
587  /// Note: if Delete_matrix_data is true the function
588  /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
589  void solve_transpose(DoubleMatrixBase* const& matrix_pt,
590  const DoubleVector& rhs,
591  DoubleVector& result);
592 
593  /// Resolve the system defined by the last assembled jacobian
594  /// and the specified rhs vector if resolve has been enabled.
595  /// Note: returns the global result Vector.
596  void resolve(const DoubleVector& rhs, DoubleVector& result);
597 
598  /// Resolve the (transposed) system defined by the last assembled
599  /// Jacobian and the specified rhs vector if resolve has been enabled.
600  void resolve_transpose(const DoubleVector& rhs, DoubleVector& result);
601 
602  /// Enable documentation of solver statistics
604  {
605  Doc_stats = true;
606  }
607 
608  /// Disable documentation of solver statistics
610  {
611  Doc_stats = false;
612  }
613 
614  /// returns the time taken to assemble the jacobian matrix and
615  /// residual vector
616  double jacobian_setup_time() const
617  {
618  return Jacobian_setup_time;
619  }
620 
621  /// return the time taken to solve the linear system (needs to be
622  /// overloaded for each linear solver)
623  virtual double linear_solver_solution_time() const
624  {
625  return Solution_time;
626  }
627 
628  /// Do the factorisation stage
629  /// Note: if Delete_matrix_data is true the function
630  /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
631  void factorise(DoubleMatrixBase* const& matrix_pt);
632 
633  /// Do the backsubstitution for SuperLU solver
634  /// Note: returns the global result Vector.
635  void backsub(const DoubleVector& rhs, DoubleVector& result);
636 
637  /// Do the back-substitution for transposed system of the SuperLU
638  /// solver Note: Returns the global result Vector.
639  void backsub_transpose(const DoubleVector& rhs, DoubleVector& result);
640 
641  /// Clean up the memory allocated by the solver
642  void clean_up_memory();
643 
644  /// Specify the solve type. Either default, serial or distributed.
645  /// See enum SuperLU_solver_type for more details.
646  void set_solver_type(const Type& t)
647  {
648  this->clean_up_memory();
649  Solver_type = t;
650  }
651 
652  // SuperLU (serial) methods
653  /// ////////////////////////
654 
655  /// Use the compressed row format in superlu serial
657  {
659  }
660 
661  /// Use the compressed column format in superlu serial
663  {
665  }
666 
667 #ifdef OOMPH_HAS_MPI
668 
669  // SuperLU Dist methods
670  /// ////////////////////
671 
672  /// Set Delete_matrix_data flag. SuperLU_dist needs its own copy
673  /// of the input matrix, therefore a copy must be made if any matrix
674  /// used with this solver is to be preserved. If the input matrix can be
675  /// deleted the flag can be set to true to reduce the amount of memory
676  /// required, and the matrix data will be wiped using its clean_up_memory()
677  /// function. Default value is false.
679  {
681  }
682 
683  /// Unset Delete_matrix_data flag. SuperLU_dist needs its own copy
684  /// of the input matrix, therefore a copy must be made if any matrix
685  /// used with this solver is to be preserved. If the input matrix can be
687  {
688  Dist_delete_matrix_data = false;
689  }
690 
691  /// Set flag so that SuperLU_DIST is allowed to permute matrix rows
692  /// and columns during factorisation. This is the default for SuperLU_DIST,
693  /// and can lead to significantly faster solves.
695  {
697  }
698 
699  /// Set flag so that SuperLU_DIST is not allowed to permute matrix
700  /// rows and columns during factorisation.
702  {
704  }
705 
706  /// Calling this method will ensure that when the problem based
707  /// solve interface is used, a global (serial) jacobian will be
708  /// assembled.
709  /// Note: calling this function will delete any distributed solve data.
711  {
713  {
714  clean_up_memory();
715  Dist_use_global_solver = true;
716  }
717  }
718 
719  /// Calling this method will ensure that when the problem based
720  /// solve interface is used, a distributed jacobian will be
721  /// assembled.
722  /// Note: calling this function will delete any global solve data.
724  {
726  {
727  clean_up_memory();
728  Dist_use_global_solver = false;
729  }
730  }
731 
732 #endif
733 
734  private:
735  /// factorise method for SuperLU (serial)
736  void factorise_serial(DoubleMatrixBase* const& matrix_pt);
737 
738  /// backsub method for SuperLU (serial)
739  void backsub_serial(const DoubleVector& rhs, DoubleVector& result);
740 
741  /// backsub method for SuperLU (serial)
742  void backsub_transpose_serial(const DoubleVector& rhs,
743  DoubleVector& result);
744 
745 #ifdef OOMPH_HAS_MPI
746  /// factorise method for SuperLU Dist
747  void factorise_distributed(DoubleMatrixBase* const& matrix_pt);
748 
749  /// backsub method for SuperLU Dist
750  void backsub_distributed(const DoubleVector& rhs, DoubleVector& result);
751 
752  /// backsub method for SuperLU Dist
754  DoubleVector& result);
755 #endif
756 
757  // SuperLUSolver member data
758  /// /////////////////////////
759 
760  /// Jacobian setup time
762 
763  /// Solution time
765 
766  /// Suppress solve?
768 
769  /// Set to true to output statistics (false by default).
770  bool Doc_stats;
771 
772  /// the solver type. see SuperLU_solver_type for details.
774 
775  /// boolean flag indicating whether superlu dist is being used
777 
778  // SuperLU (serial) member data
779  /// ////////////////////////////
780 
781  /// Storage for the LU factors as required by SuperLU
783 
784  /// Info flag for the SuperLU solver
786 
787  /// The number of unknowns in the linear system
788  unsigned long Serial_n_dof;
789 
790  /// Sign of the determinant of the matrix
792 
793  /// Use compressed row version?
795 
796  public:
797  /// How much memory do the LU factors take up? In bytes
799 
800  /// How much memory was allocated by SuperLU? In bytes
801  double get_total_needed_memory();
802 
803  private:
804 #ifdef OOMPH_HAS_MPI
805 
806  // SuperLU Dist member data
807  /// ////////////////////////
808  public:
809  /// Static flag that determines whether the warning about
810  /// incorrect distribution of RHSs will be printed or not
812 
813  private:
814  /// Flag that determines whether the MPIProblem based solve function
815  /// uses the global or distributed version of SuperLU_DIST
816  /// (default value is false).
818 
819  /// Flag is true if solve data has been generated for a global matrix
821 
822  /// Flag is true if solve data has been generated for distributed matrix
824 
825  /// Storage for the LU factors and other data required by SuperLU
827 
828  /// Number of rows for the process grid
830 
831  /// Number of columns for the process grid
833 
834  /// Info flag for the SuperLU solver
836 
837  /// If true then SuperLU_DIST is allowed to permute matrix rows
838  /// and columns during factorisation. This is the default for SuperLU_DIST,
839  /// and can lead to significantly faster solves, but has been known to
840  /// fail, hence the default value is 0.
842 
843  /// Delete_matrix_data flag. SuperLU_dist needs its own copy
844  /// of the input matrix, therefore a copy must be made if any matrix
845  /// used with this solver is to be preserved. If the input matrix can be
846  /// deleted the flag can be set to true to reduce the amount of memory
847  /// required, and the matrix data will be wiped using its clean_up_memory()
848  /// function. Default value is false.
850 
851  /// Pointer for storage of the matrix values required by SuperLU_DIST
852  double* Dist_value_pt;
853 
854  /// Pointer for storage of matrix rows or column indices required
855  /// by SuperLU_DIST
857 
858  /// Pointers for storage of matrix column or row starts
859  // required by SuperLU_DIST
861 
862 #endif
863  }; // end of SuperLUSolver
864 } // namespace oomph
865 #endif
char t
Definition: cfortran.h:568
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
Dense LU decomposition-based solve of full assembled linear system. VERY inefficient but useful to il...
~DenseLU()
Destructor, clean up the stored LU factors.
double Jacobian_setup_time
Jacobian setup time.
double jacobian_setup_time() const
returns the time taken to assemble the jacobian matrix and residual vector
void backsub(const DoubleVector &rhs, DoubleVector &result)
Do the backsubstitution step to solve the system LU result = rhs.
void clean_up_memory()
Clean up the stored LU factors.
long * Index
Pointer to storage for the index of permutations in the LU solve.
DenseLU()
Constructor, initialise storage.
void factorise(DoubleMatrixBase *const &matrix_pt)
Perform the LU decomposition of the matrix.
double * LU_factors
Pointer to storage for the LU decomposition.
void operator=(const DenseLU &)=delete
Broken assignment operator.
virtual double linear_solver_solution_time() const
return the time taken to solve the linear system (needs to be overloaded for each linear solver)
double Solution_time
Solution time.
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...
DenseLU(const DenseLU &dummy)=delete
Broken copy constructor.
int Sign_of_determinant_of_matrix
Sign of the determinant of the matrix (obtained during the LU decomposition)
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
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
void clear()
wipes the DoubleVector
Dense LU decomposition-based solve of linear system assembled via finite differencing of the residual...
FD_LU()
Constructor: empty.
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 solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results Vector which contains the solution of the li...
FD_LU(const FD_LU &dummy)=delete
Broken copy constructor.
void operator=(const FD_LU &)=delete
Broken assignment operator.
void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
Base class for all linear solvers. This merely defines standard interfaces for linear solvers,...
Definition: linear_solver.h:68
void operator=(const LinearSolver &)=delete
Broken assignment operator.
virtual double jacobian_setup_time() const
returns the time taken to assemble the Jacobian matrix and residual vector (needs to be overloaded fo...
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...
virtual 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 ...
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
virtual void enable_computation_of_gradient()
function to enable the computation of the gradient required for the globally convergent Newton method
virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the rhs vector. Solution is returned in...
virtual double linear_solver_solution_time() const
return the time taken to solve the linear system (needs to be overloaded for each linear solver)
bool is_doc_time_enabled() const
Is documentation of solve times enabled?
virtual void solve_transpose(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 Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:77
virtual void solve(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
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 solve_transpose(DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
void get_gradient(DoubleVector &gradient)
function to access the gradient, provided it has been computed
LinearSolver(const LinearSolver &dummy)=delete
Broken copy constructor.
LinearSolver()
Empty constructor, initialise the member data.
Definition: linear_solver.h:92
void disable_computation_of_gradient()
function to disable the computation of the gradient required for the globally convergent Newton metho...
virtual void solve_transpose(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
void disable_doc_time()
Disable documentation of solve times.
virtual ~LinearSolver()
Empty virtual destructor.
virtual void clean_up_memory()
Empty virtual function that can be overloaded in specific linear solvers to clean up any memory that ...
bool Gradient_has_been_computed
flag that indicates whether the gradient was computed or not
Definition: linear_solver.h:84
void enable_doc_time()
Enable documentation of solve times.
void reset_gradient()
function to reset the size of the gradient before each Newton solve
virtual void resolve_transpose(const DoubleVector &rhs, DoubleVector &result)
Solver: Resolve the system defined by the last assembled jacobian and the rhs vector....
DoubleVector Gradient_for_glob_conv_newton_solve
DoubleVector storing the gradient for the globally convergent Newton method.
Definition: linear_solver.h:88
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
bool is_resolve_enabled() const
Boolean flag indicating if resolves are enabled.
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....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:154
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
SuperLUSolver(const SuperLUSolver &dummy)=delete
Broken copy constructor.
bool Doc_stats
Set to true to output statistics (false by default).
bool Dist_delete_matrix_data
Delete_matrix_data flag. SuperLU_dist needs its own copy of the input matrix, therefore a copy must b...
void solve_transpose(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results Vector which contains the solution of the li...
static bool Suppress_incorrect_rhs_distribution_warning_in_resolve
Static flag that determines whether the warning about incorrect distribution of RHSs will be printed ...
int * Dist_index_pt
Pointer for storage of matrix rows or column indices required by SuperLU_DIST.
Type Solver_type
the solver type. see SuperLU_solver_type for details.
void operator=(const SuperLUSolver &)=delete
Broken assignment operator.
void disable_row_and_col_permutations_in_superlu_dist()
Set flag so that SuperLU_DIST is not allowed to permute matrix rows and columns during factorisation.
virtual double linear_solver_solution_time() const
return the time taken to solve the linear system (needs to be overloaded for each linear solver)
void backsub_serial(const DoubleVector &rhs, DoubleVector &result)
backsub method for SuperLU (serial)
void enable_doc_stats()
Enable documentation of solver statistics.
void use_distributed_solve_in_superlu_dist()
Calling this method will ensure that when the problem based solve interface is used,...
void backsub_transpose_distributed(const DoubleVector &rhs, DoubleVector &result)
backsub method for SuperLU Dist
double Solution_time
Solution time.
double * Dist_value_pt
Pointer for storage of the matrix values required by SuperLU_DIST.
void use_global_solve_in_superlu_dist()
Calling this method will ensure that when the problem based solve interface is used,...
bool Serial_compressed_row_flag
Use compressed row version?
Type
enum type to specify the solver behaviour. Default - will employ superlu dist if more than 1 processo...
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 factorise(DoubleMatrixBase *const &matrix_pt)
Do the factorisation stage Note: if Delete_matrix_data is true the function matrix_pt->clean_up_memor...
void factorise_serial(DoubleMatrixBase *const &matrix_pt)
factorise method for SuperLU (serial)
void backsub_distributed(const DoubleVector &rhs, DoubleVector &result)
backsub method for SuperLU Dist
bool Using_dist
boolean flag indicating whether superlu dist is being used
int * Dist_start_pt
Pointers for storage of matrix column or row starts.
unsigned long Serial_n_dof
The number of unknowns in the linear system.
void set_solver_type(const Type &t)
Specify the solve type. Either default, serial or distributed. See enum SuperLU_solver_type for more ...
void * Serial_f_factors
Storage for the LU factors as required by SuperLU.
void resolve_transpose(const DoubleVector &rhs, DoubleVector &result)
Resolve the (transposed) system defined by the last assembled Jacobian and the specified rhs vector i...
SuperLUSolver()
Constructor. Set the defaults.
void disable_delete_matrix_data_in_superlu_dist()
Unset Delete_matrix_data flag. SuperLU_dist needs its own copy of the input matrix,...
bool Suppress_solve
Suppress solve?
void enable_row_and_col_permutations_in_superlu_dist()
Set flag so that SuperLU_DIST is allowed to permute matrix rows and columns during factorisation....
void enable_delete_matrix_data_in_superlu_dist()
Set Delete_matrix_data flag. SuperLU_dist needs its own copy of the input matrix, therefore a copy mu...
~SuperLUSolver()
Destructor, clean up the stored matrices.
double jacobian_setup_time() const
returns the time taken to assemble the jacobian matrix and residual vector
int Serial_sign_of_determinant_of_matrix
Sign of the determinant of the matrix.
bool Dist_use_global_solver
Flag that determines whether the MPIProblem based solve function uses the global or distributed versi...
double get_memory_usage_for_lu_factors()
How much memory do the LU factors take up? In bytes.
void factorise_distributed(DoubleMatrixBase *const &matrix_pt)
factorise method for SuperLU Dist
bool Dist_allow_row_and_col_permutations
If true then SuperLU_DIST is allowed to permute matrix rows and columns during factorisation....
int Dist_info
Info flag for the SuperLU solver.
void use_compressed_column_for_superlu_serial()
Use the compressed column format in superlu serial.
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
bool Dist_global_solve_data_allocated
Flag is true if solve data has been generated for a global matrix.
double get_total_needed_memory()
How much memory was allocated by SuperLU? In bytes.
double Jacobian_setup_time
Jacobian setup time.
void backsub_transpose_serial(const DoubleVector &rhs, DoubleVector &result)
backsub method for SuperLU (serial)
void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the specified rhs vector if resolve has...
void use_compressed_row_for_superlu_serial()
Use the compressed row format in superlu serial.
void disable_doc_stats()
Disable documentation of solver statistics.
void * Dist_solver_data_pt
Storage for the LU factors and other data required by SuperLU.
bool Dist_distributed_solve_data_allocated
Flag is true if solve data has been generated for distributed matrix.
int Dist_npcol
Number of columns for the process grid.
void backsub(const DoubleVector &rhs, DoubleVector &result)
Do the backsubstitution for SuperLU solver Note: returns the global result Vector.
void enable_computation_of_gradient()
function to enable the computation of the gradient
void clean_up_memory()
Clean up the memory allocated by the solver.
int Dist_nprow
Number of rows for the process grid.
void backsub_transpose(const DoubleVector &rhs, DoubleVector &result)
Do the back-substitution for transposed system of the SuperLU solver Note: Returns the global result ...
int Serial_info
Info flag for the SuperLU solver.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...