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 7 : #include 8 : #include 9 : #include 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 givens_cosine, 21 : const gsl::not_null 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 34 : void apply_and_update_givens_rotation( 35 : Arg argument, const gsl::not_null*> givens_sine_history, 36 : const gsl::not_null*> 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*> orthogonalization_history, 58 : const gsl::not_null*> residual_history, 59 : const gsl::not_null*> givens_sine_history, 60 : const gsl::not_null*> 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 minimal_residual_vector( 83 : const DenseMatrix& orthogonalization_history, 84 : const DenseVector& residual_history) noexcept { 85 : const size_t length = orthogonalization_history.columns(); 86 : DenseVector 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