   v_dab_d = init_sparse_matrix_csr_gpu_float(v_dab_vals, v_dab_rowPtr, 
                   v_dab_col_ind, v_dab_shape[0], v_dab_shape[1], v_dab_nnz, v_dab_indptr_size);
 
+  // For pascal GPU, the cudaMallocManaged() could probably allow to run larger systems.
+  // Need to llok for more informations about this!!
   checkCudaErrors(cusparseCreate(&handle_cuparse));
   checkCudaErrors(cublasCreate(&handle_cublas));
 

        X4_d, norbs, &beta, nb2v_d, norbs));
  checkCudaErrors(cublasSgemm(handle_cublas, CUBLAS_OP_T, CUBLAS_OP_N, nvirt, nfermi, norbs, &alpha, &X4_d[vstart*norbs], norbs, nb2v_d,
       norbs, &beta, nm2v_re_d, nvirt));


   // imaginary part
  checkCudaErrors(cusparseScsrmv(handle_cuparse, CUSPARSE_OPERATION_NON_TRANSPOSE,

  checkCudaErrors(cublasSgemm(handle_cublas, CUBLAS_OP_T, CUBLAS_OP_N, nvirt, nfermi, norbs, &alpha, &X4_d[vstart*norbs], norbs, nb2v_d,
        norbs, &beta, nm2v_im_d, nvirt));


  // Normalization!!
  normalize_energy_gpu<<<dimGrid, dimBlock>>>(ksn2e_d, ksn2f_d, omega_re, omega_im, 
      nm2v_re_d, nm2v_im_d, nfermi, norbs, nvirt, vstart);

  // Going back: real part first
  checkCudaErrors(cublasSgemm(handle_cublas, CUBLAS_OP_N, CUBLAS_OP_N, norbs, nfermi, nvirt, &alpha, &X4_d[vstart*norbs], norbs, nm2v_re_d,
