The MAGMA-sparse Package

The MAGMA-sparse package in the MAGMA software stack contains sparse BLAS routines as well as functions to handle the complete iterative solution process of a sparse linear system of equations.

A typical application includes:

an interface passing the linear system to MAGMA
choosing the desired data structures for the respective sparse BLAS functions
sending the data structures to the device
choosing solver, eigensolver, and preconditioner
solving the respective system on the device
passing back the results

For each of these steps, multiple options are offered by the MAGMA software stack.

To start, in C or C++, include the magma_v2.h and magmasparse.h headers.

#include <magma_v2.h>
#include <magmasparse.h>

Sparse Data Structures

For a more generic programming approach, the sparse data structures (matrices and vectors) are stored in data structures containing all information necessary to access the sparse-BLAS via wrappers:

struct magma_z_matrix
{
    magma_storage_t    storage_type;            // matrix format - CSR, ELL, SELL-P
    magma_location_t   memory_location;         // CPU or DEV
    magma_symmetry_t   sym;                     // opt: indicate symmetry
    magma_diagorder_t  diagorder_type;          // opt: only needed for factorization matrices
    magma_uplo_t   fill_mode;               // fill mode full/lower/upper
    magma_int_t        num_rows;                // number of rows
    magma_int_t        num_cols;                // number of columns
    magma_int_t        nnz;                     // opt: number of nonzeros
    magma_int_t        max_nnz_row;             // opt: max number of nonzeros in one row
    magma_int_t        diameter;                // opt: max distance of entry from main diagonal
    union {
        magmaDoubleComplex      *val;           // array containing values in CPU case
        magmaDoubleComplex_ptr  dval;           // array containing values in DEV case
    };
    union {
        magmaDoubleComplex      *diag;          // opt: diagonal entries in CPU case
        magmaDoubleComplex_ptr  ddiag;          // opt: diagonal entries in DEV case
    };
    union {
        magma_index_t           *row;           // opt: row pointer CPU case
        magmaIndex_ptr          drow;           // opt: row pointer DEV case
    };
    union {
        magma_index_t           *rowidx;        // opt: array containing row indices CPU case
        magmaIndex_ptr          drowidx;        // opt: array containing row indices DEV case
    };
    union {
        magma_index_t           *col;           // opt: array containing col indices CPU case
        magmaIndex_ptr          dcol;           // opt: array containing col indices DEV case
    };
    union {
        magma_index_t           *list;          // opt: linked list pointing to next element
        magmaIndex_ptr          dlist;          // opt: linked list pointing to next element
    };
    magma_index_t      *blockinfo;              // opt: for BCSR format CPU case
    magma_int_t        blocksize;               // opt: info for SELL-P/BCSR
    magma_int_t        numblocks;               // opt: info for SELL-P/BCSR
    magma_int_t        alignment;               // opt: info for SELL-P/BCSR
    magma_order_t      major;                   // opt: row/col major for dense matrices
    magma_int_t        ld;                      // opt: leading dimension for dense
} magma_z_matrix;

The purpose of the unions (e.g. for the val array) is to support different hardware platforms, where the magmaXXX_ptr is adapted to the respective device characteristics.

For sparse matrices, the main formats are CSR, ELL and the MAGMA-specific SELL-P. Generally, the sparse-BLAS routines provide for these the best performance.

Without specifying the storage type, the memory location or the dimension of the matrix, the sparse matrix vector product can then be used via the wrapper:

magma_z_spmv(
        magmaDoubleComplex alpha,
        magma_z_matrix A,
        magma_z_matrix x,
        magmaDoubleComplex beta,
        magma_z_matrix y,
        magma_queue_t queue );

Sparse I/O

A sparse matrix stored in mtx format can be read from disk via the function:

magma_z_csr_mtx(
        magma_z_matrix *A,
        const char *filename,
        magma_queue_t queue );

A sparse linear system of dimension mxn present in main memory can be passed to MAGMA via the function:

magma_zcsrset(
        magma_int_t m,
        magma_int_t n,
        magma_index_t *row,
        magma_index_t *col,
        magmaDoubleComplex *val,
        magma_z_matrix *A,
        magma_queue_t queue );

where row, col, val contain the matrix in CSR format.

If the matrix is already present in DEV memory, the corresponding function is

magma_zcsrset_gpu(
    magma_int_t m,
    magma_int_t n,
    magmaIndex_ptr row,
    magmaIndex_ptr col,
    magmaDoubleComplex_ptr val,
    magma_z_matrix *A,
    magma_queue_t queue );

Similarly, matrices handled in MAGMA can be returned via the functions

magma_zcsrget(
    magma_z_matrix A,
    magma_int_t *m,
    magma_int_t *n,
    magma_index_t **row,
    magma_index_t **col,
    magmaDoubleComplex **val,
    magma_queue_t queue );

magma_zcsrget_gpu(
    magma_z_matrix A,
    magma_int_t *m,
    magma_int_t *n,
    magmaIndex_ptr *row,
    magmaIndex_ptr *col,
    magmaDoubleComplex_ptr *val,
    magma_queue_t queue );

respectively

write_z_csrtomtx(
    magma_z_matrix A,
    const char *filename,
    magma_queue_t queue );

Additionally, MAGMA contains routines to generate stencil discretizations of different kind.

Vectors are handled as dense matrices (Magma_DENSE) and can be initialized inside MAGMA via

magma_z_vinit(
    magma_z_matrix *x,
    magma_location_t memory_location,
    magma_int_t num_rows,
    magma_int_t num_cols,
    magmaDoubleComplex values,
    magma_queue_t queue );

where ''memory_location'' sets the location (Magma_CPU or Magma_DEV). Also, vectors can be read from file via

magma_z_vread(
    magma_z_matrix *x,
    magma_int_t length,
    char * filename,
    magma_queue_t queue );

or - in case of a block of sparse vectors stored as CSR matrix - via

magma_z_vspread(
    magma_z_matrix *x,
    const char * filename,
    magma_queue_t queue );

or passed from/to main memory:

magma_zvset(
    magma_int_t m,
    magma_int_t n,
    magmaDoubleComplex *val,
    magma_z_matrix *b,
    magma_queue_t queue );

magma_zvget(
    magma_z_matrix b
    magma_int_t *m,
    magma_int_t *n,
    magmaDoubleComplex **val,
    magma_queue_t queue );

Matrix Formats

To convert a matrix from one into another format, the CPU-based routine

magma_z_mconvert(
    magma_z_matrix A,
    magma_z_matrix *B,
    magma_storage_t old_format,
    magma_storage_t new_format,
    magma_queue_t queue );

can be used where old_format and new_format determine the specific conversion.

Memory Handling

All iterative solvers and eigensolvers included in the MAGMA-sparse package work on the device. Hence, it is required to send the respective data structures to the device for solving, and back to access the solution. The functions

magma_z_mtransfer(
    magma_z_matrix A,
    magma_z_matrix *B,
    magma_location_t src,
    magma_location_t dst,
    magma_queue_t queue );

magma_z_vtransfer(
    magma_z_matrix x,
    magma_z_matrix *y,
    magma_location_t src,
    magma_location_t dst,
    magma_queue_t queue );

allow any data copy operation - from host to device, device to device, device to host or host to host.

Linear algebra objects can be deallocated via

magma_z_mfree(
    magma_z_matrix *A,
    magma_queue_t queue );

Iterative Solvers

The MAGMA-sparse package contains a variety of linear solvers, eigensolvers, and preconditioners. The standard procedure to call a solver is to pass the linear algebra objects (located on the device) and a structure called magma_z_solver_par (respectively and magma_z_precond_par) controlling the iterative solver and collecting information during the execution:

struct magma_z_solver_par
{
    magma_solver_type  solver;                  // solver type
    magma_int_t        version;                 // sometimes there are different versions
    double             atol;                     // absolute residual stopping criterion
    double             rtol;                    // relative residual stopping criterion
    magma_int_t        maxiter;                 // upper iteration limit
    magma_int_t        restart;                 // for GMRES
    magma_ortho_t      ortho;                   // for GMRES
    magma_int_t        numiter;                 // feedback: number of needed iterations
    double             init_res;                // feedback: initial residual
    double             final_res;               // feedback: final residual
    double             iter_res;                // feedback: iteratively computed residual
    real_Double_t      runtime;                 // feedback: runtime needed
    real_Double_t      *res_vec;                // feedback: array containing residuals
    real_Double_t      *timing;                 // feedback: detailed timing
    magma_int_t        verbose;                 // print residual ever 'verbose' iterations
    magma_int_t        num_eigenvalues;         // number of EV for eigensolvers
    magma_int_t        ev_length;               // needed for framework
    double             *eigenvalues;            // feedback: array containing eigenvalues
    magmaDoubleComplex_ptr      eigenvectors;   // feedback: array containing eigenvectors on DEV
    magma_int_t        info;                    // feedback: did the solver converge etc.

the input for verbose is: 0 = production mode k>0 = convergence and timing is monitored in *res_vec and *timeing every k-th iteration

the output of info is: 0 = convergence (stopping criterion met) -1 = no convergence

-2 = convergence but stopping criterion not met within maxiter

} magma_z_solver_par;

These entities can either be initialized manually, or via the function

magma_zsolverinfo_init(
    magma_z_solver_par *solver_par,
    magma_z_preconditioner *precond_par,
    magma_queue_t queue );

setting them to some default values. For eigensolvers, the workspace needed for the eigenvectors has to be allocated consistent with the matrix dimension, which requires additionally calling

magma_zeigensolverinfo_init(
    magma_z_solver_par *solver_par,
    magma_queue_t queue );

after setting solver_par.ev_length and solver_par.num_eigenvalues to the correct numbers.

