Direct Solvers for Sparse Matrices

by Pavel Holoborodko on September 5, 2013

Newest version of toolbox includes sparse linear system solver – function mldivide or operator "\".

Solver analyzes input matrix and automatically uses the most suitable decomposition:

  1. Cholesky LDLT factorization is applied if input matrix is symmetrical/hermitian and has real positive diagonal.
  2. LU factorization is applied if LDLT failed or directly if matrix is obviously non-SPD.
    (Algorithm is based on ideas from SuperLU and uses column approximate minimum degree permutation)
  3. QR factorization for rectangular matrices to handle least squares problems (beta).
    (Left-looking QR decomposition, also uses colamd pre-ordering)

All methods are capable to compute with arbitrary precision, in real and complex cases. Tables below present timing for solving sparse linear systems with quadruple precision (34 decimal digits).

Real symmetric positive definite problems (SPD) – LDLT.
Matrix NameSizeNNZTime (sec)Relative Residual
nos1 237 x 237 1017 0.00441.55701e-38
nos2 957 x 957 4137 0.06515.2701e-38
nos3 960 x 960 15844 0.10381.12317e-33
nos4 100 x 100 594 0.00122.87428e-32
nos5 468 x 468 5172 0.04857.46582e-38
nos6 675 x 675 3255 0.03617.80187e-36
nos7 729 x 729 4617 0.06987.43915e-34
gr_30_30 900 x 900 7744 0.07158.77471e-34
bcsstk27 1224 x 1224 56126 0.21322.86669e-39
bcsstk01 48 x 48 400 0.00069.20784e-42
bcsstk02 66 x 66 4356 0.00415.87268e-36
bcsstk03 112 x 112 640 0.00121.46672e-42
bcsstk04 132 x 132 3648 0.00531.53244e-38
bcsstk05 153 x 153 2423 0.00391.1895e-38
bcsstk06 420 x 420 7860 0.02501.62052e-40
bcsstk07 420 x 420 7860 0.02481.63422e-40
bcsstk08 1074 x 1074 12960 0.13074.39824e-43
bcsstk09 1083 x 1083 18437 0.24463.80098e-39
bcsstk10 1086 x 1086 22070 0.10764.28618e-39
bcsstk11 1473 x 1473 34241 0.23581.08379e-38
bcsstk12 1473 x 1473 34241 0.23588.39512e-39
bcsstk13 2003 x 2003 83883 1.72034.76763e-42
bcsstk14 1806 x 1806 63454 0.51141.58546e-41
bcsstk15 3948 x 3948117816 5.43729.1493e-40
bcsstk16 4884 x 4884290378 6.79366.25152e-42
bcsstk1710974 x 10974428650 13.58955.52591e-40
bcsstk1811948 x 11948149090 13.82491.56103e-41
s1rmq4m1 5489 x 5489262411 7.09112.70504e-35
s2rmq4m1 5489 x 5489263351 7.90201.73581e-32
s3rmq4m1 5489 x 5489262943 7.31261.66779e-29
s1rmt3m1 5489 x 5489217651 4.50472.01389e-35
s2rmt3m1 5489 x 5489217681 4.41691.20526e-32
s3rmt3m1 5489 x 5489217669 4.31221.26389e-29
s3rmt3m3 5357 x 5357207123 3.74886.17244e-30
Real unsymmetric problems – LU.
Matrix NameSizeNNZTime (sec)Relative Residual
west0067 67 x 67 294 0.0009 4.25499e-35
west0132 132 x 132 414 0.0009 1.04691e-36
west0156 156 x 156 371 0.0010 2.90831e-21
west0167 167 x 167 507 0.0012 1.3336e-36
west0381 381 x 381 2157 0.0217 2.66493e-36
west0479 479 x 479 1910 0.0093 3.55775e-36
west0497 497 x 497 1727 0.0068 6.09013e-37
west0655 655 x 655 2854 0.0163 2.85649e-36
west0989 989 x 989 3537 0.0229 7.21709e-36
west15051505 x 1505 5445 0.0508 3.435e-36
west20212021 x 2021 7353 0.0886 3.88824e-36
bp___200 822 x 822 3802 0.0234 2.60261e-35
bp___400 822 x 822 4028 0.0307 1.58826e-34
bp___600 822 x 822 4172 0.0356 2.94056e-35
bp___800 822 x 822 4534 0.0407 1.18475e-35
bp__1000 822 x 822 4661 0.0440 2.48611e-34
bp__1200 822 x 822 4726 0.0360 4.46649e-34
bp__1400 822 x 822 4790 0.0418 2.69099e-35
bp__1600 822 x 822 4841 0.0313 2.1945e-35
impcol_a 207 x 207 572 0.0015 8.58996e-34
impcol_b 59 x 59 312 0.0005 7.71085e-34
impcol_c 137 x 137 411 0.0007 2.15967e-36
impcol_d 425 x 425 1339 0.0050 5.7882e-35
impcol_e 225 x 225 1308 0.0021 2.25724e-37
mcca 180 x 180 2659 0.0066 3.47919e-48
mcfe 765 x 765 24382 0.1163 2.99781e-48
nnc261 261 x 261 1500 0.0069 1.10753e-28
nnc666 666 x 666 4044 0.0468 1.27587e-31
nnc13741374 x 1374 8606 0.2009 3.83239e-28
orsirr_2 886 x 886 5970 0.2107 6.72416e-37
orsirr_11030 x 1030 6858 0.3204 7.12824e-37
orsreg_12205 x 2205 14133 1.4409 1.08437e-35
watt__11856 x 1856 11360 0.6724 8.83928e-32
watt__11856 x 1856 11360 0.6723 5.5764e-33
lns_39373937 x 3937 25407 1.3501 1.10669e-35
lnsp39373937 x 3937 25407 1.2554 2.56864e-35
sherman11000 x 1000 3750 0.0686 6.54415e-33
sherman21080 x 1080 23094 1.6683 2.93633e-37
sherman35005 x 5005 20033 2.0002 4.00437e-38
sherman41104 x 1104 3786 0.0583 4.35491e-34
sherman53312 x 3312 20793 0.8152 8.54008e-35
e05r0500 236 x 236 5856 0.0277 2.4959e-33
e20r05004241 x 4241131556 5.1604 2.15125e-31
Overdetermined systems (least squares problem) – QR.
Matrix NameSizeNNZTime (sec)Relative Residual
illc10331033 x 32047320.17916.49442e-06
illc18501850 x 71287584.28471.04108e-05
well10331033 x 32047320.17897.0929e-06
well18501850 x 71287584.21851.12967e-05

