Classifying sparse matrix in parfor loop

8 ビュー (過去 30 日間)
Bjorn
Bjorn 2021 年 11 月 11 日
コメント済み: Matt J 2021 年 11 月 13 日
Dear all,
The following for-loop is defined:
iend = nz - 1;
jend = nr - 1;
for i = 1 : iend
for j = 1 : jend
k = i * Nr + j;
dz = z_grid(i+1) - z_grid(i);
dr = r_grid(j+1) - r_grid(j);
A_sol(k,k) = - dz * (r_grid(j+1) * a(i+1,j+1) + r_grid(j) * b(i+1,j)) - dr * r_grid(j+1) * (c(i+1,j+1) + d(i,j+1));
end
end
with z_grid and r_grid being vectors and a, b, c and d matrices. The loop takes long to compute and I would like to use parloop to speed up the process.
From MATLAB Help I know that the code should be written in the form of:
iend = nz - 1;
jend = nr - 1;
for i = 1 : iend
for j = 1 : jend
k(j) = i * Nr + j;
dr = r_grid(j+1) - r_grid(j);
end
dz = z_grid(i+1) - z_grid(i);
A_sol(k,k) = -dz* (r_grid(j+1) * a(i+1,j+1) + r_grid(j) * b(i+1,j)) - dr * r_grid(j+1) * (c(i+1,j+1) + d(i,j+1));
end
And write A_sol(k,k) in the form A_sol(k,:). However, the rows and columns have the same indexation, and I cannot come up with a solution to still make this work. Could someone help me out?
Thank you in advance,
WIth kind regards,
Bjorn
  1 件のコメント
Matt J
Matt J 2021 年 11 月 11 日
Shouldn't k be given by,
k = (i-1) * Nr + j;
As you have it now, k will not start at k=1, but instead it will run from Nr+1 to iend*Nr+jend.

サインインしてコメントする。

採用された回答

Matt J
Matt J 2021 年 11 月 11 日
編集済み: Matt J 2021 年 11 月 11 日
The loop takes long to compute and I would like to use parloop to speed up the process.
I doubt it will, but the loop can be further optimized by pre-computing some things. Also, it will be much faster if A_sol is computed using the sparse() command.
r_grid=r_grid(:).'; %make sure these are row vectors
z_grid=z_grid(:).';
DZ=[0,diff(z_grid)].';
DR=[0,diff(r_grid)];
Term1=(-DZ .* a - DR .* c ).*r_grid;
Term2=(-DZ .* b).*r_grid;
Term3=(-DR .* d).*r_grid;
N=length(A_sol);
A_diag=zeros(N,1);
iend = nz - 1;
jend = nr - 1;
for i = 1 : iend
for j = 1 : jend
k = i * Nr + j;
A_diag(k) = Term1(i+1,j+1) +Term2(i+1,j) + Term3(i,j+1);
end
end
A_sol=spdiags(A_diag(:),0,N,N);
  2 件のコメント
Bjorn
Bjorn 2021 年 11 月 13 日
Dear Matt,
Thank you for your promt answer! Defining the matrix after the double for-loop and defining the terms before the double for loop decreased the computational time by up to a factor of ~10 for 2^10 x 2^10 matrices.
Thank you again,
WIth kind regards,
Bjorn
Matt J
Matt J 2021 年 11 月 13 日
Glad it helped but 2^10 is pretty small. You could probably make things faster just by using a non-sparse matrix.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLoops and Conditional Statements についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by