Skip to content

Feature/fgmres update #389

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
ASunil-17 wants to merge 9 commits into
base: develop
Choose a base branch
from
Draft

Feature/fgmres update #389

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
6f647a6
Zero norm breakdown handling
Oct 22, 2025
8ddafaf
Computes updated for Ax = b when residual norm is below RHS
Oct 24, 2025
87d7f56
Apply pre-commmit fixes
ASunil-17 Oct 24, 2025
242fbf7
restoring GS
Oct 24, 2025
22746d6
Update feature added to Rand and some cleanup for FGMRES
Oct 24, 2025
472083c
Apply pre-commmit fixes
ASunil-17 Oct 24, 2025
e530660
solving scaled RHS FGMRES
Nov 1, 2025
f6eaa31
Changed GramSchmidt for better error handling
Nov 19, 2025
d9347c6
Reverted CMake back to original
Nov 19, 2025
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
111 changes: 90 additions & 21 deletions resolve/LinSolverIterativeFGMRES.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,37 +108,59 @@ namespace ReSolve

// io::Logger::setVerbosity(io::Logger::EVERYTHING);

int outer_flag = 1;
int notconv = 1;
int i = 0;
int it = 0;
int j = 0;
int k = 0;
int k1 = 0;
int outer_flag = 1;
int update_flag = 1;
int notconv = 1;
int i = 0;
int it = 0;
int j = 0;
int k = 0;
int k1 = 0;

real_type t = 0.0;
real_type rnorm = 0.0;
real_type bnorm = 0.0;
real_type tolrel;
real_type rhsnorm = 0.0;
vector_type vec_v(n_);
vector_type vec_z(n_);
// V[0] = b-A*x_0
// debug

// Compute the residual
vec_R_->setToZero(memspace_);
vec_R_->copyDataFrom(rhs, memspace_, memspace_);
matrix_handler_->matvec(A_, x, vec_R_, &MINUS_ONE, &ONE, memspace_);

// Arnoldi Basis
vec_Z_->setToZero(memspace_);
vec_V_->setToZero(memspace_);

rhs->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_);
// Initializing Residual and Update
vec_Y_->setToZero(memspace_);

// Computing the first Arnodi basis vector
vec_R_->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
rnorm = 0.0;
bnorm = vector_handler_->dot(rhs, rhs, memspace_);
bnorm = vector_handler_->dot(vec_R_, vec_R_, memspace_);
rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
// rnorm = ||V_1||
rnorm = std::sqrt(rnorm);

// Checking if rnorm > norm of RHS
rhsnorm = vector_handler_->dot(rhs, rhs, memspace_);
rhsnorm = std::sqrt(rhsnorm);
if (rnorm > rhsnorm)
{
std::cout << "Initial guess is invalid. Refining Ax = b instead of Ay = r" << std::endl;
Copy link
Collaborator

@shakedregev shakedregev Oct 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the initial guess is bad, throw it out. Do not refine it.

vec_Y_->copyDataFrom(x, memspace_, memspace_);
vec_R_->copyDataFrom(rhs, memspace_, memspace_);
update_flag = 0;
}

bnorm = std::sqrt(bnorm);
io::Logger::misc() << "it 0: norm of residual "
<< std::scientific << std::setprecision(16)
<< rnorm << " Norm of rhs: " << bnorm << "\n";
initial_residual_norm_ = rnorm;

while (outer_flag)
{
if (it == 0)
Expand Down Expand Up @@ -169,6 +191,7 @@ namespace ReSolve
outer_flag = 0;
final_residual_norm_ = rnorm;
initial_residual_norm_ = rnorm;
update_flag = 1;
total_iters_ = 0;
break;
}
Expand Down Expand Up @@ -202,12 +225,19 @@ namespace ReSolve
// V_{i+1}=A*Z_i

vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_);

vec_v.setToZero(memspace_);
matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_);

// orthogonalize V[i+1], form a column of h_H_

GS_->orthogonalize(n_, vec_V_, h_H_, i);
int gs_status = GS_->orthogonalize(n_, vec_V_, h_H_, i);

