#include <stdio.h>
#include <string.h>

#include "qmp.h"
#include "sse-dwf-cg.h"

double
DWF_gauge_reader(const void *ptr, void *env, 
		 int latt_coord[4], int mu, int row, int col, int reim)
{
  double val;

  if (row == col)
    val = (reim == 0) ? 1.0 : 0.0;
  else
    val = 0.0;

  return val;
}

double
DWF_fermion_reader1(const void *ptr, void *env, 
		    int latt_coord[5], int color, int spin, int reim)
{
  return 1.0;
}

double
DWF_fermion_reader2(const void *ptr, void *env, 
		    int latt_coord[5], int color, int spin, int reim)
{
  return 0.0;
}

void
DWF_fermion_writer(void *ptr, void *env, 
		   int latt_coord[5], int color, int spin, int reim,
		   double val)
{
  return;
}


#define ND 4
#define LX  4
#define LY  4
#define LZ  4
#define LT  4
#define LS  8
#define VOL  (LX*LY*LZ*LT)

#define SU3  72
#define LF   (24*LS)

typedef float  OuterGauge[VOL*SU3];
typedef float  OuterFermion[VOL*LF];


int main(int argc, char **argv)
{
  int length[ND+1] = {LX, LY, LZ, LT, LS};
  int i,err;
  int num;
  SSE_DWF_Gauge*   g;
  SSE_DWF_Fermion* rhs;
  SSE_DWF_Fermion* x0;
  SSE_DWF_Fermion* x;
  OuterGauge u, v;
  OuterFermion psi, chi;

  QMP_init_msg_passing(&argc, &argv, QMP_SMP_ONE_ADDRESS);

  /* Let the geometry arrange itself */
  fprintf(stderr,"Dynamical machine size\n");
  err = QMP_layout_grid(length, ND);
  num = QMP_get_number_of_nodes();
  {
    const int *s = QMP_get_subgrid_dimensions();
    for(i=0; i < ND; ++i)
      fprintf(stderr,"subgrid[%d] = %d\n",i,s[i]);
  }

  err = SSE_DWF_init(length, SSE_DWF_FLOAT, NULL, NULL);
  if (err)
    QMP_error_exit("error initializing SSE_DWF_init");

  // Initialize the SSE specific fields
  if ((g = SSE_DWF_load_gauge(u, v, NULL, DWF_gauge_reader)) == NULL)
    QMP_error_exit("error initializing g");

  if ((rhs = SSE_DWF_load_fermion(chi, NULL, DWF_fermion_reader1)) == NULL)
    QMP_error_exit("error initializing rhs");

  if ((x0 = SSE_DWF_load_fermion(psi, NULL, DWF_fermion_reader2)) == NULL)
    QMP_error_exit("error initializing x0");

  if ((x = SSE_DWF_allocate_fermion()) == NULL)
    QMP_error_exit("error initializing x");

  {
    // Call the solver
    double out_epsilon = 0.0;
    double M_0 = 1.0;    // Warning check sign convention
    double m_q = 0.3;
    double Rsd = 1.0e-5;
    int iter = 0;
    int MaxCG = 1000;

    printf("M_0=%g  m_q=%g  Rsd=%g  out_eps=%g\n",
	   M_0, m_q, Rsd, out_epsilon);

    err = SSE_DWF_cg_solver(x, &out_epsilon, &iter,
			    g, M_0, m_q,
			    x0, rhs, 
			    Rsd, MaxCG);
    if (err)
      QMP_error_exit("error in cg solver");
  }

  /* Save the result */
  SSE_DWF_save_fermion(psi, NULL, DWF_fermion_writer, x);

  /* Cleanup */
  SSE_DWF_delete_fermion(x);
  SSE_DWF_delete_fermion(x0);
  SSE_DWF_delete_fermion(rhs);
  SSE_DWF_delete_gauge(g);

  /* Bolt */
  SSE_DWF_fini();
  QMP_finalize_msg_passing();

  return 0;
}