#include "corrX.h"

using namespace std;
using namespace scl;

int main(int argc, char** argp, char** envp)
{
  cudaError_t err;

  int deviceCount;
  cudaGetDeviceCount(&deviceCount);

  int device;
  cudaDeviceProp deviceProp;
  bool device_has_double = false;

  for (device = 0; device < deviceCount; ++device) {
    cudaGetDeviceProperties(&deviceProp, device);
    if (deviceProp.major == 1 && deviceProp.minor == 3) {
      device_has_double = true;
      err = cudaSetDevice(device);
      if (err) error("Error, cudaSetDevice, err ="+fmt('d',3,err).get_ostr());
      break;
    }
  }

  if (!device_has_double) {
    for (device = 0; device < deviceCount; ++device) {
      cudaGetDeviceProperties(&deviceProp, device);
      if (deviceProp.major != 9999 && deviceProp.minor != 9999) {
        err = cudaSetDevice(device);
        if (err) error("Error, cudaSetDevice, err ="+fmt('d',3,err).get_ostr());
        break;
      }
      else {
        error("Error, no non_emulation cuda device available");
      }
    }
  }

  cout << '\n';
  cout << "Device set is " << deviceProp.name << ' ' << "\n";
  cout << "Device number is " << device << '\n';
  if (device_has_double) cout << "Device has double" << '\n';

  realmat X;    // everything is presumed float, not double

  vecread("yields.dat",X);

  size_t sizeX = X.size();
  if (sizeX>deviceProp.maxThreadsDim[0]*deviceProp.maxGridSize[0]) {
    error("Error, corrX, X is too large");
  }

  size_t bytes = sizeX*sizeof(X[0]);

  float* devX;
  err = cudaMalloc((void**)&devX, bytes);
  if (err != cudaSuccess) error("Error, cudaMalloc failed");

  cudaMemcpyKind kind = cudaMemcpyHostToDevice;

  err = cudaMemcpy(devX, X.begin(), bytes, kind);
  if (err != cudaSuccess) error("Error, cudaMemcpy failed");

  realmat C(X.ncol(),X.ncol(),0.0);

  size_t sizeC = C.size();
  bytes = sizeC*sizeof(X[0]);

  float* devC;

  err = cudaMalloc((void**)&devC, bytes);
  if (err != cudaSuccess) error("Error, cudaMalloc failed");
 
  bool success = corrXdev(devX, (size_t)X.nrow(), (size_t)X.ncol(), &devC);
  if (!success) error("Error, corrXdev failed");

  bytes = sizeC*sizeof(X[0]);

  kind = cudaMemcpyDeviceToHost;

  err = cudaMemcpy(C.begin(), devC, bytes, kind);
  if (err != cudaSuccess) error("Error, cudaMemcpy failed");

  cout << C << '\n';

  cudaFree(devX);
  cudaFree(devC);

  cudaThreadExit();

  return 0;
}