topic Iterative solvers are most in Intel® Fortran Compiler

Fast Small Dense Matrix Solver

ScottBoyce — Wed, 06 Mar 2013 21:29:55 GMT

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

Iterative solvers are most

TimP — Thu, 07 Mar 2013 01:35:25 GMT

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.

I am linking to the math

ScottBoyce — Thu, 07 Mar 2013 02:21:37 GMT

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

Scott,

John_Campbell — Thu, 07 Mar 2013 04:13:57 GMT

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

Yes, B is non-symmetric and P

ScottBoyce — Thu, 07 Mar 2013 04:19:05 GMT

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.

I have a generic question.

SergeyKostrov — Thu, 07 Mar 2013 14:11:23 GMT

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.

It only takes a few seconds,

ScottBoyce — Thu, 07 Mar 2013 15:58:34 GMT

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.

Scott,

John_Campbell — Mon, 11 Mar 2013 05:33:37 GMT

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