Are there any matrix-free solvers?

[medakk]

Hi! Does ndarray support any matrix-free solvers?

In matrix-free solver, instead of x = solve(A, b), we run x = solve(f, b), where f is a function f(y) = Ay that can evaluate the product of A and any vector. This is useful in some situations where we don't need to store the entire matrix. https://en.wikipedia.org/wiki/Matrix-free_methods

This is Eigen's version: https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html

[bytesnake]

For eigenvalue problems there is #184 in the pipeline. It implements a blocked conjugated gradient algorithms in matrix-free fashion by accepting the linear operator as a closure here. A similar conjugated gradient algorithm could be implemented for equation solver.

[bytesnake]

Do you have a specific use-case?

[medakk]

I'm working on a physics simulator. During the time integration step, I have to solve for a large matrix A Ax=b. The matrix A is quite sparse however and I could just use the sparse solver. But most of A is calculating the forces and the gradients of the forces, which could be done with explicitly storing it.

[bytesnake]

You could do the following for symmetric positive definite matrix: * Use cholesky decomposition from `sprs` crate and solve manually * Implement conjugated gradient, which is quite simple. You could search for `painless conjugated gradient` (they even have canned versions at the end which can be directly implemented without understanding the details) Otherwise you need to implement GMRES (see issue #107) or BiCGSTAB. 13.04.2020 20:23:43 Karthik Karanth <notifications@github.com>:

…

I'm working on a physics simulator. During the time integration step, I have to solve for a large matrix A Ax=b . The matrix A is quite sparse however and I could just use the sparse solver. But most of A is calculating the forces and the gradients of the forces, which could be done with explicitly storing it. — You are receiving this because you commented. Reply to this email directly, view it on GitHub [#190 (comment)] , or unsubscribe [https://github.com/notifications/unsubscribe-auth/AAHRRKP3LZDLXU67R6SA4ELRMNKCFANCNFSM4MGDI4YQ] . [https://github.com/notifications/beacon/AAHRRKLVFRVYRLKLB327J5DRMNKCFA5CNFSM4MGDI4Y2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOESFAXCI.gif]

[bytesnake]

The first version is obviously not matrix free 13.04.2020 20:23:43 Karthik Karanth <notifications@github.com>:

…

I'm working on a physics simulator. During the time integration step, I have to solve for a large matrix A Ax=b . The matrix A is quite sparse however and I could just use the sparse solver. But most of A is calculating the forces and the gradients of the forces, which could be done with explicitly storing it. — You are receiving this because you commented. Reply to this email directly, view it on GitHub [#190 (comment)] , or unsubscribe [https://github.com/notifications/unsubscribe-auth/AAHRRKP3LZDLXU67R6SA4ELRMNKCFANCNFSM4MGDI4YQ] . [https://github.com/notifications/beacon/AAHRRKLVFRVYRLKLB327J5DRMNKCFA5CNFSM4MGDI4Y2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOESFAXCI.gif]

[medakk]

Thanks! I was reading that exact document(painless conjugate gradients). I think I'll use sprs for my work now.

Otherwise you need to implement GMRES (see issue #107) or BiCGSTAB.

By implement, do you mean a pure rust implementation or calling into a BLAS library to solve it? If you could give me some pointers on which part of the codebase this would touch, I could work on it.

[bytesnake]

You still want to call into a BLAS library for optimized matrix-vector product, but the implementation itself is in Rust. If you're matrix is not symmetric, then the BiCGSTAB algorithm is IMO more easy to implement. You could also try to adapt previous attempts (these are two years old and probably won't compile anymore)
https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method
https://users.rust-lang.org/t/prior-work-on-krylov-subspace-methods/13261/5

[medakk]

EDIT: Disregard, I hadn't read the docs completely.