Matrix formats for lobpcg

[blattms]

Hi,

I wanted to use lobpcg to compute some eigenvalues and eigenvectors. Assuming that all matrix types should be supported I chose Magma_ELL as sparse matrix format.
Unfortunately, due to the implementation of magma_d_spmv in magma_d_blaswrapper.cpp this format is only supported for real vectors (num_cols==1). As I am searching for multiple eigenvectors my vector should be a dense matrix with several columns. Yet, Magma_SELLP is supported for my case. That seems a bit strange. Is there an algorithmic reason for this?

It would be nice to have a list of supported matrix formats somewhere.

Please also note that there is some code commented out in function magma_d_spmv that makes it a noop for some formats (e.g. MAGMA_CSR with MagmaRowMajor).

Cheers,

Markus

[Stan Tomov]

We can easily enable SpMM for any sparse matrix, e.g., just do SpMV for each vector. However, special optimizations for SpMM can do much better than SpMV. We have done this only for the Magma_SELLP format, and have disable the others - just as a reminder to add them (but we can also enable the slow implementation while we add the improvements for the others, e.g., starting from Magma_ELL).
Stan

[blattms]

Thanks for your answer Stan.
Unfortunately that would mean patching Magma, as I want to use the lobcpg solver for an ELLP matrix (and not do sparse matrix-vector product directly. They are just where it breaks..). I also assume that Magma makes use on BLAS level 3 operations (dense matriy-matrix multiplications) for multiple right hand sides and we would loose performance by this approach.
I am going for SELLP now and hope that it works out.

[blattms]

Ok, I did move to using SELLP, but it is not working yet. I am getting an error message "error: number of vectors has to be multiple of 2" in the sparse matrix vector multiplication.
No sure what this is due to.

So I tested lobpcg using testing_dsolver. To make it run through without error mesage I had to use blocksize 32 (for higher one "error: too many threads requested.
error: too many threads requested." but no error code was returned). But it seems like the results differ quite a lot when using SELLP:

Wen requesting the use of SELLP with blocksize=32 and alignment=8 I get this

Code: Select all

matrixinfo = [
%   size   (m x n)     ||   nonzeros (nnz)   ||   nnz/m   ||   stored nnz
%============================================================================%
       132       132            4347               32              5888
%============================================================================%
];
Eigenvalues:
-1.617020e+01  -1.205513e+01  -9.965225e+00  -8.563851e+00  -3.259683e+00  -2.923920e+00  -5.317256e-01  -1.116960e-03  9.093301e+00  5.388978e+01  6.490582e+01  7.867743e+01  8.791892e+01  9.520502e+01  1.079634e+02  1.094998e+02  1.163479e+02  1.246189e+02  1.305209e+02  1.355360e+02  1.502891e+02  1.649066e+02  1.672114e+02  1.879650e+02  2.056203e+02  2.406668e+02  2.580640e+02  2.768265e+02  3.229031e+02  3.427514e+02  3.691790e+02  4.366912e+02  

Final residuals:
[   48.5361  39.8537  38.7237  31.7085  23.0854  22.6478   9.4308   0.4406  19.1769  77.2629  80.1439  92.7044 115.0919 103.4186 107.4072 117.2513 106.8727 118.3167 143.3277 156.7057 113.8573 155.3311 149.2352 149.7174 173.7212 226.6655 261.0226 226.0232 319.8264 303.9063 318.7915 373.8634 ];


Residuals are stored in file residualNorms
Plot the residuals using: myplot 

convergence = [
%   iter   ||   residual-nrm2    ||   runtime    ||   SpMV-count   ||   info
%=================================================================================%
        1       0.000000e+00          0.009271                2            0
%=================================================================================%

%=================================================================================%
% LOBPCG iteration solver summary:
%    initial residual: 0.000000e+00
%    preconditioner setup: 0.0000 sec
%    iterations:    1
%    SpMV-count:    2
%    exact final residual: 0.000000e+00
%    runtime: 0.0093 sec
%    preconditioner runtime: 0.0000 sec
%=================================================================================%
];

solverinfo = [
%   iter   ||   residual-nrm2    ||   runtime    ||   SpMV-count   ||   info
%=================================================================================%
        1       0.000000e+00          0.009271                2            0
%=================================================================================%

%=================================================================================%
% LOBPCG iteration solver summary:
%    initial residual: 0.000000e+00
%    preconditioner setup: 0.0000 sec
%    iterations:    1
%    SpMV-count:    2
%    exact final residual: 0.000000e+00
%    runtime: 0.0093 sec
%    preconditioner runtime: 0.0000 sec
%=================================================================================%
];

With CSR using the same matrix

Code: Select all

matrixinfo = [
%   size   (m x n)     ||   nonzeros (nnz)   ||   nnz/m   ||   stored nnz
%============================================================================%
       132       132            4347               32              4347
%============================================================================%
];
restart at step #6
restart at step #19
restart at step #20
restart at step #21
restart at step #22
restart at step #23
restart at step #24
restart at step #25
restart at step #26
restart at step #27
restart at step #28
restart at step #29
Eigenvalues:
-9.938945e-14  -4.197015e-14  -3.383334e-14  -3.243467e-14  -2.336373e-14  -2.892851e-15  6.606877e-01  6.606877e-01  1.568701e+00  4.108731e+00  4.108731e+00  6.121250e+00  1.221696e+01  1.274937e+01  1.274937e+01  1.321201e+01  2.214689e+01  2.792283e+01  2.792283e+01  3.205128e+01  4.195567e+01  4.743933e+01  4.983593e+01  4.983593e+01  5.089055e+01  5.798131e+01  6.253386e+01  6.410256e+01  7.745958e+01  7.745958e+01  9.381312e+01  9.382367e+01  

Final residuals:
[    0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0002   0.0003 ];


Residuals are stored in file residualNorms
Plot the residuals using: myplot 
convergence = [
%   iter   ||   residual-nrm2    ||   runtime    ||   SpMV-count   ||   info
%=================================================================================%
       29       3.506094e-09          1.185907               30            0
%=================================================================================%

%=================================================================================%
% LOBPCG iteration solver summary:
%    initial residual: 3.288210e+01
%    preconditioner setup: 0.0000 sec
%    iterations:   29
%    SpMV-count:   30
%    exact final residual: 2.015266e-09
%    runtime: 1.1859 sec
%    preconditioner runtime: 0.0009 sec
%=================================================================================%
];

solverinfo = [
%   iter   ||   residual-nrm2    ||   runtime    ||   SpMV-count   ||   info
%=================================================================================%
       29       3.506094e-09          1.185907               30            0
%=================================================================================%

%=================================================================================%
% LOBPCG iteration solver summary:
%    initial residual: 3.288210e+01
%    preconditioner setup: 0.0000 sec
%    iterations:   29
%    SpMV-count:   30
%    exact final residual: 2.015266e-09
%    runtime: 1.1859 sec
%    preconditioner runtime: 0.0009 sec
%=================================================================================%
];

Seems a bit fishy to me.

[blattms]

Well, if I only request one eigenvalue, then the results seem similar. But this seems to be the only case. In all others lobpcg with SELLP thinks the initial residual is already small enough. Tested with up to 9 Eigen values.