/**

********************************************************************************
@page MAGMA-sparse Sparse Overview


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 {#datastructures}
=================================
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 {#io}
=================================
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 {#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 {#memorhandling}
=================================
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 {#sparsesolvers}
=================================
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 {#sparseexample}
=================================
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();

*/