For the preconditioner configuration, a similar structure is used:

typedef struct magma_z_preconditioner { magma_solver_type solver; magma_solver_type trisolver; magma_int_t levels; magma_int_t sweeps; magma_int_t pattern; magma_int_t bsize; magma_int_t offset; magma_precision format; double atol; double rtol; magma_int_t maxiter; magma_int_t restart; magma_int_t numiter; magma_int_t spmv_count; double init_res; double final_res; real_Double_t runtime; // feedback: preconditioner runtime real_Double_t setuptime; // feedback: preconditioner setup time magma_z_matrix M; magma_z_matrix L; magma_z_matrix LT; magma_z_matrix U; magma_z_matrix UT; magma_z_matrix LD; magma_z_matrix UD; magma_z_matrix d; magma_z_matrix d2; magma_z_matrix work1; magma_z_matrix work2; magma_int_t* int_array_1; magma_int_t* int_array_2; cusparseSolveAnalysisInfo_t cuinfo; // for cuSPARSE ILU cusparseSolveAnalysisInfo_t cuinfoL; // for cuSPARSE ILU cusparseSolveAnalysisInfo_t cuinfoLT; // for cuSPARSE ILU cusparseSolveAnalysisInfo_t cuinfoU; // for cuSPARSE ILU cusparseSolveAnalysisInfo_t cuinfoUT; // for cuSPARSE ILU } magma_z_preconditioner;

An easy way to access the data collected during a solver execution is given by the function

magma_zsolverinfo(
    magma_z_solver_par *solver_par,
    magma_z_preconditioner *precond_par,
    magma_queue_t queue );

After completion,

magma_zsolverinfo_free(
    magma_z_solver_par *solver_par,
    magma_z_preconditioner *precond,
    magma_queue_t queue );

deallocates all memory used within the solver and preconditioner structure.

A solver can be called via the wrapper

magma_z_solver(
    magma_z_matrix A, magma_z_matrix b,
    magma_z_matrix *x, magma_zopts *zopts,
    magma_queue_t queue );

where zopts is a structure containing both, the solver and the preconditioner information: struct magma_zopts{ magma_z_solver_par solver_par; magma_z_preconditioner precond_par; magma_storage_t input_format; int blocksize; int alignment; magma_storage_t output_format; magma_location_t input_location; magma_location_t output_location; magma_scale_t scaling; }magma_zopts;

All entities of this structure can be initialized from command line by calling

magma_zparse_opts(
    int argc,
    char** argv,
    magma_zopts *opts,
    int *matrices,
    magma_queue_t queue );

Example

Especially when using MAGMA-sparse for the first time, the easiest way to get familiar with the package is to use and modify one of the predefined testers.

In the following example we assume to have an application coded in C and running in double precision that at some instance requires solving a linear system of the form Ax=b, where A and b are generated within the application. Furthermore, we assume that the matrix A is in CSR format and stored in row/col/val, and the vector b is stored in valb. These entities are passed to MAGMA-sparse. In the end, we want the solution passed back to the application in the vector valx.

initialize MAGMA-sparse and create some LA objects magma_dopts dopts; magma_queue_t queue; magma_queue_create( 0, &queue ); magma_d_matrix A, A_d, x, x_d, b, b_d;

pass linear system to MAGMA-sparse magma_dcsrset( m, n, row, col, val, &A, queue ); magma_dvset( m, n, valb, &b, queue );

copy the linear system to the device magma_d_vtransfer( b, &b_d, Magma_CPU, Magma_DEV, queue ); magma_d_mtransfer( A, &A_d, Magma_CPU, Magma_DEV, queue ); allocate solution vector - on device magma_d_vinit( &x_d, Magma_DEV, A.num_cols, 1, one, queue );

configure solver dopts.solver_par.solver = Magma_PGMRES; dopts.solver_par.restart = 30; dopts.solver_par.rtol = 1e-10; magma_dsolverinfo_init( &dopts.solver_par, &dopts.precond_par, queue );

configure the preconditioner dopts.precond_par.solver = Magma_ILU; dopts.precond_par.levels = 0; dopts.precond_par.trisolver = Magma_CUSOLVE; magma_d_precondsetup( A, b, &dopts.solver_par, &dopts.precond_par, queue );

solve the linear system magma_d_solver( A_d, b_d, &x_d, &dopts, queue );

copy the solution vector back to the host and pass it back to the application magma_d_mtransfer( x_d, &x, Magma_DEV, Magma_CPU, queue ); magma_dvget( x, &m, &n, &valx, queue );

clean up the memory magma_dsolverinfo_free( &dopts.solver_par, &dopts.precond_par, queue ); magma_d_mfree(&A_d, queue ); magma_d_mfree(&A, queue ); magma_d_mfree(&x_d, queue ); magma_d_mfree(&b_d, queue ); magma_d_mfree(&b, queue );

finalize MAGMA magma_queue_destroy( queue ); magma_finalize();