Translated by

このページのコンテンツは英語から自動翻訳されています。自動翻訳をオフにする場合は「<a class="turn_off_mt" href="#" onclick="window._kiq.push(['set', { 'event': 'Turn off MT' }]);">ここ</a>」をクリックしてください。

Solving A{k} * x + b = 0 for large numbers of A{k} with same structure/filling.

Cedric Wannaz

Cedric Wannaz (view profile)

さんによって質問されました 2017 年 12 月 14 日

Cedric Wannaz (view profile)

さんによって 編集されました 2017 年 12 月 25 日
Christine Tobler

Christine Tobler (view profile)

さんの 回答が採用されました
Dear all,
A{k} are square, non-symmetrical, large ( NxN with N in 1E5-5E5) sparse matrices (typical density <1E-4), column strictly diagonally dominant. They all have the same (block, see EDIT 1) structure/filling :
and b is Nx1 with nnz(b) typically in 1-1E4.
---- Digression start ----
I discussed the case A*x+B=0 already in another context, for B NxM with M large, where I can take advantage of a single LU factorization of a single A and solve in parallel:
[L, U, P, Q, R] = lu( A ) ;
parfor ...
X{k} = Q * (U \ (L \ (P * (R \ B_cell{k}))))
end
where B_cell splits B in blocks of appropriate size, in this thread.
---- Digression end ----
My question in the current context is: is there a transformation/factorization that I could use, that exploits the fact that all A{k} have strictly the same filling/structure (i.e. non-zero elements are always at the same place), that I would apply e.g. to A{1} and then re-use to accelerate the solving of A{k}*x+b=0 for all other A{k} (e.g. by providing MLDIVIDE (or LINSOLVE with smaller, dense blocks) with matrices that shortcut part of the factorization)?
Also, could the block-structure and the partial filling of b be used in any manner?
EDIT 1, 2017/12/19 @ 17:32 UTC - On this line, if I perform a LU decomposition as shown in the digression above, for example, I get matrices P (row permutation) and Q (column reordering) that will (or should) be the same for all A{k}. Noting that:
[L, U, P, Q, R] = lu( A{k} ) ;
x{k} = Q * (U \ (L \ (P * (R \ b)))) ;
is always more efficient than:
x{k} = A{k} \ b ;
because we avoid the analysis of A, is there a way to exploit the fact that P and Q are known after the first iteration (i.e. P and Q computed for k=1 will remain the same for all ks)?
Typical A{k} and b can be downloaded here as a MAT File.
Cheers,
Cedric
EDIT 1, 2017/12/14 @ 19:39 UTC
Here is another sparse matrix where I delineated the block structure. It is associated with another simulation, hence the different size, but in this simulation all A{k} are 99149x99149 with the same structure:

0 件のコメント

サインイン to comment.

2 件の回答

Christine Tobler (view profile)

2017 年 12 月 19 日
採用された回答

In Tim Davis' C++ library UMFPACK, there are two phases to computing the decomposition of a matrix: The first one only uses the sparsity structure, while the second one also uses the numbers. Unfortunately, such separate phases are not accessible through MATLAB directly.
Another approach, if the numbers in A{k} do not change too strongly, would be to use the LU factorization of A{1} as a preconditioner for an iterative solver (for example gmres, pcg).

Steven Lord

Steven Lord (view profile)

2017 年 12 月 19 日
And if you expect that the solution for a particular A{k} should be "close" to the solution for the next A{k+1}, most if not all of the iterative solvers let you specify a starting point for the iterative method. Specifying the previous solution may help speed the next iteration. Check the documentation for the iterative method of your choice for the x0 input argument.
Cedric Wannaz

Cedric Wannaz (view profile)

2017 年 12 月 19 日
Dear Christine and Steven, thank you very much for your inputs. This is exactly what I was after, as I was stalling and this gives me a series of other steps/approaches to test.
Cheers, Cedric

サインイン to comment.

Bjorn Gustavsson (view profile)

2017 年 12 月 14 日

While we are waiting for good replies I can only suggest that you look at Tim Davis' contributions at the file exchange (in the case you haven't done so yet): FACTORIZE, sparseinv. Perhaps he has implemented something serving your needs.
HTH

Cedric Wannaz

Cedric Wannaz (view profile)

2017 年 12 月 14 日
Thank you Bjorn! The factorization that I mention for the other context actually comes from Tim Davis' contribution (I went through suitesparse and then through his FEX posts). I have not found any approach yet based on this material for the current question though.

サインイン to comment.

Translated by