Skip to content

Use modified Gram-Schmidt orthogonalization in GMRES #600

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
michelebucelli wants to merge 5 commits into
base: main
Choose a base branch
from
Draft

Use modified Gram-Schmidt orthogonalization in GMRES #600

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
08a5d0c
GMRES uses modified Gram-Schmidt
michelebucelli Aug 13, 2026
e4e285d
Use increment operators in FourierInterpolation, where possible
michelebucelli Aug 13, 2026
25a7b8a
Merge branch 'main' into fix/modified-gram-schmidt-gmres
michelebucelli Aug 21, 2026
b9b78f1
Fix typo in switch default case statement
michelebucelli Aug 21, 2026
d3c66cc
Simplify GMRES orthogonalization functions by removing the check for …
michelebucelli Aug 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Code/Source/linear_solver/dot.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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++) {
Expand Down
134 changes: 69 additions & 65 deletions Code/Source/linear_solver/gmres.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -20,11 +19,14 @@
#include "Array3.h"
#include "DebugMsg.h"

#include <limits>
#include <math.h>

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;
Expand Down Expand Up @@ -58,15 +60,72 @@ void bc_pre(fsi_linear_solver::FSILS_lhsType& lhs, fsi_linear_solver::FSILS_subL
}
}

/// @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.
///
/// @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<double> &u,
Array<double> &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);
omp_la::omp_mul_s(nNo, 1.0 / h(i + 1, i), w);
}

/// @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.
///
/// @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<double> &u, Array<double> &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));
Comment thread
michelebucelli marked this conversation as resolved.
}

h(i + 1, i) = norm::fsi_ls_normv(dof, mynNo, lhs.commu, w);
omp_la::omp_mul_v(dof, nNo, 1.0 / h(i + 1, i), w);
}

} // 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<double>& Val, const Array<double>& R, Array<double>& 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<double> &Val, const Array<double> &R, Array<double> &X) {
#define n_debug_gmres
#ifdef debug_gmres
DebugMsg dmsg(__func__, lhs.commu.task);
dmsg.banner();
#endif
Expand Down Expand Up @@ -168,25 +227,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);
Comment thread
michelebucelli marked this conversation as resolved.

for (int j = 0; j <= i-1; j++) {
double tmp = c(j)*h(j,i) + s(j)*h(j+1,i);
Expand Down Expand Up @@ -337,25 +378,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);
Expand Down Expand Up @@ -527,26 +550,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);
Expand Down
17 changes: 4 additions & 13 deletions Code/Source/solver/FourierInterpolation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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];
}
}

Expand Down Expand Up @@ -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] -=
Expand Down
Loading