if (gs_status != 0) // checking for successful breakdown
{
notconv = 0; // exiting outer loop
break;
}

if (i != 0)
{
for (index_type k = 1; k <= i; k++)
Expand Down Expand Up @@ -271,7 +301,7 @@ namespace ReSolve
for (j = 0; j <= i; j++)
{
vec_z.setData(vec_Z_->getData(j, memspace_), memspace_);
vector_handler_->axpy(&h_rs_[j], &vec_z, x, memspace_);
vector_handler_->axpy(&h_rs_[j], &vec_z, vec_Y_, memspace_);
}
}
else
Expand All @@ -287,20 +317,20 @@ namespace ReSolve

vec_v.setData(vec_V_->getData(memspace_), memspace_);
this->precV(&vec_z, &vec_v);
// and add to x
vector_handler_->axpy(&ONE, &vec_v, x, memspace_);
// and add to update
vector_handler_->axpy(&ONE, &vec_v, vec_Y_, memspace_);
}

/* test solution */
/* test update */

if (rnorm <= tolrel || it >= maxit_)
{
// rnorm_aux = rnorm;
outer_flag = 0;
}

rhs->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_);
vec_R_->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
matrix_handler_->matvec(A_, vec_Y_, vec_V_, &MINUS_ONE, &ONE, memspace_);
rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
// rnorm = ||V_1||
rnorm = std::sqrt(rnorm);
Expand All @@ -314,6 +344,36 @@ namespace ReSolve
<< rnorm << "\n";
}
} // outer while

if (update_flag)
{
vector_handler_->axpy(&ONE, x, vec_Y_, memspace_);
vec_R_->copyDataFrom(rhs, memspace_, memspace_);
matrix_handler_->matvec(A_, vec_Y_, vec_R_, &MINUS_ONE, &ONE, memspace_);
final_residual_norm_ = vector_handler_->dot(vec_R_, vec_R_, memspace_);
final_residual_norm_ = std::sqrt(final_residual_norm_);

// Compare this with initial and update x accordingly
if (final_residual_norm_ < initial_residual_norm_)
{
std::cout << "Update to intial guess is successful, final residual (solution plus update) "
<< std::scientific << std::setprecision(16)
<< final_residual_norm_ / rhsnorm << "\n";

x->copyDataFrom(vec_Y_, memspace_, memspace_);
}
else
{
std::cout << "Update to intial guess is not successful, final residual greater than initial residual "
<< std::scientific << std::setprecision(16)
<< final_residual_norm_ / rhsnorm << "\n";
}
}
else
{
x->copyDataFrom(vec_Y_, memspace_, memspace_);
}

return 0;
}

Expand Down Expand Up @@ -559,6 +619,11 @@ namespace ReSolve
{
vec_V_ = new vector_type(n_, restart_ + 1);
vec_V_->allocate(memspace_);
vec_Y_ = new vector_type(n_);
vec_Y_->allocate(memspace_);
vec_R_ = new vector_type(n_);
vec_R_->allocate(memspace_);

if (flexible_)
{
vec_Z_ = new vector_type(n_, restart_ + 1);
Expand All @@ -585,13 +650,17 @@ namespace ReSolve
delete[] h_rs_;
delete vec_V_;
delete vec_Z_;
delete vec_R_;
delete vec_Y_;

h_H_ = nullptr;
h_c_ = nullptr;
h_s_ = nullptr;
h_rs_ = nullptr;
vec_V_ = nullptr;
vec_Z_ = nullptr;
vec_Y_ = nullptr;
vec_R_ = nullptr;

return 0;
}
Expand Down
2 changes: 2 additions & 0 deletions resolve/LinSolverIterativeFGMRES.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,8 @@ namespace ReSolve

vector_type* vec_V_{nullptr};
vector_type* vec_Z_{nullptr};
vector_type* vec_R_{nullptr};
vector_type* vec_Y_{nullptr};

real_type* h_H_{nullptr};
real_type* h_c_{nullptr};
Expand Down
Loading