#include <stdlib.h>
   #include <stdio.h>
   #include <string.h>
   #include <math.h>
   #include "magma_v2.h"
   #include "magmasparse.h"

  double c_gmres_gpu(
     int m,
     int n,
     magma_index_t *rowPointerA,
     magma_index_t *colIndA,
     double *valA,
     double *valB, double **psol )
     {


    *psol = (double*) calloc(m, sizeof(double));

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

     magma_d_matrix A={Magma_CSR}, dA={Magma_CSR};
     magma_d_vector B={Magma_CSR}, dB={Magma_CSR}, dx={Magma_CSR},x={Magma_CSR};



//   pass linear system to MAGMA-sparse
     magma_dcsrset( m, m, rowPointerA, colIndA, valA,&A, queue );

//   charger les vecteurs RHS
     magma_dvset( m, n, valB, &B, queue );
     magma_dvset( m, n, *psol, &x, queue );

//  configure solver
    dopts.solver_par.solver     = Magma_GMRES;
    dopts.solver_par.restart    = 50;
    dopts.solver_par.maxiter    = 1000;
    dopts.solver_par.rtol       = 1e-10;
    dopts.solver_par.atol       = 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 );
    magma_dsolverinfo( &dopts.solver_par, &dopts.precond_par, queue );


//   copier la matrice A sur le GPU
     magma_dmtransfer( A, &dA, Magma_CPU, Magma_DEV, queue );
     magma_dmtransfer( B, &dB, Magma_CPU, Magma_DEV, queue );
   

//   allocate solution vector - on device
//   magma_d_vinit( &dx, Magma_DEV, m, 1.0, queue );
     magma_dmtransfer( x, &dx, Magma_CPU, Magma_DEV, queue );


//  solve the linear system
//    magma_d_solver (dA, dB, &dx, &dopts, queue);
      magma_dfgmres(dA,dB,&dx,&dopts.solver_par, &dopts.precond_par,queue);
 
//  copy the solution vector back to the host and pass it back to the application
   
    magma_dmtransfer( dx, &x, Magma_DEV, Magma_CPU, queue );
    magma_dvget( x, &m, &n, psol, queue );

//  clean up the memory

    magma_dsolverinfo_free( &dopts.solver_par, &dopts.precond_par, queue );
    magma_d_mfree( &dA, queue );
    magma_d_mfree( &dx, queue );
    magma_d_mfree( &dB, queue );
    magma_d_mfree( &A, queue );
    magma_d_mfree( &B, queue );
  
//  and finalize MAGMA.
    magma_queue_destroy( queue );
    magma_finalize();


    //printf("valAB\n");
    //for (i = 0; i < 10; ++i) {
    //    printf("%.10f\n", (*psol)[i]);
    //};

return 0;
}