The tests were run by André Van Moer on his system: Core i7-3960X Extreme Edition, 3900 MHz, 64GB of system memory with Windows 7 installed. (André, thank you for your constant support, suggestions and enthusiastic help!)

Test matrices are from MatrixMarket.

To run the tests in your environment please use following instructions:

  1. Download and install the latest version of toolbox (green button on the top right).
  2. Download and unpack test matrices and scripts: sp_solvers_test.zip (37MB).
  3. Start MATLAB. Add toolbox into MATLAB’s search path. For example, on Windows:

    >> addpath('C:\Users\[UserName]\Documents\Multiprecision Computing Toolbox\')
  4. Run tests by command:

    >> sp_solver_test

Wait for tests to complete.

{ 2 comments… read them below or add one }

André Van Moer September 6, 2013 at 3:59 pm

Dear Pavel,

I tried the latest version (3.5.4.4586) of MP and ran the test “sp_solver_test.m”

I downloaded all the data via FTP from: ftp://math.nist.gov//pub/MatrixMarket2

Another file is necessary: mmread.m , which can be downloaded from:
http://math.nist.gov/MatrixMarket/mmio/matlab/mmread.m

I had some timings<0.00, thus I modified the format i.e.:
function print_results(fname, rows, cols, entries, timing, error)
fprintf('%15s\t\t%5d x %5d\t\tnnz = %7d\t\ttime = %8.4f sec\terror = %g\n', fname, rows, cols, entries, timing, error);
end

I will send my results in attachment via a normal email.
Best regards,

André

Reply

Pavel Holoborodko September 6, 2013 at 4:10 pm

Dear André,

Thank you for your correction!
Indeed we need mmread.m to load matrices from MatrixMarket format to MATLAB.

Actually today I have (partially) ported solvers to quadruple precision.
As a result, speed is much higher.

For example, matrix bcsstk18.mtx can now be solved in 18 seconds instead of 255!

And I had to change output format in order to see timings of the rest matrices ;).
Will release updated version shortly.

Btw, I am using old CPU: Core i7 930, timings for newer CPUs will be much better.

Reply

Leave a Comment

Previous post:

Next post: