Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include "Utilities/Gsl.hpp" 7 : 8 : /// \cond 9 : class DataVector; 10 : class Matrix; 11 : /// \endcond 12 : 13 : /// LAPACK wrappers 14 1 : namespace lapack { 15 : /// @{ 16 : /*! 17 : * \ingroup LinearSolverGroup 18 : * \brief Wrapper for LAPACK dgesv, which solves the general linear equation 19 : * \f$A x = b\f$ by LUP (Lower-triangular, upper-triangular, and permutation) 20 : * decomposition. 21 : * \details Several interfaces are provided with combinations of 22 : * input parameters due to the versatility of the LAPACK utility. 23 : * - `rhs_in_solution_out` or `solution` : If the `rhs` `DataVector` is not 24 : * separately supplied, this single pass-by-pointer `DataVector` acts as the 25 : * input \f$ b\f$, and is overwritten with the solution \f$x\f$; in this case 26 : * `rhs_in_solution_out` must be sufficiently large to contain the solution. If 27 : * `rhs` is separately supplied, `solution` will contain \f$x\f$ at the end of 28 : * the algorithm, and must be sufficiently large to contain the solution. 29 : * - `rhs` (optional) : An input `DataVector` representing \f$b\f$ that is not 30 : * overwritten. 31 : * - `matrix_operator`: The `Matrix` \f$A\f$. If passed by pointer, this matrix 32 : * will be 'destroyed' and overwritten with LU decomposition information. 33 : * Optionally, a `pivots` vector may be passed by pointer so that the full LUP 34 : * decomposition information can be recovered from the LAPACK output. If passed 35 : * by const reference, the wrapper will make a copy so that LAPACK does not 36 : * overwrite the supplied matrix. 37 : * - `pivots` (optional) : a `std::vector<int>`, passed by pointer, which will 38 : * contain the permutation information necessary to reassemble the full matrix 39 : * information if LAPACK is permitted to modify the input matrix in-place. 40 : * - `number_of_rhs` : The number of columns of \f$x\f$ and \f$b\f$. This is 41 : * thought of as the 'number of right-hand-sides' to be solved for. In most 42 : * cases, that additional dimension can be inferred from the dimensionality of 43 : * the matrix and the input vector(s), which occurs if `number_of_rhs` is set to 44 : * 0 (default). If the provided `DataVector` is an inappropriate size for that 45 : * inference (i.e. is not a multiple of the largest dimension of the matrix), or 46 : * if you only want to solve for the first `number_of_rhs` columns, the 47 : * parameter `number_of_rhs` must be supplied. 48 : * 49 : * The function return `int` is the value provided by the `INFO` field of the 50 : * LAPACK call. It is 0 for a successful linear solve, and nonzero values code 51 : * for types of failures of the algorithm. 52 : * See LAPACK documentation for further details about `dgesv`: 53 : * http://www.netlib.org/lapack/ 54 : */ 55 1 : int general_matrix_linear_solve(gsl::not_null<DataVector*> rhs_in_solution_out, 56 : gsl::not_null<Matrix*> matrix_operator, 57 : int number_of_rhs = 0); 58 : 59 1 : int general_matrix_linear_solve(gsl::not_null<DataVector*> rhs_in_solution_out, 60 : gsl::not_null<std::vector<int>*> pivots, 61 : gsl::not_null<Matrix*> matrix_operator, 62 : int number_of_rhs = 0); 63 : 64 1 : int general_matrix_linear_solve(gsl::not_null<DataVector*> solution, 65 : gsl::not_null<Matrix*> matrix_operator, 66 : const DataVector& rhs, int number_of_rhs = 0); 67 : 68 1 : int general_matrix_linear_solve(gsl::not_null<DataVector*> solution, 69 : gsl::not_null<std::vector<int>*> pivots, 70 : gsl::not_null<Matrix*> matrix_operator, 71 : const DataVector& rhs, int number_of_rhs = 0); 72 : 73 1 : int general_matrix_linear_solve(gsl::not_null<DataVector*> rhs_in_solution_out, 74 : const Matrix& matrix_operator, 75 : int number_of_rhs = 0); 76 : 77 1 : int general_matrix_linear_solve(gsl::not_null<DataVector*> solution, 78 : const Matrix& matrix_operator, 79 : const DataVector& rhs, int number_of_rhs = 0); 80 : /// @} 81 : } // namespace lapack