-
Notifications
You must be signed in to change notification settings - Fork 89
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.
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
Iterative refinement (IR) #243
Conversation
|
||
solver_->apply(lend(one_op), lend(residual), lend(one_op), dense_x); | ||
residual->copy_from(dense_b); | ||
system_matrix_->apply(lend(neg_one_op), dense_x, lend(one_op), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's quite simple, but nonetheless I believe we in general put as comments the general algorithm of the solver, maybe you could also do that here? (or was that only for the step kernels?)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
it was only for the step kernels
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
LGTM
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
LGTM!
size_type num_cols, stopping_status *stop_status) | ||
{ | ||
const auto tidx = | ||
static_cast<size_type>(blockDim.x) * blockIdx.x + threadIdx.x; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is the first time I see it handled that way (using explicit casting to size_type
), what is the normal auto
deduction? Is it not uint64
?
I did a quick implementation of the iterative refinement method after getting back from work today (yesterday for you). It seemed like to much of a low hanging fruit not to have it in the first release.
Needles, to say, this does not try to do mixed precision IR (MPIR), but just the plain old IR, where you can specify what you want to use as the inner solver. By default I set it to the simplest possible thing - the identity operator, which results in Richardson iteration with scaling parameter 1.
To add mixed precision, we would have to solve the problem of conversions between value types - for that we would need to implement #247 - but even without it it's still useful, as it gives us Richardson iteration support, and if the default solver is replaced by the block-Jacobi precondtiioner, we actually get the block-Jacobi relaxation method (or what would be the "adaptive block-Jacobi relaxation method").
While writing the unit tests, I also figured out that
stopping_status
did not have comparisons, so I added them, and while adding them I noticed that we had 2 methods for resetting the status (clear()
andreset()
), so I removedclear()
and usedreset()
whereverclear()
was used previously.