decomposition used in parfor

13 次查看（过去 30 天）
Michal
Michal 2022-8-18
评论： Matt J 2022-8-19
In the following problem, where I need to solve many times same linear system "A" with different right sides "b" is very benefitial to use decomposition, which in this case works typically much faster, than simple A\b.
See this example:
clear all
N = 300;
M = 2000;
A = rand(N);
b = rand(N,M);
x = zeros(N,M);
tic;
for ii=1:M
x(:,ii) = A \ b(:,ii);
end
toc
tic;
dA = decomposition(A);
dAc = parallel.pool.Constant(dA);
parfor ii=1:M
x(:,ii) = dAc.Value \ b(:,ii);
end
toc
As I learned from here, TMW expert proposed to use parallel.pool.Constant, but this does not work at all!!!
Warning: Saving a decomposition is not supported.
Is there any way how to use decomposition together with parfor or any other parallel evaluation construct???
  7 个评论
显示 5更早的评论
Edric Ellis
Edric Ellis 2022-8-19
Ran in:
Ah OK, I should have been clearer. This case needs the function_handle constructor for parallel.pool.Constant.
N = 10;
M = 4;
A = rand(N);
b = rand(N,M);
x = zeros(N,M);
dAc = parallel.pool.Constant(@() decomposition(A));
parfor ii=1:M
x(:,ii) = dAc.Value \ b(:,ii);
end
disp('done.')
done.
This causes each worker to build the decomposition for itself from the value of A.
Bruno Luong
Bruno Luong 2022-8-19
So this work around makes the same decomposition performed by all worker instead of once, right?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2022-8-18
编辑：Matt J 2022-8-18
Your example doesn't describe a situation where either decomposition() or parfor will be beneficial. You should instead just use x=A\b.
N = 300;
M = 2000;
A = rand(N);
b = rand(N,M);
x = zeros(N,M);
tic;
x=A\b;
toc%Elapsed time is 0.010876 seconds.
[L,U,P]=lu(A);
tic;
parfor ii=1:M
x(:,ii) = U\(L\(P *b(:,ii)));
end
toc%Elapsed time is 0.284230 seconds.
If the problem is that the b(:,ii) are not simultaneously available, you can do the LU decomposition explicitly, as in my example above.
  15 个评论
显示 13更早的评论
Bruno Luong
Bruno Luong 2022-8-18
I have another take, It is very simple
each column iA(:,j) where iA designate inv(A) is solution of
A*xj = ej
at the numerical precision error.
If you believe the algorithm behind inv is a "bad" one, the simply compute
iA = A \ eye(size(A)).
Personaly I use INV many many time and I nevr encount precision issue more than "\".
If you don't believe that the solution of
A*x = b
can be approximate by
x = sum(b(j)*xj)
from the most fundamental property required by the linear operator, then you don't believe the whole algebra computation with finite precision, or the problem is elsewhere (your system is ill-conditioning to start with).
Matt J
Matt J 2022-8-19
This is my primary problem: I would like to use precomputed dA = decomposition(A) at parfor loop, where each b(:,ii) -> b_ii is computed independently, but this approach is impossible due to the fact, that parfor does not work in this case.
If the b are computed independently, it would be better to postpone the inversion until after the parfor loop
parfor ii=1:M
b_ii = rand(N,1); % independent generation of the ii-th rhs b
B(:,ii) = b_ii;
end
x=A\B;

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by