# Sparse eigensolver example

59 views

### An Pa

Dec 9, 2020, 10:43:06 AM12/9/20
to MAGMA User
Hello everyone,

I need to calculate the smallest eigenvalues of a huge  sparse matrix. I did not find any examples.  Could you guide me a bit?

I used the standard example and also some notes from Magma User Forum.

My "code":
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include  <cuda.h>

#include "magma_v2.h"
#include "magmasparse.h"
#include "magma_lapack.h"
#include <iostream>

using namespace std;

int main()
{
int n = 3;
int *col;
int *row;

row = (int*) calloc(n+1, sizeof(int));
col = (int*) calloc(n,   sizeof(int));

for (int i = 0; i < n; ++i) {
col[i] = i;
row[i] = i;
val[i] = MAGMA_Z_MAKE(1.0,0.0);
}

if (val == NULL) {
fprintf( stderr, "malloc failed\n" );
return 0;
}
else {
fprintf( stderr, "malloc succeeded\n" );
}

//Initialize MAGMA
magma_init();
magma_zopts opts;
magma_queue_t queue;
magma_queue_create(0, &queue);

memset(&opts, 0, sizeof(magma_zopts));

magma_z_matrix A={Magma_CSR}, dA={Magma_CSR};
magma_zcsrset(n,n,row,col,val,&A,queue);

opts.solver_par.solver = Magma_LOBPCG; // choose an LOBPCG solver
opts.solver_par.num_eigenvalues = 3; // number of eigenvalues you want to compute
opts.solver_par.maxiter = 100; // max number of iterations
opts.solver_par.rtol = 1e-8; // stopping criterion - relative accuracy of first eigenvalue
opts.precond_par.solver = Magma_ILU; // preconditioner
opts.precond_par.levels = 0; // ILU(0) - no fill-in
opts.precond_par.trisolver = Magma_CUSOLVE; //exact triangular solves

magma_zsolverinfo_init(&opts.solver_par,&opts.precond_par, queue );
magma_zeigensolverinfo_init(&opts.solver_par, queue );

magma_zmtransfer(A, &dA, Magma_CPU, Magma_DEV, queue);

magma_zlobpcg(dA, &opts.solver_par, &opts.precond_par, queue);

std::cout<<opts.solver_par.eigenvalues<<std::endl;

magma_zmfree( &dA, queue );
magma_zmfree( &A, queue );

magma_queue_destroy( queue );
magma_finalize();

return 0;
}

Output:
Eigenvalues:
0.000000e+00  0.000000e+00  0.000000e+00
What is wrong?

### Stanimire Tomov

Dec 10, 2020, 3:50:25 AM12/10/20
to An Pa, MAGMA User
Hi,

The examples are in the ‘testing’ directory. If you go there, you can do for example
make testing_dsolver
./testing_dsolver   --solver LOBPCG --ev 3 LAPLACE2D 100

This will generate 5 point stencil sparse matrix on domain 100x100, so matrix
size will be 10,000 x 10,000.

You can also test on matrices from the Florida matrix collection, e.g.,
./testing_dsolver --solver LOBPCG --ev 8 test_matrices/Trefethen_2000.mtx

Related to your test, I think you forgot to initialize where the last row ends, i.e., add

row[n]=n;

Also, I had to set this
opts.solver_par.ev_length = n;

For some reason, we don’t deduce that the length of the vectors is n and have to set it.
After this I get

Eigenvalues:
1.000000e+00  1.000000e+00  1.000000e+00

Hope this helps.

Best regards,
Stan

--
You received this message because you are subscribed to the Google Groups "MAGMA User" group.
To unsubscribe from this group and stop receiving emails from it, send an email to magma-user+...@icl.utk.edu.
To view this discussion on the web visit https://groups.google.com/a/icl.utk.edu/d/msgid/magma-user/da3c2405-1580-443c-86e5-3fdac26aa07fn%40icl.utk.edu.

### An Pa

Dec 10, 2020, 8:29:54 AM12/10/20
to MAGMA User, to...@icl.utk.edu, MAGMA User, An Pa
Dear Stan,

Thank you very much for your help. It helped me!

If I may, I ask one more question. In which way you could suggest me to modify this code (say, to improve performance) or everything more or less is Ok? Of course, I understand, that ideas regarding the modifications depend on the matrix form, but may be it is possible to say something in general?

четверг, 10 декабря 2020 г. в 03:50:25 UTC-5, to...@icl.utk.edu:

### An Pa

Dec 10, 2020, 9:24:56 AM12/10/20
to MAGMA User, An Pa, to...@icl.utk.edu, MAGMA User
Dear Stan,

I've just faced with another problem:

If I start to fill matrix by different values I obtain: Floating point exception (core dumped). What is wrong?

Thank you again!

Filling part:

for (int i = 0; i < 900; ++i) {

col[i] = i;
row[i] = i;
val[i] = MAGMA_Z_MAKE(1.0,0.0);
}

for (int i = 900; i < 1800; ++i) {

col[i] = i;
row[i] = i;
val[i] = MAGMA_Z_MAKE(-3.0,0.0);
}

for (int i = 1800; i < 2700; ++i) {

col[i] = i;
row[i] = i;
val[i] = MAGMA_Z_MAKE(2.0,0.0);
}

row[n] = n;
четверг, 10 декабря 2020 г. в 08:29:54 UTC-5, An Pa:

### Stanimire Tomov

Dec 10, 2020, 12:05:42 PM12/10/20
to An Pa, MAGMA User
Hello again,

Are you referring to modifying testing_dsolver or the LOBPCG code?
Also, if it is for LOBPCG, probably algorithmically some parts can be modified based on the problems.
Performance of kernels probably also can be improved but you would have to profile them and compare with
some theoretical expectations (like roofline models).

Stan

### Stanimire Tomov

Dec 10, 2020, 12:11:27 PM12/10/20
to An Pa, MAGMA User
I didn’t run this case. Is n = 2700 for this example.
One thing that I notice is that you put some negative numbers and some positive,
so the matrix is not positive definite (or negative definite). LOBPCG requires definiteness.
We have to look at the code though to see is we can detect this at runtime and exit nicely
and notify, instead of crashing like this, if this is the problem.

We are open to contribution - if you are interested to investigate more cases like this and
improve the code, as you asked about performance improvement, that will be great.

Stan

### An Pa

Dec 11, 2020, 2:51:04 AM12/11/20
to MAGMA User, to...@icl.utk.edu, MAGMA User, An Pa
Dear Stan,

I fill the array, for example, as diag{1,....1,3,....3,1,...,1}(positive definiteness), but dump appeared again.

I am sure that my "code" is wrong :) Probably it is better idea(in order to calculate my needs) to take as a basis the standard example ( testing_zsolver). It is appropriate within the present topic to discuss the structure of the "testing_zsolver.cpp" in order to remake it for my purposes, if some questions appear?

четверг, 10 декабря 2020 г. в 12:11:27 UTC-5, to...@icl.utk.edu:

### An Pa

Dec 11, 2020, 5:27:38 AM12/11/20
to MAGMA User, An Pa, to...@icl.utk.edu, MAGMA User
Dear Stan,

I've taken the source from the testing example:

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "magma_v2.h"
#include "magmasparse.h"

#include <iostream>

using namespace std;

int main(  int argc, char** argv )
{
int n = 10000;
int *col;
int *row;

double gamma0, omega, eta;
gamma0 = 0.5;
omega=1.0;
eta=1.0;

row = (int*) calloc(n+1, sizeof(int));
col = (int*) calloc(n,   sizeof(int));

double *val;

val = (double*) calloc(n, sizeof(double));

for (int i = 0; i < 5000; ++i) {

col[i] = i;
row[i] = i;
val[i] = MAGMA_D_MAKE(2.0,0.0);
}

for (int i = 5000; i < 10000; ++i) {

col[i] = i;
row[i] = i;
val[i] = MAGMA_D_MAKE(1.0,0.0);
}

row[n] = n;

if (val == NULL) {
fprintf( stderr, "malloc failed\n" );
return 0;
}
else {
fprintf( stderr, "malloc succeeded\n" );
}

magma_int_t info = 0;
magma_init();
magma_print_environment();

magma_dopts zopts;
magma_queue_t queue;
magma_queue_create( 0, &queue );

// double zero = MAGMA_D_MAKE(0.0, 0.0);
magma_d_matrix A={Magma_CSR}, B={Magma_CSR}, dB={Magma_CSR};
magma_d_matrix x={Magma_CSR}, b={Magma_CSR};

int i=1;
magma_dparse_opts( argc, argv, &zopts, &i, queue );
B.blocksize = zopts.blocksize;
B.alignment = zopts.alignment;

magma_dsolverinfo_init( &zopts.solver_par, &zopts.precond_par, queue );

magma_dcsrset(n,n,row,col,val,&A,queue);

// for the eigensolver case
zopts.solver_par.ev_length = A.num_cols;
magma_deigensolverinfo_init( &zopts.solver_par, queue );

// scale matrix
magma_dmscale( &A, zopts.scaling, queue );

// preconditioner
if ( zopts.solver_par.solver != Magma_ITERREF ) {
magma_d_precondsetup( A, b, &zopts.solver_par, &zopts.precond_par, queue  );
}

magma_dmconvert( A, &B, Magma_CSR, zopts.output_format, queue );

printf( "\n%% matrix info: %lld-by-%lld with %lld nonzeros\n\n",
(long long) A.num_rows, (long long) A.num_cols, (long long) A.nnz );

printf("matrixinfo = [\n");
printf("%%   size   (m x n)     ||   nonzeros (nnz)   ||   nnz/m   ||   stored nnz\n");
printf("%%============================================================================%%\n");
printf("  %8lld  %8lld      %10lld             %4lld        %10lld\n",
(long long) B.num_rows, (long long) B.num_cols, (long long) B.true_nnz,
(long long) (B.true_nnz/B.num_rows), (long long) B.nnz );
printf("%%============================================================================%%\n");
printf("];\n");

magma_dmtransfer( B, &dB, Magma_CPU, Magma_DEV, queue );

magma_dvinit_rand( &b, Magma_DEV, A.num_rows, 1, queue );

magma_dvinit_rand( &x, Magma_DEV, A.num_cols, 1, queue );

info = magma_d_solver( dB, b, &x, &zopts, queue );
if( info != 0 ) {
printf("%%error: solver returned: %s (%lld).\n",
magma_strerror( info ), (long long) info );
}

zopts.solver_par.verbose = 0;

magma_dmfree(&dB, queue );
magma_dmfree(&B, queue );
magma_dmfree(&A, queue );
magma_dmfree(&x, queue );
magma_dmfree(&b, queue );

magma_queue_destroy( queue );
magma_finalize();
return info;
}

Everything works well!

However, if I start to add non-diagonal elements as:

int n = 10000;
int *col;
int *row;

double gamma0, omega, eta;
gamma0 = 0.5;
omega=1.0;
eta=1.0;

row = (int*) calloc(n+1+1, sizeof(int));
col = (int*) calloc(n+1,   sizeof(int));

double *val;

val = (double*) calloc(n+1, sizeof(double));

for (int i = 0; i < 5000; ++i) {

col[i] = i;
row[i] = i;
val[i] = MAGMA_D_MAKE(2.0,0.0);
}

for (int i = 5000; i < 10000; ++i) {

col[i] = i;
row[i] = i;
val[i] = MAGMA_D_MAKE(1.0,0.0);
}

col[10000] = 100;
row[10000] = 10;
val[10000] = MAGMA_D_MAKE(2.0,0.0);

// It is assumed that symmetric part is automatically added? (I did with it and it dod not help me)
//    col[10001] = 10;
//    row[10001] = 100;
//    val[10001] = MAGMA_D_MAKE(2.0,0.0);

row[n+1] = n+1;

In this case I again obtain memory dump, could you correct me? Thank you in advance!

пятница, 11 декабря 2020 г. в 02:51:04 UTC-5, An Pa:
Message has been deleted