Fast Small Dense Matrix Solver

[ScottBoyce]

I have a general square dense matrix A (not symmetric) which is formed by A=P^TBP 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

[TimP]

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.

[ScottBoyce]

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

[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 - P^TBPx if this can be managed better.

John

[ScottBoyce]

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.

[SergeyKostrov]

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.

[ScottBoyce]

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=P^TBP, but I am unsure if there is anything faster than using DGEMM.

[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=P^TBP 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