Fast Small Dense Matrix Solver

Fast Small Dense Matrix Solver

Portrait de Scott B.

I have a general square dense matrix A (not symmetric) which is formed by A=PTBP where B was in a compressed storage scheme and P is a rectangular matrix. The size of A ranges from 10x10 to 500x500, where B can be 150,000x150,000 and is sparse.

What would be the best way to solve for x given b (system of linear equations) that result from

Ax=b  =>  x=A-1b

Right now I am just using LAPACK DGESV that is linked to MKL (so assume I am using their solver). Is there any benifit to going to a interative solver or any recomendations as to how to best solve this system of equations as fast as possible.

Thanks for any comments

8 posts / 0 nouveau(x)
Dernière contribution
Reportez-vous à notre Notice d'optimisation pour plus d'informations sur les choix et l'optimisation des performances dans les produits logiciels Intel.
Portrait de Tim Prince

Iterative solvers are most likely to be used where there is high sparsity within the band so that banded or full matrix solvers are inefficient.   Fully 3 dimensional field problems (finite element etc.) are likely to produce such sparse structure, while 2D or nearly 2D problems may be well handled by banded or skyline storage direct solvers.  As you are considering MKL, this is probably more topical on that forum.

Portrait de Scott B.

I am linking to the math kernel library and solving the system of equations by calling DGESV() .

Portrait de John Campbell

Scott,

As A is dense and non-symmetric ( B non-symmetric also?), I would expect that a direct, rather than itterative solver would be more robust. You might also want to check how "well-conditioned" the solution is. The solution might require a more general direct solver with row pivoting if this becomes an issue. I would calculate b - A.x to check the solution accuracy, or even b - PTBPx if this can be managed better.

John

 

Portrait de Scott B.

Yes, B is non-symmetric and P is orthonormal rectangular matrix. The solution is very stable with no near zero eigenvalues. LAPACK DGESV solves Ax=b through LU decomposition with partial pivoting and row interchanges.

Portrait de Sergey Kostrov

I have a generic question.

>>...The size of A ... 500x500, where B can be 150,000x150,000...

How long does it take to solve it on your computer? Thanks in advance.

Portrait de Scott B.

It only takes a few seconds, but for each solution of A creates a new version of B and which is then matrix multiplied by P to build a new version of A which then needs a new solution. I like to speed up, even by a fraction of a second, solving the system of equations. There also is of course a slow down do to the A=PTBP, but I am unsure if there is anything faster than using DGEMM.

Portrait de John Campbell

Scott,

You should record the times for the different stages of the calculation to determine where the effort should be placed.
If the calculation :  A=PTBP is a significant time, then utilising the sparsity of B might be significant.
In a large finite element formation, the calculation of force in f = K.x can be carried out much quicker if K is considered as f = Sum (K_element.x) when the element K matrices are easily (quickly) available.

John

Connectez-vous pour laisser un commentaire.