SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/LinearSolver - Gmres.cpp Hit Total Coverage
Commit: ebec864322c50bab8dca0a90baf8d01875114261 Lines: 1 2 50.0 %
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/Gmres.hpp"
       5             : 
       6             : #include <blaze/math/Column.h>
       7             : #include <blaze/math/Submatrix.h>
       8             : #include <blaze/math/Subvector.h>
       9             : #include <blaze/math/lapack/trsv.h>
      10             : 
      11             : #include "DataStructures/DenseMatrix.hpp"
      12             : #include "DataStructures/DenseVector.hpp"
      13             : #include "Utilities/ConstantExpressions.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : 
      16           1 : namespace LinearSolver::gmres::detail {
      17             : 
      18             : namespace {
      19             : 
      20             : void update_givens_rotation(const gsl::not_null<double*> givens_cosine,
      21             :                             const gsl::not_null<double*> givens_sine,
      22             :                             const double rho, const double sigma) noexcept {
      23             :   if (UNLIKELY(rho == 0.)) {
      24             :     *givens_cosine = 0.;
      25             :     *givens_sine = 1.;
      26             :   } else {
      27             :     const double tmp = sqrt(square(rho) + square(sigma));
      28             :     *givens_cosine = abs(rho) / tmp;
      29             :     *givens_sine = *givens_cosine * sigma / rho;
      30             :   }
      31             : }
      32             : 
      33             : template <typename Arg>
      34             : void apply_and_update_givens_rotation(
      35             :     Arg argument, const gsl::not_null<DenseVector<double>*> givens_sine_history,
      36             :     const gsl::not_null<DenseVector<double>*> givens_cosine_history,
      37             :     const size_t iteration) noexcept {
      38             :   const size_t k = iteration + 1;
      39             :   for (size_t i = 0; i < k - 1; ++i) {
      40             :     const double tmp = (*givens_cosine_history)[i] * argument[i] +
      41             :                        (*givens_sine_history)[i] * argument[i + 1];
      42             :     argument[i + 1] = (*givens_cosine_history)[i] * argument[i + 1] -
      43             :                       (*givens_sine_history)[i] * argument[i];
      44             :     argument[i] = tmp;
      45             :   }
      46             :   update_givens_rotation(make_not_null(&(*givens_cosine_history)[k - 1]),
      47             :                          make_not_null(&(*givens_sine_history)[k - 1]),
      48             :                          argument[k - 1], argument[k]);
      49             :   argument[k - 1] = (*givens_cosine_history)[k - 1] * argument[k - 1] +
      50             :                     (*givens_sine_history)[k - 1] * argument[k];
      51             :   argument[k] = 0.;
      52             : }
      53             : 
      54             : }  // namespace
      55             : 
      56             : void solve_minimal_residual(
      57             :     const gsl::not_null<DenseMatrix<double>*> orthogonalization_history,
      58             :     const gsl::not_null<DenseVector<double>*> residual_history,
      59             :     const gsl::not_null<DenseVector<double>*> givens_sine_history,
      60             :     const gsl::not_null<DenseVector<double>*> givens_cosine_history,
      61             :     const size_t iteration) noexcept {
      62             :   residual_history->resize(iteration + 2);
      63             :   givens_sine_history->resize(iteration + 1);
      64             :   givens_cosine_history->resize(iteration + 1);
      65             :   // Givens-rotate the `orthogonalization_history` to eliminate its lower-right
      66             :   // entry, making it an upper triangular matrix.
      67             :   // Thus, we iteratively update the QR decomposition of the Hessenberg matrix
      68             :   // that is built through orthogonalization of the Krylov basis vectors.
      69             :   apply_and_update_givens_rotation(
      70             :       blaze::column(*orthogonalization_history, iteration), givens_sine_history,
      71             :       givens_cosine_history, iteration);
      72             :   // Also Givens-rotate the `residual_history`. The vector excluding its last
      73             :   // entry represents the accumulated Givens-rotation applied to the vector
      74             :   // `(initial_residual, 0, 0, ...)` and the last entry represents the
      75             :   // remaining residual.
      76             :   (*residual_history)[iteration + 1] =
      77             :       -(*givens_sine_history)[iteration] * (*residual_history)[iteration];
      78             :   (*residual_history)[iteration] =
      79             :       (*givens_cosine_history)[iteration] * (*residual_history)[iteration];
      80             : }
      81             : 
      82             : DenseVector<double> minimal_residual_vector(
      83             :     const DenseMatrix<double>& orthogonalization_history,
      84             :     const DenseVector<double>& residual_history) noexcept {
      85             :   const size_t length = orthogonalization_history.columns();
      86             :   DenseVector<double> minres = blaze::subvector(residual_history, 0, length);
      87             :   blaze::trsv(blaze::submatrix(orthogonalization_history, 0, 0, length, length),
      88             :               minres, 'U', 'N', 'N');
      89             :   return minres;
      90             : }
      91             : 
      92             : }  // namespace LinearSolver::gmres::detail

Generated by: LCOV version 1.14