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
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
48namespace 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
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?
123 {
124 return Doc_time;
125 }
126
127 /// Boolean flag indicating if resolves are enabled
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
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 {
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
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 {
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 {
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 {
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 {
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 {
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)
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).
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
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
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:151
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...