linear_solver.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 //This header defines a class for linear solvers
32 
33 //Include guards
34 #ifndef OOMPH_LINEAR_SOLVER_HEADER
35 #define OOMPH_LINEAR_SOLVER_HEADER
36 
37 
38 // Config header generated by autoconfig
39 #ifdef HAVE_CONFIG_H
40  #include <oomph-lib-config.h>
41 #endif
42 
43 #ifdef OOMPH_HAS_MPI
44 #include "mpi.h"
45 #endif
46 
47 //oomph-lib headers
48 #include "Vector.h"
49 #include "double_vector.h"
50 #include "matrices.h"
51 
52 namespace oomph
53 {
54 
55 
56 // Forward declaration of problem class
57 class Problem;
58 
59 //====================================================================
60 /// Base class for all linear solvers. This merely defines standard
61 /// interfaces for linear solvers, so that different solvers can be
62 /// used in a clean and transparent manner. Note that LinearSolvers
63 /// are primarily used to solve the linear systems arising in
64 /// oomph-lib's Newton iteration. Their primary solve function
65 /// therefore takes a pointer to the associated problem, construct its
66 /// Jacobian matrix and residual vector, and return the solution
67 /// of the linear system formed by the Jacobian and the residual vector.
68 /// We also provide broken virtual interfaces
69 /// to a linear-algebra-type solve function in which the matrix
70 /// and the rhs can be specified, but this are not guaranteed to
71 /// implemented for all linear solvers (e.g. for frontal solvers).
72 //====================================================================
74 {
75  protected:
76 
77  ///\short Boolean that indicates whether the matrix (or its factors, in
78  ///the case of direct solver) should be stored so that the resolve function
79  ///can be used.
81 
82  /// \short Boolean flag that indicates whether the time taken
83  // for the solves should be documented
84  bool Doc_time;
85 
86  /// \short flag that indicates whether the gradient required for the
87  /// globally convergent Newton method should be computed or not
89 
90  /// \short flag that indicates whether the gradient was computed or not
92 
93  /// \short DoubleVector storing the gradient for the globally convergent
94  /// Newton method
96 
97  public:
98 
99  /// Empty constructor, initialise the member data
102 
103  /// Broken copy constructor
104  LinearSolver(const LinearSolver& dummy)
105  {
106  BrokenCopy::broken_copy("LinearSolver");
107  }
108 
109  /// Broken assignment operator
110  void operator=(const LinearSolver&)
111  {
112  BrokenCopy::broken_assign("LinearSolver");
113  }
114 
115  /// Empty virtual destructor
116  virtual ~LinearSolver() {}
117 
118  /// \short Enable documentation of solve times
119  void enable_doc_time() {Doc_time=true;}
120 
121  /// \short Disable documentation of solve times
122  void disable_doc_time() {Doc_time=false;}
123 
124  /// \short Is documentation of solve times enabled?
125  bool is_doc_time_enabled() const {return Doc_time;}
126 
127  /// \short Boolean flag indicating if resolves are enabled
128  bool is_resolve_enabled() const {return Enable_resolve;}
129 
130  /// \short Enable resolve (i.e. store matrix and/or LU decomposition, say)
131  /// Virtual so it can be overloaded to perform additional tasks
132  virtual void enable_resolve() {Enable_resolve=true;}
133 
134  /// \short Disable resolve (i.e. store matrix and/or LU decomposition, say)
135  /// This function simply resets an internal flag. It's virtual so
136  /// it can be overloaded to perform additional tasks such as
137  /// cleaning up memory that is only required for the resolve.
138  virtual void disable_resolve() {Enable_resolve=false;}
139 
140  /// \short Solver: Takes pointer to problem and returns the results vector
141  /// which contains the solution of the linear system defined by
142  /// the problem's fully assembled Jacobian and residual vector.
143  virtual void solve(Problem* const &problem_pt, DoubleVector &result)=0;
144 
145  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
146  /// vector and returns the solution of the linear system.
147  virtual void solve(DoubleMatrixBase* const &matrix_pt,
148  const DoubleVector &rhs,
149  DoubleVector &result)
150  {
151  throw OomphLibError(
152  "DoubleVector based solve function not implemented for this solver",
153  OOMPH_CURRENT_FUNCTION,
154  OOMPH_EXCEPTION_LOCATION);
155  }
156 
157  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
158  /// vector and returns the solution of the linear system.
159  virtual void solve(DoubleMatrixBase* const &matrix_pt,
160  const Vector<double> &rhs,
161  Vector<double> &result)
162  {
163  throw OomphLibError(
164  "Vector<double> based solve function not implemented for this solver",
165  OOMPH_CURRENT_FUNCTION,
166  OOMPH_EXCEPTION_LOCATION);
167  }
168 
169  /// \short Resolve the system defined by the last assembled jacobian
170  /// and the rhs vector. Solution is returned in the vector result.
171  /// (broken virtual)
172  virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
173  {
174  throw OomphLibError(
175  "Resolve function not implemented for this linear solver",
176  OOMPH_CURRENT_FUNCTION,
177  OOMPH_EXCEPTION_LOCATION);
178  }
179 
180  /// \short Empty virtual function that can be overloaded in specific
181  /// linear solvers to clean up any memory that may have been
182  /// allocated (e.g. when preparing for a re-solve).
183  virtual void clean_up_memory() {}
184 
185  /// \short returns the time taken to assemble the Jacobian matrix and
186  /// residual vector (needs to be overloaded for each solver)
187  virtual double jacobian_setup_time() const
188  {
189  throw OomphLibError(
190  "jacobian_setup_time has not been implemented for this linear solver",
191  OOMPH_CURRENT_FUNCTION,
192  OOMPH_EXCEPTION_LOCATION);
193  return 0;
194  }
195 
196  /// \short return the time taken to solve the linear system (needs to be
197  /// overloaded for each linear solver)
198  virtual double linear_solver_solution_time() const
199  {
200  throw OomphLibError(
201  "linear_solver_solution_time() not implemented for this linear solver",
202  OOMPH_CURRENT_FUNCTION,
203  OOMPH_EXCEPTION_LOCATION);
204  return 0;
205  }
206 
207  /// \short function to enable the computation of the gradient required
208  /// for the globally convergent Newton method
210  {
211  throw OomphLibError(
212  "enable_computation_of_gradient() not implemented for "
213  "this linear solver",
214  OOMPH_CURRENT_FUNCTION,
215  OOMPH_EXCEPTION_LOCATION);
216  }
217 
218  /// \short function to disable the computation of the gradient required
219  /// for the globally convergent Newton method
221  {
222  Compute_gradient=false;
223  }
224 
225  /// \short function to reset the size of the gradient before each Newton
226  /// solve
228  {
230  }
231 
232  /// \short function to access the gradient, provided it has been computed
233  void get_gradient(DoubleVector& gradient)
234  {
235 #ifdef PARANOID
237  {
238 #endif
240 #ifdef PARANOID
241  }
242  else
243  {
244  throw OomphLibError(
245  "The gradient has not been computed for this linear solver!",
246  OOMPH_CURRENT_FUNCTION,
247  OOMPH_EXCEPTION_LOCATION);
248  }
249 #endif
250  }
251 
252 };
253 
254 //=============================================================================
255 /// \short Dense LU decomposition-based solve of full assembled linear system.
256 /// VERY inefficient but useful to illustrate the principle.
257 /// Only suitable for use with Serial matrices and vectors.
258 /// This solver will only work with non-distributed matrices and vectors
259 /// (note: DenseDoubleMatrix is not distributable)
260 //============================================================================
261 class DenseLU : public LinearSolver
262 {
263  ///The DenseDoubleMatrix class is a friend
264  friend class DenseDoubleMatrix;
265 
266  public:
267 
268  /// Constructor, initialise storage
270  Solution_time(0.0),
272  Index(0), LU_factors(0)
273  {
274  // Shut up!
275  Doc_time=false;
276  }
277 
278  /// Broken copy constructor
279  DenseLU(const DenseLU& dummy)
280  {
281  BrokenCopy::broken_copy("DenseLU");
282  }
283 
284  /// Broken assignment operator
285  void operator=(const DenseLU&)
286  {
287  BrokenCopy::broken_assign("DenseLU");
288  }
289 
290  /// Destructor, clean up the stored LU factors
292 
293  /// \short Solver: Takes pointer to problem and returns the results Vector
294  /// which contains the solution of the linear system defined by
295  /// the problem's fully assembled Jacobian and residual Vector.
296  void solve(Problem* const &problem_pt, DoubleVector &result);
297 
298  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
299  /// vector and returns the solution of the linear system.
300  void solve(DoubleMatrixBase* const &matrix_pt,const DoubleVector &rhs,
301  DoubleVector &result);
302 
303  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
304  /// vector and returns the solution of the linear system.
305  void solve(DoubleMatrixBase* const &matrix_pt,const Vector<double> &rhs,
306  Vector<double> &result);
307 
308  /// \short returns the time taken to assemble the jacobian matrix and
309  /// residual vector
310  double jacobian_setup_time() const
311  {
312  return Jacobian_setup_time;
313  }
314 
315  /// \short return the time taken to solve the linear system (needs to be
316  /// overloaded for each linear solver)
317  virtual double linear_solver_solution_time() const
318  {
319  return Solution_time;
320  }
321 
322  protected:
323 
324  /// Perform the LU decomposition of the matrix
325  void factorise(DoubleMatrixBase* const &matrix_pt);
326 
327  /// Do the backsubstitution step to solve the system LU result = rhs
328  void backsub(const DoubleVector &rhs,
329  DoubleVector &result);
330 
331  /// perform back substitution using Vector<double>
332  void backsub(const Vector<double> &rhs, Vector<double> &result);
333 
334  /// Clean up the stored LU factors
335  void clean_up_memory();
336 
337  /// Jacobian setup time
339 
340  /// Solution time
342 
343  /// \short Sign of the determinant of the matrix
344  /// (obtained during the LU decomposition)
346 
347  private:
348 
349  /// Pointer to storage for the index of permutations in the LU solve
350  long* Index;
351 
352  /// Pointer to storage for the LU decomposition
353  double* LU_factors;
354 };
355 
356 
357 //====================================================================
358 /// \short Dense LU decomposition-based solve of linear system
359 /// assembled via finite differencing of the residuals Vector.
360 /// Even more inefficient than DenseLU but excellent sanity check!
361 //====================================================================
362 class FD_LU : public DenseLU
363 {
364  public:
365 
366  /// Constructor: empty
367  FD_LU() : DenseLU() {}
368 
369  /// Broken copy constructor
370  FD_LU(const FD_LU& dummy)
371  {
372  BrokenCopy::broken_copy("FD_LU");
373  }
374 
375  /// Broken assignment operator
376  void operator=(const FD_LU&)
377  {
378  BrokenCopy::broken_assign("FD_LU");
379  }
380 
381  /// \short Solver: Takes pointer to problem and returns the results Vector
382  /// which contains the solution of the linear system defined by
383  /// the problem's residual Vector (Jacobian computed by FD approx.)
384  void solve(Problem* const &problem_pt, DoubleVector &result);
385 
386  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
387  /// vector and returns the solution of the linear system.
388  void solve(DoubleMatrixBase* const &matrix_pt,
389  const DoubleVector &rhs,
390  DoubleVector &result)
391  {DenseLU::solve(matrix_pt,rhs,result);}
392 
393  /// \short Linear-algebra-type solver: Takes pointer to a matrix
394  /// and rhs vector and returns the solution of the linear system
395  /// Call the broken base-class version. If you want this, please
396  /// implement it
397  void solve(DoubleMatrixBase* const &matrix_pt,
398  const Vector<double> &rhs,
399  Vector<double> &result)
400  {LinearSolver::solve(matrix_pt,rhs,result);}
401 
402 };
403 
404 
405 
406 
407 ///////////////////////////////////////////////////////////////////////////////
408 ///////////////////////////////////////////////////////////////////////////////
409 ///////////////////////////////////////////////////////////////////////////////
410 
411 
412 
413 
414 //=============================================================================
415 /// \short. SuperLU Project Solver class. This is a combined wrapper for both
416 /// SuperLU and SuperLU Dist.
417 /// See http://crd.lbl.gov/~xiaoye/SuperLU/
418 /// Default Behaviour: If this solver is distributed over more than one
419 /// processor then SuperLU Dist is used.
420 /// Member data naming convention: member data associated with the SuperLU
421 /// Dist solver begins Dist_... and member data associated with the serial
422 /// SuperLU solver begins Serial_... .
423 //=============================================================================
425 {
426 
427  public:
428 
429  /// \short enum type to specify the solver behaviour.
430  /// Default - will employ superlu dist if more than 1 processor.
431  /// Serial - will always try use superlu (serial).
432  /// Distributed - will always try to use superlu dist.
434 
435  /// Constructor. Set the defaults.
437  : Serial_f_factors(0)
438  {
439  // Set solver wide default values and null pointers
440  Doc_stats=false;
441  Suppress_solve=false;
442  Using_dist = false;
444 
445 #ifdef OOMPH_HAS_MPI
446  // Set default values and nullify pointers for SuperLU Dist
453  Dist_value_pt=0;
454  Dist_index_pt=0;
455  Dist_start_pt=0;
456 #endif
457 
458  // Set default values and null pointers for SuperLU (serial)
461  Serial_n_dof=0;
462  }
463 
464  /// Broken copy constructor
466  {
467  BrokenCopy::broken_copy("SuperLUSolver");
468  }
469 
470  /// Broken assignment operator
471  void operator=(const SuperLUSolver&)
472  {
473  BrokenCopy::broken_assign("SuperLUSolver");
474  }
475 
476  ///Destructor, clean up the stored matrices
478  {
479  clean_up_memory();
480  }
481 
482  /// function to enable the computation of the gradient
484  {
485  Compute_gradient=true;
486  }
487 
488  // SuperLUSolver methods
489  ////////////////////////
490 
491  /// Overload disable resolve so that it cleans up memory too
493  {
495  clean_up_memory();
496  }
497 
498  /// \short Solver: Takes pointer to problem and returns the results Vector
499  /// which contains the solution of the linear system defined by
500  /// the problem's fully assembled Jacobian and residual Vector.
501  void solve(Problem* const &problem_pt, DoubleVector &result);
502 
503  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
504  /// vector and returns the solution of the linear system.
505  /// The function returns the global result Vector.
506  /// Note: if Delete_matrix_data is true the function
507  /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
508  void solve(DoubleMatrixBase* const &matrix_pt,
509  const DoubleVector &rhs,
510  DoubleVector &result);
511 
512 
513 /* /// \short Linear-algebra-type solver: Takes pointer to a matrix */
514 /* /// and rhs vector and returns the solution of the linear system */
515 /* /// Call the broken base-class version. If you want this, please */
516 /* /// implement it */
517 /* void solve(DoubleMatrixBase* const &matrix_pt, */
518 /* const Vector<double> &rhs, */
519 /* Vector<double> &result) */
520 /* {LinearSolver::solve(matrix_pt,rhs,result);} */
521 
522 
523  /// \short Resolve the system defined by the last assembled jacobian
524  /// and the specified rhs vector if resolve has been enabled.
525  /// Note: returns the global result Vector.
526  void resolve(const DoubleVector &rhs, DoubleVector &result);
527 
528  /// Enable documentation of solver statistics
529  void enable_doc_stats() {Doc_stats = true;}
530 
531  /// Disable documentation of solver statistics
532  void disable_doc_stats() {Doc_stats = false;}
533 
534  /// \short returns the time taken to assemble the jacobian matrix and
535  /// residual vector
536  double jacobian_setup_time() const
537  {
538  return Jacobian_setup_time;
539  }
540 
541  /// \short return the time taken to solve the linear system (needs to be
542  /// overloaded for each linear solver)
543  virtual double linear_solver_solution_time() const
544  {
545  return Solution_time;
546  }
547 
548  /// \short Do the factorisation stage
549  /// Note: if Delete_matrix_data is true the function
550  /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
551  void factorise(DoubleMatrixBase* const &matrix_pt);
552 
553  /// \short Do the backsubstitution for SuperLU solver
554  /// Note: returns the global result Vector.
555  void backsub(const DoubleVector &rhs,
556  DoubleVector &result);
557 
558  /// Clean up the memory allocated by the solver
559  void clean_up_memory();
560 
561  /// \short Specify the solve type. Either default, serial or distributed.
562  /// See enum SuperLU_solver_type for more details.
563  void set_solver_type(const Type& t)
564  {
565  this->clean_up_memory();
566  Solver_type = t;
567  }
568 
569  // SuperLU (serial) methods
570  ///////////////////////////
571 
572  /// Use the compressed row format in superlu serial
575 
576  /// Use the compressed column format in superlu serial
579 
580 #ifdef OOMPH_HAS_MPI
581 
582  // SuperLU Dist methods
583  ///////////////////////
584 
585  /// \short Set Delete_matrix_data flag. SuperLU_dist needs its own copy
586  /// of the input matrix, therefore a copy must be made if any matrix
587  /// used with this solver is to be preserved. If the input matrix can be
588  /// deleted the flag can be set to true to reduce the amount of memory
589  /// required, and the matrix data will be wiped using its clean_up_memory()
590  /// function. Default value is false.
593 
594  /// \short Unset Delete_matrix_data flag. SuperLU_dist needs its own copy
595  /// of the input matrix, therefore a copy must be made if any matrix
596  /// used with this solver is to be preserved. If the input matrix can be
598  {Dist_delete_matrix_data=false;}
599 
600  /// \short Set flag so that SuperLU_DIST is allowed to permute matrix rows
601  /// and columns during factorisation. This is the default for SuperLU_DIST,
602  /// and can lead to significantly faster solves.
605 
606  /// \short Set flag so that SuperLU_DIST is not allowed to permute matrix rows
607  /// and columns during factorisation.
610 
611  /// \short Calling this method will ensure that when the problem based
612  /// solve interface is used, a global (serial) jacobian will be
613  /// assembled.
614  /// Note: calling this function will delete any distributed solve data.
616  {
618  {
619  clean_up_memory();
621  }
622  }
623 
624  /// \short Calling this method will ensure that when the problem based
625  /// solve interface is used, a distributed jacobian will be
626  /// assembled.
627  /// Note: calling this function will delete any global solve data.
629  {
631  {
632  clean_up_memory();
634  }
635  }
636 
637 #endif
638 
639  private:
640 
641  /// factorise method for SuperLU (serial)
642  void factorise_serial(DoubleMatrixBase* const &matrix_pt);
643 
644  /// backsub method for SuperLU (serial)
645  void backsub_serial(const DoubleVector &rhs,
646  DoubleVector &result);
647 
648 #ifdef OOMPH_HAS_MPI
649  /// factorise method for SuperLU Dist
650  void factorise_distributed(DoubleMatrixBase* const &matrix_pt);
651 
652  /// backsub method for SuperLU Dist
653  void backsub_distributed(const DoubleVector &rhs,
654  DoubleVector &result);
655 #endif
656 
657  // SuperLUSolver member data
658  ////////////////////////////
659 
660  /// Jacobian setup time
662 
663  /// Solution time
665 
666  /// Suppress solve?
668 
669  /// Set to true to output statistics (false by default).
670  bool Doc_stats;
671 
672  /// the solver type. see SuperLU_solver_type for details.
674 
675  /// boolean flag indicating whether superlu dist is being used
677 
678  // SuperLU (serial) member data
679  ///////////////////////////////
680 
681  /// Storage for the LU factors as required by SuperLU
683 
684  /// Info flag for the SuperLU solver
686 
687  /// The number of unknowns in the linear system
688  unsigned long Serial_n_dof;
689 
690  /// Sign of the determinant of the matrix
692 
693  /// Use compressed row version?
695 
696 #ifdef OOMPH_HAS_MPI
697 
698  // SuperLU Dist member data
699  ///////////////////////////
700 public:
701  /// \short Static flag that determines whether the warning about
702  /// incorrect distribution of RHSs will be printed or not
704 
705 private:
706  /// \short Flag that determines whether the MPIProblem based solve function
707  /// uses the global or distributed version of SuperLU_DIST
708  /// (default value is false).
710 
711  /// Flag is true if solve data has been generated for a global matrix
713 
714  /// Flag is true if solve data has been generated for distributed matrix
716 
717  /// Storage for the LU factors and other data required by SuperLU
719 
720  /// Number of rows for the process grid
722 
723  /// Number of columns for the process grid
725 
726  /// Info flag for the SuperLU solver
728 
729  /// \short If true then SuperLU_DIST is allowed to permute matrix rows
730  /// and columns during factorisation. This is the default for SuperLU_DIST,
731  /// and can lead to significantly faster solves, but has been known to
732  /// fail, hence the default value is 0.
734 
735  /// \short Delete_matrix_data flag. SuperLU_dist needs its own copy
736  /// of the input matrix, therefore a copy must be made if any matrix
737  /// used with this solver is to be preserved. If the input matrix can be
738  /// deleted the flag can be set to true to reduce the amount of memory
739  /// required, and the matrix data will be wiped using its clean_up_memory()
740  /// function. Default value is false.
742 
743  /// Pointer for storage of the matrix values required by SuperLU_DIST
744  double *Dist_value_pt;
745 
746  /// \short Pointer for storage of matrix rows or column indices required
747  /// by SuperLU_DIST
749 
750  /// \short Pointers for storage of matrix column or row starts
751  // required by SuperLU_DIST
753 
754 #endif
755 }; // end of SuperLUSolver
756 } // end of oomph namespace
757 #endif
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...
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...
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
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 ...
void clean_up_memory()
Clean up the memory allocated by the solver.
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 use_global_solve_in_superlu_dist()
Calling this method will ensure that when the problem based solve interface is used, a global (serial) jacobian will be assembled. Note: calling this function will delete any distributed solve data.
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 broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void * Dist_solver_data_pt
Storage for the LU factors and other data required by SuperLU.
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_distributed(DoubleMatrixBase *const &matrix_pt)
factorise method for SuperLU Dist
void clear()
wipes the DoubleVector
int Sign_of_determinant_of_matrix
Sign of the determinant of the matrix (obtained during the LU decomposition)
bool Using_dist
boolean flag indicating whether superlu dist is being used
void set_solver_type(const Type &t)
Specify the solve type. Either default, serial or distributed. See enum SuperLU_solver_type for more ...
bool Dist_global_solve_data_allocated
Flag is true if solve data has been generated for a global matrix.
int Serial_info
Info flag for the SuperLU solver.
void backsub_distributed(const DoubleVector &rhs, DoubleVector &result)
backsub method for SuperLU Dist
bool Dist_delete_matrix_data
Delete_matrix_data flag. SuperLU_dist needs its own copy of the input matrix, therefore a copy must b...
~SuperLUSolver()
Destructor, clean up the stored matrices.
DenseLU()
Constructor, initialise storage.
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:80
The Problem class.
Definition: problem.h:152
int * Dist_start_pt
Pointers for storage of matrix column or row starts.
int Serial_sign_of_determinant_of_matrix
Sign of the determinant of the matrix.
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...
bool Gradient_has_been_computed
flag that indicates whether the gradient was computed or not
Definition: linear_solver.h:91
char t
Definition: cfortran.h:572
bool Dist_distributed_solve_data_allocated
Flag is true if solve data has been generated for distributed matrix.
bool is_resolve_enabled() const
Boolean flag indicating if resolves are enabled.
virtual double linear_solver_solution_time() const
return the time taken to solve the linear system (needs to be overloaded for each linear solver) ...
. SuperLU Project Solver class. This is a combined wrapper for both SuperLU and SuperLU Dist...
void enable_doc_stats()
Enable documentation of solver statistics.
double * LU_factors
Pointer to storage for the LU decomposition.
int * Dist_index_pt
Pointer for storage of matrix rows or column indices required by SuperLU_DIST.
void operator=(const LinearSolver &)
Broken assignment operator.
bool Doc_stats
Set to true to output statistics (false by default).
void backsub_serial(const DoubleVector &rhs, DoubleVector &result)
backsub method for SuperLU (serial)
double jacobian_setup_time() const
returns the time taken to assemble the jacobian matrix and residual vector
void disable_delete_matrix_data_in_superlu_dist()
Unset Delete_matrix_data flag. SuperLU_dist needs its own copy of the input matrix, therefore a copy must be made if any matrix used with this solver is to be preserved. If the input matrix can be.
void get_gradient(DoubleVector &gradient)
function to access the gradient, provided it has been computed
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(const DoubleVector &rhs, DoubleVector &result)
Do the backsubstitution step to solve the system LU result = rhs.
void reset_gradient()
function to reset the size of the gradient before each Newton solve
double jacobian_setup_time() const
returns the time taken to assemble the jacobian matrix and residual vector
Dense LU decomposition-based solve of full assembled linear system. VERY inefficient but useful to il...
int Dist_nprow
Number of rows for the process grid.
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
void use_distributed_solve_in_superlu_dist()
Calling this method will ensure that when the problem based solve interface is used, a distributed jacobian will be assembled. Note: calling this function will delete any global solve data.
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:84
FD_LU()
Constructor: empty.
void factorise(DoubleMatrixBase *const &matrix_pt)
Perform the LU decomposition of the matrix.
DoubleVector Gradient_for_glob_conv_newton_solve
DoubleVector storing the gradient for the globally convergent Newton method.
Definition: linear_solver.h:95
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
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...
Dense LU decomposition-based solve of linear system assembled via finite differencing of the residual...
double Jacobian_setup_time
Jacobian setup time.
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 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 Dist_allow_row_and_col_permutations
If true then SuperLU_DIST is allowed to permute matrix rows and columns during factorisation. This is the default for SuperLU_DIST, and can lead to significantly faster solves, but has been known to fail, hence the default value is 0.
int Dist_info
Info flag for the SuperLU solver.
void enable_computation_of_gradient()
function to enable the computation of the gradient
void disable_doc_time()
Disable documentation of solve times.
void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the specified rhs vector if resolve has...
Type Solver_type
the solver type. see SuperLU_solver_type for details.
LinearSolver(const LinearSolver &dummy)
Broken copy constructor.
LinearSolver()
Empty constructor, initialise the member data.
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 enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
double Solution_time
Solution time.
double * Dist_value_pt
Pointer for storage of the matrix values required by SuperLU_DIST.
int Dist_npcol
Number of columns for the process grid.
void clean_up_memory()
Clean up the stored LU factors.
void disable_doc_stats()
Disable documentation of solver statistics.
void operator=(const FD_LU &)
Broken assignment operator.
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
virtual void clean_up_memory()
Empty virtual function that can be overloaded in specific linear solvers to clean up any memory that ...
virtual void enable_computation_of_gradient()
function to enable the computation of the gradient required for the globally convergent Newton method...
virtual ~LinearSolver()
Empty virtual destructor.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
DenseLU(const DenseLU &dummy)
Broken copy constructor.
unsigned long Serial_n_dof
The number of unknowns in the linear system.
bool Suppress_solve
Suppress solve?
bool is_doc_time_enabled() const
Is documentation of solve times enabled?
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...
Type
enum type to specify the solver behaviour. Default - will employ superlu dist if more than 1 processo...
FD_LU(const FD_LU &dummy)
Broken copy constructor.
void backsub(const DoubleVector &rhs, DoubleVector &result)
Do the backsubstitution for SuperLU solver Note: returns the global result Vector.
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 operator=(const DenseLU &)
Broken assignment operator.
void operator=(const SuperLUSolver &)
Broken assignment operator.
double Jacobian_setup_time
Jacobian setup time.
static bool Suppress_incorrect_rhs_distribution_warning_in_resolve
Static flag that determines whether the warning about incorrect distribution of RHSs will be printed ...
void * Serial_f_factors
Storage for the LU factors as required by SuperLU.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
SuperLUSolver(const SuperLUSolver &dummy)
Broken copy constructor.
bool Compute_gradient
flag that indicates whether the gradient required for the globally convergent Newton method should be...
Definition: linear_solver.h:88
void use_compressed_row_for_superlu_serial()
Use the compressed row format in superlu serial.
void enable_doc_time()
Enable documentation of solve times.
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 ...
double Solution_time
Solution time.
bool Dist_use_global_solver
Flag that determines whether the MPIProblem based solve function uses the global or distributed versi...
Abstract base class for matrices of doubles – adds abstract interfaces for solving, LU decomposition and multiplication by vectors.
Definition: matrices.h:275
~DenseLU()
Destructor, clean up the stored LU factors.
void use_compressed_column_for_superlu_serial()
Use the compressed column format in superlu serial.
void disable_computation_of_gradient()
function to disable the computation of the gradient required for the globally convergent Newton metho...
bool Serial_compressed_row_flag
Use compressed row version?
void solve(DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the ...
long * Index
Pointer to storage for the index of permutations in the LU solve.
SuperLUSolver()
Constructor. Set the defaults.
void factorise_serial(DoubleMatrixBase *const &matrix_pt)
factorise method for SuperLU (serial)