3 ms·

Good a place as any to ask I guess. Do any of these GPU libraries have a BiCGStab (or similar) that handles multiple right hand sides? CuPy seems to have GMRES,

Good a place as any to ask I guess. Do any of these GPU libraries have a BiCGStab (or similar) that handles multiple right hand sides? CuPy seems to have GMRES, which would be fine, but as far as I can tell it just does one right hand side.

trostaft 2y ago
IIRC jax's `scipy.sparse.linalg.bicgstab` does support multiple right hand sides. EDIT: Or rather, all the solvers under jax's `scipy.sparse.linalg` all support multiple right hand sides.
- bee_rider 2y ago
  Oh dang, that’s pretty awesome, thanks. “array or tree of arrays” sounds very general, probably even better than an old fashioned 2D array.
  - trostaft 2y ago
    'tree of arrays' Ahh, that's just Jax's concept of pytrees. It was something that they invented to make it easier (this is how I view it, not complete) to pass complex objects to function but still be able to easily consider them as a concatenated vector for AD etc.. E.g. a common pattern is to pass parameters `p` to a function and then internally break them into their physical interpretations, e.g. `mass = p[0]`, `velocity = p[1]`. Pytrees let you just use something like a dictionary `p = {'mass' = 1.0, 'velocity = 1.0'}`, which is a stylistically more natural structure to pass around, and then jax is structured to understand later when AD'ing or otherwise that you're doing so with respect to the 'leaves' of the tree, or the values of the mass and velocity. Hopefully someone corrects me if I'm not right about this. I'm hardly 100% on Jax's vision on PyTrees. As an aside, just a list of right hand sides `[b1, b2, ..., bm]` is valid.
johndough 2y ago
If you have many right hand sides, you could also compute an LU factorization and then solve the right hand sides via back-substitution. https://docs.cupy.dev/en/stable/reference/generated/cupyx.scipy.linalg.lu_factor.html https://docs.cupy.dev/en/stable/reference/generated/cupyx.sc... or https://docs.cupy.dev/en/stable/reference/generated/cupyx.scipy.sparse.linalg.factorized.html#cupyx.scipy.sparse.linalg.factorized https://docs.cupy.dev/en/stable/reference/generated/cupyx.sc... if your linear system is sparse. But whether that works well depends on the problem you are trying to solve.
- bee_rider 2y ago
  My systems are sparse, but might not fit on the GPU when factorized. Actually, usually I do CPU stuff with lots of ram, and Pardiso, so it isn’t an issue. But I was hoping to try out something like ILU+bicgstab on the GPU and the python-verse seems like it has the lowest barrier-to-entry for just playing around.
  - johndough 2y ago
    For my tasks, I had some success with algebraic multigrid solvers as preconditioner, for example from AMGCL or PyAMG. They are also reasonably easy to get started with. https://github.com/pyamg/pyamg https://github.com/pyamg/pyamg https://github.com/ddemidov/amgcl https://github.com/ddemidov/amgcl But I only have to deal with positive definite systems, so YMMV. I am not sure whether those libraries can deal with multiple right-hand sides, but most complexity is in the preconditioners anyway.

trostaft 2y ago

IIRC jax's `scipy.sparse.linalg.bicgstab` does support multiple right hand sides. EDIT: Or rather, all the solvers under jax's `scipy.sparse.linalg` all support multiple right hand sides.

bee_rider 2y ago
Oh dang, that’s pretty awesome, thanks. “array or tree of arrays” sounds very general, probably even better than an old fashioned 2D array.
- trostaft 2y ago
  'tree of arrays' Ahh, that's just Jax's concept of pytrees. It was something that they invented to make it easier (this is how I view it, not complete) to pass complex objects to function but still be able to easily consider them as a concatenated vector for AD etc.. E.g. a common pattern is to pass parameters `p` to a function and then internally break them into their physical interpretations, e.g. `mass = p[0]`, `velocity = p[1]`. Pytrees let you just use something like a dictionary `p = {'mass' = 1.0, 'velocity = 1.0'}`, which is a stylistically more natural structure to pass around, and then jax is structured to understand later when AD'ing or otherwise that you're doing so with respect to the 'leaves' of the tree, or the values of the mass and velocity. Hopefully someone corrects me if I'm not right about this. I'm hardly 100% on Jax's vision on PyTrees. As an aside, just a list of right hand sides `[b1, b2, ..., bm]` is valid.

johndough 2y ago

If you have many right hand sides, you could also compute an LU factorization and then solve the right hand sides via back-substitution. https://docs.cupy.dev/en/stable/reference/generated/cupyx.scipy.linalg.lu_factor.html https://docs.cupy.dev/en/stable/reference/generated/cupyx.sc... or https://docs.cupy.dev/en/stable/reference/generated/cupyx.scipy.sparse.linalg.factorized.html#cupyx.scipy.sparse.linalg.factorized https://docs.cupy.dev/en/stable/reference/generated/cupyx.sc... if your linear system is sparse. But whether that works well depends on the problem you are trying to solve.

bee_rider 2y ago
My systems are sparse, but might not fit on the GPU when factorized. Actually, usually I do CPU stuff with lots of ram, and Pardiso, so it isn’t an issue. But I was hoping to try out something like ILU+bicgstab on the GPU and the python-verse seems like it has the lowest barrier-to-entry for just playing around.
- johndough 2y ago
  For my tasks, I had some success with algebraic multigrid solvers as preconditioner, for example from AMGCL or PyAMG. They are also reasonably easy to get started with. https://github.com/pyamg/pyamg https://github.com/pyamg/pyamg https://github.com/ddemidov/amgcl https://github.com/ddemidov/amgcl But I only have to deal with positive definite systems, so YMMV. I am not sure whether those libraries can deal with multiple right-hand sides, but most complexity is in the preconditioners anyway.