How parallelize the solution of sparse matrices using mldivide
显示 更早的评论
I am trying to parallelize the solution of
x= A\B (mldivide.
My variables are: x = V(I*J*K), A = A(I*J*K x I*J*K) sparse matrix, vec = u (I,J,K) +V/constant +Bswitch (I*J*K x I*J*K)*V
To do this without parallelization, my code currently does this:
V_stacked = reshape(V,I*J*L,1);
vec = u_stacked + V_stacked/Delta + Bswitch*V_stacked;
V_stacked = A\vec;
To parallelize I have tried
u_stacked = reshape(u,I*J,L);
V_stacked = reshape(V,I*J,L);
BswitchTimesVstacked = Bswitch*reshape(V,I*J*L,1);
BswitchTimesVstacked = reshape(BswitchTimesVstacked,I*J,L);
vec = u_stacked + V_stacked/Delta + BswitchTimesVstacked;
tic
parfor l = 1:L
V_stacked(:,l) = A(:,:,l)\vec(:,l);
end
But as A is still I*J*L times I*J*L, it wont work. I am not sure if 1. what I am doing so far is correct and 2. how to reshape B appropriately.
Mathematical info is here: http://www.princeton.edu/~moll/HACTproject/two_asset_kinked.pdf (section 4)
Any help is highly appreciated :-)
4 个评论
Heiko Weichelt
2019-3-1
Hi Emil
In order to give you the most helpful answer, let me ask you a few clarifying questions:
- What is your motivation behind trying to paralellize the sparse MLDIVIDE step?
- What is your expected benefit?
- What sizes will you use (smallest and largest)?
- Do you have a GPU or a cluster available?
More specific to the code you provided:
- Is x = V(I*J*K x I*J*K) and Bswitch (I*J*K x I*J*K) dense or sparse?
- If V is I*J*K x I*J*K and you use V_stacked = reshape(V,I*J*L,1), I assume you meant V_stacked = reshape(V,(I*J*L)^2,1);
In general, notice that sparse MLDIVIDE already uses multi-threading, which means the algorithm itself performs some steps parallel. Usually, as long as the matrix fits in memory and MLDIVIDE's additional memory usage (which can be large compared to the original sparse matrix due to fill-in) can be handled by the avaliavle RAM, it is hard to write a faster algorithm yourself.
Also see https://www.mathworks.com/matlabcentral/answers/95958-which-matlab-functions-benefit-from-multithreaded-computation for more information about multi-threading.
Best,
Heiko
Emil Partsch
2019-3-1
Matt J
2019-3-1
V_stacked(:,l) = A(:,:,l)\vec(:,l);
If A is sparse, how in the above line did you reshape it into a 3D array so as to make it 3D indexable?
Emil Partsch
2019-3-1
编辑：Emil Partsch
2019-3-1
采纳的回答
Matt J
2019-3-4
0 个投票
You can make a 3D stack of sparse matrices A(:,:,l) by converting them to ndSparse type. Then, the parfor construct will work. Alternatively, you can make a block diagonal matrix where all the A(:,:,l) form the diagonal blocks. Then you can solve all systems simultaneously, which would take advantage of Matlab's internal parallelization. I don't know which would be faster.
5 个评论
Emil Partsch
2019-3-4
编辑：Emil Partsch
2019-3-4
Matt J
2019-3-4
No, you missed the point. If you have a block diagonal matrix containing all your equations, there is no need to loop anymore. You can solve the whole system in a single call to mldivide,
V=A\vec(:)
Emil Partsch
2019-3-4
Matt J
2019-3-4
Was it faster when you used parfor? I would not expect so.
Emil Partsch
2019-3-5
编辑：Emil Partsch
2019-3-5
更多回答（1 个）
I recommend using
V_stacked = pagefun(@mldivide,gpuArray(A),gpuArray(vec));
类别
在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
