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