Skip to content
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

Implement BiCGSTAB #340

Merged
merged 1 commit into from
Oct 8, 2023
Merged

Implement BiCGSTAB #340

merged 1 commit into from
Oct 8, 2023

Conversation

jlogan03
Copy link
Collaborator

@jlogan03 jlogan03 commented Oct 8, 2023

Some thoughts about this implementation -

Pros

  • It runs and is even fairly correct and a bit documented
  • Implemented on f32 and f64 with generic indices

Cons

  • Doesn't represent any particular reusable abstract pattern of iterative solver interface, and only makes a little progress toward the goals of Linear solver interfaces #339
  • Solver call signature requires the turbof***
    • I haven't found any way around this - even when it's obvious from the inputs and outputs which variant should be selected, it still fails to select the correct variant without explicit type annotation
  • Macros make things worse for the users in terms of readability, understanding errors, etc
  • Not fully generic on value type
  • Not generic on sparse vs. dense RHS and initial guess
  • Can't be implemented directly on Complex types
  • Allocates then drops several vectors on each iteration which have the same sparsity at every iteration and could reuse the same memory if we had a mechanism for this
    • In practice, this isn't as bad as it sounds since the first few iterations will clear out a correctly sized region of heap s.t. the subsequent allocations are cheap to perform... as long as there is not much else happening on the machine

I spent a solid 8 hours attempting to make a version that is generic over the value type. This did not go well.

Some of the implementations of the underlying matrix-vector and vector-scalar operations are not fully generic, so the trait bounds on the impl end up being more than 10 lines long and incredibly unclear as to which types will and will not be compatible, and the type system ultimately fails on overflow due to inferring generic outputs across iterations in the solver loop.

This could all be fixed to some extent, but would require quite a bit of work to go back into the core of the crate and make a single definition of what trait bounds constitute a mathematical field over which we can do linear algebra, then propagate that definition to every struct and impl in the crate, which I believe is out of scope for this PR.

@jlogan03 jlogan03 marked this pull request as ready for review October 8, 2023 18:19
jlogan03 added a commit to jlogan03/sprs that referenced this pull request Oct 8, 2023
@jlogan03 jlogan03 changed the title BiCGSTAB solver Implement BiCGSTAB Oct 8, 2023
Copy link
Collaborator

@mulimoen mulimoen left a comment

Choose a reason for hiding this comment

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

This is some fantastic work and it's great to see a proper general solver. This will be immensely useful to users.

Making everything generic has many drawbacks, the ergonomics can get much worse and it takes a lot longer to write for dubious benefit, the macro is great for implementing the two most useful impls and we can deal with Complex if anybody asks for it.

I only have some minor nits on this PR, I consider it good for a merge in the current state

sprs/src/sparse/linalg/bicgstab.rs Outdated Show resolved Hide resolved
sprs/src/sparse/linalg/bicgstab.rs Show resolved Hide resolved
sprs/src/sparse/linalg/bicgstab.rs Show resolved Hide resolved
sprs/src/sparse/linalg/bicgstab.rs Outdated Show resolved Hide resolved
sprs/src/sparse/linalg/bicgstab.rs Outdated Show resolved Hide resolved
sprs/src/sparse/linalg/bicgstab.rs Show resolved Hide resolved
@jlogan03
Copy link
Collaborator Author

jlogan03 commented Oct 8, 2023

@mulimoen I like the approx usage for testing, but had to roll it back because we don't keep approx as a default feature so it crashes the tests. I'd be happy to bring it back if we want to update either default features to include approx, or to include the approx feature for testing

@mulimoen
Copy link
Collaborator

mulimoen commented Oct 8, 2023

I don't feel strongly about approx, bit of a pain to enable it for doctests. A bit of resolver history: The resolver previously merged features across build types, so having approx on during tests would enable it also for regular builds. Since we are using resolver v2 we can now use it for tests only

@jlogan03 jlogan03 merged commit d784194 into sparsemat:master Oct 8, 2023
18 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants