Sparse eigensolver example

[An Pa]

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));

    magmaDoubleComplex *val;

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

    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?

Thank you in advance!

[Stanimire Tomov]

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]

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]

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]

Hello again,

Glad everything works!

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]

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

To view this discussion on the web visit https://groups.google.com/a/icl.utk.edu/d/msgid/magma-user/a4c8270e-fd6d-4f82-862e-4dd99029ee9en%40icl.utk.edu.

[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]

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:

[An Pa]

Thank you, I've found where I am mistaken. 

пятница, 11 декабря 2020 г. в 13:27:38 UTC+3, An Pa: