#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cuda_runtime_api.h>
#include <cublas.h>

// includes, project
#include "magma.h"
#include "magma_lapack.h"
#include "testings.h"



void calc_magma(int n, double* A, double* w){

  double  *h_R, *d_R, *h_work;
  magma_int_t *iwork;
    
  /* Matrix size */
  magma_int_t N=0, n2;

  magma_int_t i, j, info;
  magma_int_t ione = 1, izero = 0;
  magma_int_t ISEED[4] = {0,0,0,1};

  const char *uplo = MagmaLowerStr;
  const char *jobz = MagmaVectorsStr;
  
  N = n;
  
  /* Query for workspace sizes */
  printf("Quering for workspace\n");

  double aux_work[1];
  magma_int_t aux_iwork[1];

  magma_dsyevd_gpu( jobz[0], uplo[0],
                    N, d_R, N, w,
                    h_R, N,
                    aux_work,  -1,
                    aux_iwork, -1,
                    &info );

  magma_int_t lwork, liwork;
  lwork  = (magma_int_t) aux_work[0];
  liwork = aux_iwork[0];

  /* Allocate host memory for the matrix */

  magma_int_t ldda = N;

  TESTING_HOSTALLOC( h_R, double, N*N );
  TESTING_DEVALLOC(  d_R, double, N*ldda );
  TESTING_HOSTALLOC( h_work, double,      lwork  );
  TESTING_MALLOC(    iwork,  magma_int_t, liwork );

  n2 = N*N;
  
  /* transfer the matrix */
  magma_dsetmatrix( N, N, A, N, d_R, ldda );

  /* warm up run */
  printf("Warming up...\n");
  
  magma_dsyevd_gpu( jobz[0], uplo[0],
                    N, d_R, ldda, w,
                    h_R, N,
                    h_work, lwork,
                    iwork, liwork,
                    &info );

  magma_dsetmatrix( N, N, A, N, d_R, ldda );

  /* query for optimal workspace sizes */
  printf("Quering optimal sizes...\n");

  magma_dsyevd_gpu( jobz[0], uplo[0],
                    N, d_R, ldda, w,
                    h_R, N,
                    h_work, -1,
                    iwork,  -1,
                    &info );
  
  int lwork_save  = lwork;
  int liwork_save = liwork;

  lwork  = min( lwork,  (magma_int_t) h_work[0] );
  liwork = min( liwork, iwork[0] );

  /* ====================================================================
     Performs operation using MAGMA
     =================================================================== */
  printf("Performing MAGMA operation...\n");

  magma_dsetmatrix( N, N, A, N, d_R, ldda );
  
  
  magma_dsyevd_gpu( jobz[0], uplo[0],
                    N, d_R, ldda, w,
                    h_R, N,
                    h_work, lwork,
                    iwork, liwork,
                    &info );
  if(info == 0){
    printf("Successful calculation of eigenvectors/eigenvalues\n");
  }
  
  magma_dgetmatrix( N, N, d_R, ldda, A, N );
  
}
