{"action":"created","issue":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/19500","repository_url":"https://api.github.com/repos/JuliaLang/julia","labels_url":"https://api.github.com/repos/JuliaLang/julia/issues/19500/labels{/name}","comments_url":"https://api.github.com/repos/JuliaLang/julia/issues/19500/comments","events_url":"https://api.github.com/repos/JuliaLang/julia/issues/19500/events","html_url":"https://github.com/JuliaLang/julia/issues/19500","id":193313167,"number":19500,"title":"Speed up UMFpack LU for multiple (sparse) right hand sides","user":{"login":"lruthotto","id":4249607,"avatar_url":"https://avatars.githubusercontent.com/u/4249607?v=3","gravatar_id":"","url":"https://api.github.com/users/lruthotto","html_url":"https://github.com/lruthotto","followers_url":"https://api.github.com/users/lruthotto/followers","following_url":"https://api.github.com/users/lruthotto/following{/other_user}","gists_url":"https://api.github.com/users/lruthotto/gists{/gist_id}","starred_url":"https://api.github.com/users/lruthotto/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/lruthotto/subscriptions","organizations_url":"https://api.github.com/users/lruthotto/orgs","repos_url":"https://api.github.com/users/lruthotto/repos","events_url":"https://api.github.com/users/lruthotto/events{/privacy}","received_events_url":"https://api.github.com/users/lruthotto/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":5,"created_at":"2016-12-03T20:59:34Z","updated_at":"2016-12-04T15:51:30Z","closed_at":null,"body":"The LU factorization (through UMFpack) currently implemented in julia 0.5 is rather slow when using it for multiple right hand sides. Looking at umfpack.jl lines 259 you can see that there is a loop over all right hand sides. \r\n\r\nAs far as I could see UMFpack does not seem to have a built-in option for multiple rhs. So, i followed the MATLAB code [here](https://github.com/dealii/dealii/blob/master/bundled/umfpack/UMFPACK/MATLAB/umfpack_solve.m) to build one. \r\n\r\nIt works quite nicely and I'd be pleased to contribute that to Base. To this end, it would be great to get some help and advice. For example, I'm not sure where to put the code. Also, I'm not sure why it takes so long to get the parts of the lufactorization.\r\n\r\nSee this example here: \r\n\r\n```jl\r\nfunction myLUsolve(A::Base.SparseArrays.UMFPACK.UmfpackLU,B)\r\n\r\n    # get parts of LUfactorization\r\n    println(\"time to extract LUfactors\")\r\n    @time (L,U,p,q,R) = A[:(:)];\r\n\r\n    println(\"time to solve for $(size(B,2)) right hand sides\")\r\n    @time begin\r\n        X = zeros(eltype(B),size(B))\r\n        X[q,:] = U\\(L\\((Diagonal(R) * B)[p,:]))\r\n    end\r\n    return X\r\nend\r\n\r\nnrhs = 300\r\nn    = 10000\r\nA    = UniformScaling(10) + sprandn(n,n,0.0005);\r\nrhs  = randn(n,nrhs);\r\n\r\nprintln(\"=== lufact(A)\\\\rhs for  real matrix ===\")\r\nprintln(\"\\t factorize\")\r\n@time luA = lufact(A)\r\nprintln(\"\\t julia solve\")\r\n@time t1 =  luA\\rhs\r\nprintln(\"\\t my solve\")\r\n@time t2 = myLUsolve(luA,rhs)\r\nprintln(\"err: $(norm(t1-t2)/norm(t1))\")\r\n```\r\nOn my machine, Ubuntu 14, julia 0.5.1, I get the following timings\r\n```jl\r\n=== lufact(A)\\rhs for  real matrix ===\r\n\t factorize\r\n  5.479585 seconds (66 allocations: 536.308 MB, 0.64% gc time)\r\n\t julia solve\r\n 29.411826 seconds (613 allocations: 183.115 MB, 0.82% gc time)\r\n\t my solve\r\ntime to extract LUfactors\r\n  1.517049 seconds (42 allocations: 562.236 MB, 5.83% gc time)\r\ntime to solve for 300 right hand sides\r\n 14.070643 seconds (31 allocations: 114.442 MB, 0.96% gc time)\r\n 15.619639 seconds (7.67 k allocations: 677.021 MB, 1.44% gc time)\r\nerr: 1.5411555266919556e-12\r\n```\r\n\r\n"},"comment":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/comments/264711512","html_url":"https://github.com/JuliaLang/julia/issues/19500#issuecomment-264711512","issue_url":"https://api.github.com/repos/JuliaLang/julia/issues/19500","id":264711512,"user":{"login":"lruthotto","id":4249607,"avatar_url":"https://avatars.githubusercontent.com/u/4249607?v=3","gravatar_id":"","url":"https://api.github.com/users/lruthotto","html_url":"https://github.com/lruthotto","followers_url":"https://api.github.com/users/lruthotto/followers","following_url":"https://api.github.com/users/lruthotto/following{/other_user}","gists_url":"https://api.github.com/users/lruthotto/gists{/gist_id}","starred_url":"https://api.github.com/users/lruthotto/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/lruthotto/subscriptions","organizations_url":"https://api.github.com/users/lruthotto/orgs","repos_url":"https://api.github.com/users/lruthotto/repos","events_url":"https://api.github.com/users/lruthotto/events{/privacy}","received_events_url":"https://api.github.com/users/lruthotto/received_events","type":"User","site_admin":false},"created_at":"2016-12-04T15:51:30Z","updated_at":"2016-12-04T15:51:30Z","body":"Thanks for the hint. I see a similar speedup now when deactivating iterative refinements. I do see additional room for speedup here (or elsewhere) though if we could have an OpenMP type parallelism for the upper- and lower-triangular solves. UMFpack seems to be using only one thread on my machine on the solve stage. "}}
{"actor":{"display_login":"lruthotto"}}