SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/LinearSolver - Lapack.cpp Hit Total Coverage
Commit: ebec864322c50bab8dca0a90baf8d01875114261 Lines: 1 9 11.1 %
Date: 2020-11-25 20:28:50
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #include "NumericalAlgorithms/LinearSolver/Lapack.hpp"
       5             : 
       6             : #include <cstddef>
       7             : #include <vector>
       8             : 
       9             : #include "DataStructures/DataVector.hpp"
      10             : #include "DataStructures/Matrix.hpp"
      11             : #include "ErrorHandling/Assert.hpp"
      12             : #include "Utilities/Gsl.hpp"
      13             : 
      14             : extern "C" {
      15             : #pragma GCC diagnostic push
      16             : #pragma GCC diagnostic ignored "-Wredundant-decls"
      17           0 : extern void dgesv_(int*, int*, double*, int*, int*, double*, int*,  // NOLINT
      18             :                    int*);
      19             : #pragma GCC diagnostic pop
      20             : }
      21             : 
      22           1 : namespace lapack {
      23           0 : int general_matrix_linear_solve(
      24             :     const gsl::not_null<DataVector*> rhs_in_solution_out,
      25             :     const gsl::not_null<Matrix*> matrix_operator,
      26             :     const int number_of_rhs) noexcept {
      27             :   const size_t rhs_vector_size = matrix_operator->rows();
      28             :   std::vector<int> ipiv(rhs_vector_size);
      29             :   return general_matrix_linear_solve(rhs_in_solution_out, make_not_null(&ipiv),
      30             :                                      matrix_operator, number_of_rhs);
      31             : }
      32             : 
      33           0 : int general_matrix_linear_solve(
      34             :     const gsl::not_null<DataVector*> rhs_in_solution_out,
      35             :     const gsl::not_null<std::vector<int>*> pivots,
      36             :     const gsl::not_null<Matrix*> matrix_operator, int number_of_rhs) noexcept {
      37             :   int output_vector_size = matrix_operator->columns();
      38             :   int rhs_vector_size = matrix_operator->rows();
      39             :   int matrix_spacing = matrix_operator->spacing();
      40             :   ASSERT(output_vector_size == rhs_vector_size,
      41             :          "The LAPACK-based general linear solve requires a square matrix "
      42             :          "input, not "
      43             :              << rhs_vector_size << " by " << output_vector_size);
      44             : 
      45             :   if (number_of_rhs == 0) {
      46             :     ASSERT(
      47             :         (rhs_in_solution_out->size() % matrix_operator->rows()) ==
      48             :             0_st,
      49             :         "The provided DataVector does not have size equal to (number of "
      50             :         "equations) * (larger matrix dimension), so the number of right-hand "
      51             :         "sides cannot be inferred and must be provided explicitly");
      52             :     number_of_rhs = static_cast<int>(
      53             :         rhs_in_solution_out->size() /
      54             :         std::max(matrix_operator->rows(), matrix_operator->columns()));
      55             :   }
      56             :   int info = 0;
      57             :   ASSERT(static_cast<size_t>(rhs_vector_size * number_of_rhs) <=
      58             :                  static_cast<size_t>(rhs_in_solution_out->size()) and
      59             :              static_cast<size_t>(output_vector_size * number_of_rhs) <=
      60             :                  static_cast<size_t>(rhs_in_solution_out->size()),
      61             :          "The single DataVector passed to the LAPACK call must be sufficiently "
      62             :          "large to contain x and b in A x = b");
      63             :   dgesv_(&rhs_vector_size, &number_of_rhs, matrix_operator->data(),
      64             :          &matrix_spacing, pivots->data(), rhs_in_solution_out->data(),
      65             :          &output_vector_size, &info);
      66             :   return info;
      67             : }
      68             : 
      69           0 : int general_matrix_linear_solve(const gsl::not_null<DataVector*> solution,
      70             :                                 const gsl::not_null<Matrix*> matrix_operator,
      71             :                                 const DataVector& rhs,
      72             :                                 const int number_of_rhs) noexcept {
      73             :   const size_t rhs_vector_size = matrix_operator->rows();
      74             :   std::vector<int> ipiv(rhs_vector_size);
      75             :   return general_matrix_linear_solve(solution, make_not_null(&ipiv),
      76             :                                      matrix_operator, rhs, number_of_rhs);
      77             : }
      78             : 
      79           0 : int general_matrix_linear_solve(const gsl::not_null<DataVector*> solution,
      80             :                                 const gsl::not_null<std::vector<int>*> pivots,
      81             :                                 const gsl::not_null<Matrix*> matrix_operator,
      82             :                                 const DataVector& rhs,
      83             :                                 int number_of_rhs) noexcept {
      84             :   // NOLINTNEXTLINE(clang-analyzer-deadcode)
      85             :   const int output_vector_size = matrix_operator->columns();
      86             :   const int rhs_vector_size = matrix_operator->rows();
      87             :   if (number_of_rhs == 0) {
      88             :     ASSERT(rhs.size() % matrix_operator->rows() == 0,
      89             :            "The provided DataVector does not have size equal to (number of "
      90             :            "equations) * (number_of_matrix_rows), so the number of right-hand "
      91             :            "sides cannot be inferred and must be provided explicitly");
      92             :     number_of_rhs = static_cast<int>(rhs.size() / matrix_operator->rows());
      93             :   }
      94             :   ASSERT(solution->size() >=
      95             :              static_cast<size_t>(number_of_rhs) * matrix_operator->columns(),
      96             :          "The provided pointer for the output of the LAPACK linear solve is "
      97             :          "too small for the operation. Solution size is: "
      98             :              << solution->size() << " and should be: "
      99             :              << number_of_rhs * output_vector_size << ".");
     100             :   ASSERT(rhs.size() >=
     101             :              static_cast<size_t>(number_of_rhs) * matrix_operator->rows(),
     102             :          "The provided pointer for the input right-hand side vector of the "
     103             :          "LAPACK linear solve is too small for the operation. Vector size is: "
     104             :              << rhs.size()
     105             :              << " and should be: " << number_of_rhs * rhs_vector_size << ".");
     106             :   std::copy(rhs.begin(), rhs.begin() + number_of_rhs * rhs_vector_size,
     107             :             solution->begin());
     108             :   return general_matrix_linear_solve(solution, pivots, matrix_operator,
     109             :                                      number_of_rhs);
     110             : }
     111             : 
     112           0 : int general_matrix_linear_solve(
     113             :     const gsl::not_null<DataVector*> rhs_in_solution_out,
     114             :     const Matrix& matrix_operator, const int number_of_rhs) noexcept {
     115             :   // LAPACK is permitted to modify the matrix in-place, so we copy before
     116             :   // providing the operator if the original must be preserved.
     117             :   Matrix copied_matrix_operator = matrix_operator;
     118             :   return general_matrix_linear_solve(rhs_in_solution_out,
     119             :                                      make_not_null(&copied_matrix_operator),
     120             :                                      number_of_rhs);
     121             : }
     122             : 
     123           0 : int general_matrix_linear_solve(const gsl::not_null<DataVector*> solution,
     124             :                                 const Matrix& matrix_operator,
     125             :                                 const DataVector& rhs,
     126             :                                 const int number_of_rhs) noexcept {
     127             :   // LAPACK is permitted to modify the matrix in-place, so we copy before
     128             :   // providing the operator if the original must be preserved.
     129             :   Matrix copied_matrix_operator = matrix_operator;
     130             :   return general_matrix_linear_solve(
     131             :       solution, make_not_null(&copied_matrix_operator), rhs, number_of_rhs);
     132             : }
     133             : 
     134             : }  // namespace lapack

Generated by: LCOV version 1.14