From 08a5d0cdc8180c8c1a58107b4224aee0b683d8e6 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 13 Aug 2026 11:58:46 -0500 Subject: [PATCH 1/4] GMRES uses modified Gram-Schmidt --- Code/Source/linear_solver/gmres.cpp | 179 ++++++++++++++++++---------- 1 file changed, 114 insertions(+), 65 deletions(-) diff --git a/Code/Source/linear_solver/gmres.cpp b/Code/Source/linear_solver/gmres.cpp index 4f36cad27..0ab0be855 100644 --- a/Code/Source/linear_solver/gmres.cpp +++ b/Code/Source/linear_solver/gmres.cpp @@ -11,7 +11,6 @@ #include "fsils_api.hpp" #include "add_bc_mul.h" -#include "bcast.h" #include "dot.h" #include "norm.h" #include "omp_la.h" @@ -20,11 +19,14 @@ #include "Array3.h" #include "DebugMsg.h" +#include #include namespace gmres { -void bc_pre(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLsType& ls, const int dof, +namespace { + +void bc_pre(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLsType& ls, const int dof, const int mynNo, const int nNo) { int nsd = dof - 1; @@ -58,15 +60,117 @@ void bc_pre(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subL } } +/// @brief Norm of the newest Krylov vector before it was orthogonalized. +/// +/// Recovered from column 'i' of the Hessenberg matrix: h(j,i) is the component +/// of the vector along the j-th orthonormal basis vector, and h(i+1,i) is the +/// part left orthogonal to all of them. Only squares are summed, so the result +/// is free of cancellation. +/// +/// @param[in] h Hessenberg matrix, with column 'i' already filled. +/// @param[in] i Index of the current Arnoldi step. +/// +/// @return The norm the vector had before orthogonalization. +double pre_orth_norm(const Array &h, const int i) { + double w_norm = h(i + 1, i) * h(i + 1, i); + + for (int j = 0; j <= i; j++) { + w_norm = w_norm + h(j, i) * h(j, i); + } + + return sqrt(w_norm); +} + +/// @brief Orthogonalize the newest Krylov vector against the preceding ones, +/// for one unknown per node. +/// +/// Applies modified Gram-Schmidt to u(:,i+1), filling column 'i' of the +/// Hessenberg matrix and normalizing u(:,i+1) in place. +/// +/// A subdiagonal entry that is negligible against the norm the vector had +/// before orthogonalization is a lucky breakdown: the Krylov space is invariant +/// and already contains the solution. In that case h(i+1,i) is set to zero and +/// the vector is left unnormalized, which drives err(i+1) to zero through the +/// Givens rotations the caller applies next. +/// +/// @param[in] lhs FSILS left-hand side structure, used for its communicator. +/// @param[in] nNo Number of nodes stored on this process, ghost nodes included. +/// @param[in] mynNo Number of nodes owned by this process. +/// @param[in] i Index of the current Arnoldi step. +/// @param[in,out] u Krylov basis. Column i+1 holds the vector to orthogonalize +/// and is overwritten with the new orthonormal basis vector. +/// @param[in,out] h Hessenberg matrix. Column 'i' is overwritten. +void orthogonalize_s(fsi_linear_solver::FSILS_lhsType &lhs, const int nNo, + const int mynNo, const int i, Array &u, + Array &h) { + auto w = u.rcol(i + 1); + + for (int j = 0; j <= i; j++) { + h(j, i) = dot::fsils_dot_s(mynNo, lhs.commu, u.rcol(j), w); + omp_la::omp_sum_s(nNo, -h(j, i), w, u.rcol(j)); + } + + h(i + 1, i) = norm::fsi_ls_norms(mynNo, lhs.commu, w); + + if (h(i + 1, i) > + std::numeric_limits::epsilon() * pre_orth_norm(h, i)) { + omp_la::omp_mul_s(nNo, 1.0 / h(i + 1, i), w); + } else { + h(i + 1, i) = 0.0; + } +} + +/// @brief Orthogonalize the newest Krylov vector against the preceding ones, +/// for 'dof' unknowns per node. +/// +/// Applies modified Gram-Schmidt to u(:,:,i+1), filling column 'i' of the +/// Hessenberg matrix and normalizing u(:,:,i+1) in place. +/// +/// A subdiagonal entry that is negligible against the norm the vector had +/// before orthogonalization is a lucky breakdown: the Krylov space is invariant +/// and already contains the solution. In that case h(i+1,i) is set to zero and +/// the vector is left unnormalized, which drives err(i+1) to zero through the +/// Givens rotations the caller applies next. +/// +/// @param[in] lhs FSILS left-hand side structure, used for its communicator. +/// @param[in] dof Number of unknowns per node. +/// @param[in] nNo Number of nodes stored on this process, ghost nodes included. +/// @param[in] mynNo Number of nodes owned by this process. +/// @param[in] i Index of the current Arnoldi step. +/// @param[in,out] u Krylov basis. Slice i+1 holds the vector to orthogonalize +/// and is overwritten with the new orthonormal basis vector. +/// @param[in,out] h Hessenberg matrix. Column 'i' is overwritten. +void orthogonalize_v(fsi_linear_solver::FSILS_lhsType &lhs, const int dof, + const int nNo, const int mynNo, const int i, + Array3 &u, Array &h) { + auto w = u.rslice(i + 1); + + for (int j = 0; j <= i; j++) { + h(j, i) = dot::fsils_dot_v(dof, mynNo, lhs.commu, u.rslice(j), w); + omp_la::omp_sum_v(dof, nNo, -h(j, i), w, u.rslice(j)); + } + + h(i + 1, i) = norm::fsi_ls_normv(dof, mynNo, lhs.commu, w); + + if (h(i + 1, i) > + std::numeric_limits::epsilon() * pre_orth_norm(h, i)) { + omp_la::omp_mul_v(dof, nNo, 1.0 / h(i + 1, i), w); + } else { + h(i + 1, i) = 0.0; + } +} + +} // namespace + /// @brief Solver the system Val * X = R. /// /// Reproduces the Fortran 'GMRES' subroutine. // -void gmres(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLsType& ls, const int dof, - const Array& Val, const Array& R, Array& X) -{ - #define n_debug_gmres - #ifdef debug_gmres +void gmres(fsi_linear_solver::FSILS_lhsType &lhs, + fsi_linear_solver::FSILS_subLsType &ls, const int dof, + const Array &Val, const Array &R, Array &X) { +#define n_debug_gmres +#ifdef debug_gmres DebugMsg dmsg(__func__, lhs.commu.task); dmsg.banner(); #endif @@ -168,25 +272,7 @@ void gmres(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subLs } } - for (int j = 0; j <= i+1; j++) { - h(j,i) = dot::fsils_nc_dot_v(dof, mynNo, u.rslice(j), u.rslice(i+1)); - } - - // h_col is modofied here so don't use 'rcol() method'. - auto h_col = h.col(i); - bcast::fsils_bcast_v(i+2, h_col, lhs.commu); - h.set_col(i, h_col); - - for (int j = 0; j <= i; j++) { - auto u_slice_1 = u.rslice(i+1); - omp_la::omp_sum_v(dof, nNo, -h(j,i), u_slice_1, u.rslice(j)); - h(i+1,i) = h(i+1,i) - h(j,i)*h(j,i); - } - - h(i+1,i) = sqrt(fabs(h(i+1,i))); - - u_slice_1 = u.rslice(i+1); - omp_la::omp_mul_v(dof, nNo, 1.0/h(i+1,i), u_slice_1); + orthogonalize_v(lhs, dof, nNo, mynNo, i, u, h); for (int j = 0; j <= i-1; j++) { double tmp = c(j)*h(j,i) + s(j)*h(j+1,i); @@ -337,25 +423,7 @@ void gmres_s(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_sub spar_mul::fsils_spar_mul_ss(lhs, lhs.rowPtr, lhs.colPtr, Val, u_col, u_col_1); u.set_col(i+1, u_col_1); - for (int j = 0; j <= i+1; j++) { - h(j,i) = dot::fsils_nc_dot_s(mynNo, u.col(j), u.col(i+1)); - } - - auto h_col = h.col(i); - bcast::fsils_bcast_v(i+2, h_col, lhs.commu); - h.set_col(i, h_col); - - for (int j = 0; j <= i; j++) { - auto u_col_1 = u.col(i+1); - omp_la::omp_sum_s(nNo, -h(j,i), u_col_1, u.col(j)); - u.set_col(i+1, u_col_1); - h(i+1,i) = h(i+1,i) - h(j,i)*h(j,i); - } - h(i+1,i) = sqrt(fabs(h(i+1,i))); - - u_col_1 = u.col(i+1); - omp_la::omp_mul_s(nNo, 1.0/h(i+1,i), u_col_1); - u.set_col(i+1, u_col_1); + orthogonalize_s(lhs, nNo, mynNo, i, u, h); for (int j = 0; j <= i-1; j++) { double tmp = c(j)*h(j,i) + s(j)*h(j+1,i); @@ -527,26 +595,7 @@ void gmres_v(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_sub } } - for (int j = 0; j <= i+1; j++) { - h(j,i) = dot::fsils_nc_dot_v(dof, mynNo, u.rslice(j), u.rslice(i+1)); - #ifdef debug_gmres_v - dmsg << "h(j,i): " << h(j,i); - #endif - } - - auto h_col = h.col(i); - bcast::fsils_bcast_v(i+2, h_col, lhs.commu); - h.set_col(i, h_col); - - for (int j = 0; j <= i; j++) { - auto u_slice_1 = u.rslice(i+1); - omp_la::omp_sum_v(dof, nNo, -h(j,i), u_slice_1, u.rslice(j)); - h(i+1,i) = h(i+1,i) - h(j,i)*h(j,i); - } - h(i+1,i) = sqrt(fabs(h(i+1,i))); - - u_slice_1 = u.rslice(i+1); - omp_la::omp_mul_v(dof, nNo, 1.0/h(i+1,i), u_slice_1); + orthogonalize_v(lhs, dof, nNo, mynNo, i, u, h); for (int j = 0; j <= i-1; j++) { double tmp = c(j)*h(j,i) + s(j)*h(j+1,i); From e4e285db41d09fd787c7981ca440f36ba7722635 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 13 Aug 2026 14:31:46 -0500 Subject: [PATCH 2/4] Use increment operators in FourierInterpolation, where possible --- Code/Source/solver/FourierInterpolation.cpp | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/Code/Source/solver/FourierInterpolation.cpp b/Code/Source/solver/FourierInterpolation.cpp index 88022c917..32b292dd2 100644 --- a/Code/Source/solver/FourierInterpolation.cpp +++ b/Code/Source/solver/FourierInterpolation.cpp @@ -118,9 +118,8 @@ FourierInterpolation FourierInterpolation::from_time_series( for (unsigned int i = 0; i < n_time_points; ++i) { times_shifted[i] -= result.initial_time; for (unsigned int j = 0; j < n_components; ++j) { - values_shifted(j, i) = values_shifted(j, i) - - result.linear_trend_initial_values[j] - - result.linear_trend_slopes[j] * times_shifted[i]; + values_shifted(j, i) -= result.linear_trend_initial_values[j] + + result.linear_trend_slopes[j] * times_shifted[i]; } } @@ -581,16 +580,8 @@ void FourierInterpolation::evaluate_internal(double time, const double K = t * dk; for (int j = 0; j < n_components; ++j) { - // Using value[j] = value[j] + ... instead of value[j] += ..., because - // the latter changes the order of operations enough to break some of - // the tests. - // @todo[michelebucelli] This seems pretty fragile! It happens in a few - // other places in this file as well, where using an increment - // operator (*= or +=) changes the order of operations resulting in - // changes in the test values. This should be investigated and the - // best (i.e. most accurate) operation order should be chosen. - value[j] = value[j] + fourier_coefficients_real(j, i) * std::cos(K) - - fourier_coefficients_imaginary(j, i) * std::sin(K); + value[j] += fourier_coefficients_real(j, i) * std::cos(K) - + fourier_coefficients_imaginary(j, i) * std::sin(K); if (evaluate_derivative) { derivative[j] -= From b9b78f17e55477bdedfcf95b548cde415ccee6cc Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 21 Aug 2026 11:22:48 -0500 Subject: [PATCH 3/4] Fix typo in switch default case statement --- Code/Source/linear_solver/dot.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/linear_solver/dot.cpp b/Code/Source/linear_solver/dot.cpp index 897d5b2bc..2d6f4ad34 100644 --- a/Code/Source/linear_solver/dot.cpp +++ b/Code/Source/linear_solver/dot.cpp @@ -68,7 +68,7 @@ double fsils_dot_v(const int dof, const int nNo, FSILS_commuType& commu, const A } break; - defualt: + default: for (int i = 0; i < nNo; i++) { double sum{0.0}; for (int j = 0; j < U.nrows(); j++) { From d3c66cc4d24869426870fe91ee2582f827e4d910 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 21 Aug 2026 11:37:23 -0500 Subject: [PATCH 4/4] Simplify GMRES orthogonalization functions by removing the check for the "lucky breakdown" --- Code/Source/linear_solver/gmres.cpp | 49 ++--------------------------- 1 file changed, 2 insertions(+), 47 deletions(-) diff --git a/Code/Source/linear_solver/gmres.cpp b/Code/Source/linear_solver/gmres.cpp index 0ab0be855..a94104885 100644 --- a/Code/Source/linear_solver/gmres.cpp +++ b/Code/Source/linear_solver/gmres.cpp @@ -60,39 +60,12 @@ void bc_pre(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subL } } -/// @brief Norm of the newest Krylov vector before it was orthogonalized. -/// -/// Recovered from column 'i' of the Hessenberg matrix: h(j,i) is the component -/// of the vector along the j-th orthonormal basis vector, and h(i+1,i) is the -/// part left orthogonal to all of them. Only squares are summed, so the result -/// is free of cancellation. -/// -/// @param[in] h Hessenberg matrix, with column 'i' already filled. -/// @param[in] i Index of the current Arnoldi step. -/// -/// @return The norm the vector had before orthogonalization. -double pre_orth_norm(const Array &h, const int i) { - double w_norm = h(i + 1, i) * h(i + 1, i); - - for (int j = 0; j <= i; j++) { - w_norm = w_norm + h(j, i) * h(j, i); - } - - return sqrt(w_norm); -} - /// @brief Orthogonalize the newest Krylov vector against the preceding ones, /// for one unknown per node. /// /// Applies modified Gram-Schmidt to u(:,i+1), filling column 'i' of the /// Hessenberg matrix and normalizing u(:,i+1) in place. /// -/// A subdiagonal entry that is negligible against the norm the vector had -/// before orthogonalization is a lucky breakdown: the Krylov space is invariant -/// and already contains the solution. In that case h(i+1,i) is set to zero and -/// the vector is left unnormalized, which drives err(i+1) to zero through the -/// Givens rotations the caller applies next. -/// /// @param[in] lhs FSILS left-hand side structure, used for its communicator. /// @param[in] nNo Number of nodes stored on this process, ghost nodes included. /// @param[in] mynNo Number of nodes owned by this process. @@ -111,13 +84,7 @@ void orthogonalize_s(fsi_linear_solver::FSILS_lhsType &lhs, const int nNo, } h(i + 1, i) = norm::fsi_ls_norms(mynNo, lhs.commu, w); - - if (h(i + 1, i) > - std::numeric_limits::epsilon() * pre_orth_norm(h, i)) { - omp_la::omp_mul_s(nNo, 1.0 / h(i + 1, i), w); - } else { - h(i + 1, i) = 0.0; - } + omp_la::omp_mul_s(nNo, 1.0 / h(i + 1, i), w); } /// @brief Orthogonalize the newest Krylov vector against the preceding ones, @@ -126,12 +93,6 @@ void orthogonalize_s(fsi_linear_solver::FSILS_lhsType &lhs, const int nNo, /// Applies modified Gram-Schmidt to u(:,:,i+1), filling column 'i' of the /// Hessenberg matrix and normalizing u(:,:,i+1) in place. /// -/// A subdiagonal entry that is negligible against the norm the vector had -/// before orthogonalization is a lucky breakdown: the Krylov space is invariant -/// and already contains the solution. In that case h(i+1,i) is set to zero and -/// the vector is left unnormalized, which drives err(i+1) to zero through the -/// Givens rotations the caller applies next. -/// /// @param[in] lhs FSILS left-hand side structure, used for its communicator. /// @param[in] dof Number of unknowns per node. /// @param[in] nNo Number of nodes stored on this process, ghost nodes included. @@ -151,13 +112,7 @@ void orthogonalize_v(fsi_linear_solver::FSILS_lhsType &lhs, const int dof, } h(i + 1, i) = norm::fsi_ls_normv(dof, mynNo, lhs.commu, w); - - if (h(i + 1, i) > - std::numeric_limits::epsilon() * pre_orth_norm(h, i)) { - omp_la::omp_mul_v(dof, nNo, 1.0 / h(i + 1, i), w); - } else { - h(i + 1, i) = 0.0; - } + omp_la::omp_mul_v(dof, nNo, 1.0 / h(i + 1, i), w); } } // namespace