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

Jacobi, Gauss-Seidel & SOR iteration for SparseMatrixCSC using iterators #156

Merged
merged 5 commits into from Aug 11, 2017

Conversation

Projects
None yet
2 participants
@haampie
Collaborator

haampie commented Jul 31, 2017

Unfortunately the simple stationary methods like Gauss-Seidel and SOR are a bit awkward to implement for column-major (sparse) matrices.

My goal with this PR is to support stationary methods without copying or transforming the underlying SparseMatrixCSC, while still exploiting in-place updates as much as possible.

Todo:

  • Tests
  • Combine this with current implementation for dense matrices
  • Symmetric SOR
  • Sensible stopping criterion*

* It should be avoided to compute norm(b - Ax) since that's as expensive as an iteration. Another option is ||x_next - x_curr|| / ||x_cur|| < tol, but that might be complicated in the current Gauss-Seidel implementation as it updates x in-place. At the moment the stopping criterion is just the maximum number of iterations. That might be good enough as well: stationary methods are often used as smoothers in multigrid methods where just a small, fixed number of iterations is carried out.

Fixes #147

@haampie

This comment has been minimized.

Show comment
Hide comment
@haampie

haampie Aug 9, 2017

Collaborator

At this point the sparse version on full matrices is significantly faster... :p

A = rand(30, 30) + 30I
b_sparse = @benchmark IterativeSolvers.jacobi!(x, mat, b, maxiter = 10) setup=(x=zeros(30); mat=sparse($A); b=$A*ones(30))
b_full = @benchmark IterativeSolvers.jacobi!(x, mat, b, maxiter = 10) setup=(x=zeros(30); mat=$A; b=$A*ones(30))
> b_sparse
BenchmarkTools.Trial:
  memory estimate:  1.05 KiB
  allocs estimate:  9
  minimum time:     11.584 μs (0.00% GC)
  samples:          10000
> b_full
BenchmarkTools.Trial:
  memory estimate:  126.45 KiB
  allocs estimate:  957
  minimum time:     71.540 μs (0.00% GC)
  samples:          10000

The current / full matrix version computes norm(A*x-b) which is allocating and as expensive as an iteration. But it might also have to do with row-wise matrix access.

Collaborator

haampie commented Aug 9, 2017

At this point the sparse version on full matrices is significantly faster... :p

A = rand(30, 30) + 30I
b_sparse = @benchmark IterativeSolvers.jacobi!(x, mat, b, maxiter = 10) setup=(x=zeros(30); mat=sparse($A); b=$A*ones(30))
b_full = @benchmark IterativeSolvers.jacobi!(x, mat, b, maxiter = 10) setup=(x=zeros(30); mat=$A; b=$A*ones(30))
> b_sparse
BenchmarkTools.Trial:
  memory estimate:  1.05 KiB
  allocs estimate:  9
  minimum time:     11.584 μs (0.00% GC)
  samples:          10000
> b_full
BenchmarkTools.Trial:
  memory estimate:  126.45 KiB
  allocs estimate:  957
  minimum time:     71.540 μs (0.00% GC)
  samples:          10000

The current / full matrix version computes norm(A*x-b) which is allocating and as expensive as an iteration. But it might also have to do with row-wise matrix access.

@haampie

This comment has been minimized.

Show comment
Hide comment
@haampie

haampie Aug 9, 2017

Collaborator

For dense matrices it makes sense to use the same column-wise approach as in the sparse implementation. See https://gist.github.com/haampie/6d24c91df406e9574f3d4cd451db2260

For matrices of size 300x300 you get a 2x speedup (Jacobi):

improvement

Or about 1.5x in SOR:

sor

Collaborator

haampie commented Aug 9, 2017

For dense matrices it makes sense to use the same column-wise approach as in the sparse implementation. See https://gist.github.com/haampie/6d24c91df406e9574f3d4cd451db2260

For matrices of size 300x300 you get a 2x speedup (Jacobi):

improvement

Or about 1.5x in SOR:

sor

@haampie haampie changed the title from WIP: Jacobi, Gauss-Seidel & SOR iteration for SparseMatrixCSC using iterators to Jacobi, Gauss-Seidel & SOR iteration for SparseMatrixCSC using iterators Aug 11, 2017

@andreasnoack andreasnoack merged commit c2a6698 into JuliaMath:master Aug 11, 2017

2 checks passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.6%) to 88.457%
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment