#include <stdio.h>

#define QDP_Precision 'F'
#define QDP_Nc 3
#include <qdp.h>

int lattice_size[4] = { 4,4,4,4 };
void set_mat(QDP_ColorMatrix *mat);
void set_vec(QDP_ColorVector *vec);

int
main(int argc, char *argv[])
{
  QLA_Complex z;
  QDP_ColorMatrix *mat;
  QDP_ColorVector *v1, *v2;

  /* start QDP */
  QDP_initialize(argc, argv);

  /* set lattice size and create layout */
  QDP_set_latsize(4, lattice_size);
  QDP_create_layout();

  /* create some fields */
  mat = QDP_create_M();
  v1 = QDP_create_V();
  v2 = QDP_create_V();

  /* set values */
  set_mat(mat);
  set_vec(v1);

  /* v2[i] = mat[i] * v1[i] for all sites i */
  QDP_V_eq_M_times_V(v2, mat, v1, QDP_all);

  /* z = v1[i] . v2[i] for even sites i */
  QDP_c_eq_V_dot_V(&z, v1, v2, QDP_even);

  /* print result on node 0 */
  if(QDP_this_node==0)
    printf("result = %f%+fi\n", QLA_real(z), QLA_imag(z));

  /* free fields */
  QDP_destroy_M(mat);
  QDP_destroy_V(v1);
  QDP_destroy_V(v2);

  /* shutdown QDP */
  QDP_finalize();

  return 0;
}

void
set_vec(QDP_ColorVector *vec)
{
  void fill_vec(QLA_ColorVector *v, int coords[]);

  /* set values through a function */
  QDP_V_eq_func(vec, fill_vec, QDP_all);
}

/* function used for setting color vector */
void
fill_vec(QLA_ColorVector *v, int coords[])
{
  /* if along some wall */
  if(coords[3]==0) {
    int i;
    for(i=0; i<QDP_Nc; ++i) {
      if(i==0) {
	/* set the 0th color component to 1 */
	QLA_c_eq_r(QLA_elem_V(*v,i), 1);
      } else {
	/* otherwise set to 0 */
	QLA_c_eq_r(QLA_elem_V(*v,i), 0);
      }
    }
  } else {
    QLA_V_eq_zero(v);
  }
}

void
set_mat(QDP_ColorMatrix *mat)
{
  QLA_Int sseed;
  QDP_Int *fseed;
  QDP_RandomState *rs;
  void set_fseed(QLA_Int *seed, int coords[]);

  fseed = QDP_create_I();
  rs = QDP_create_S();

  sseed = 987654321;
  QDP_I_eq_func(fseed, set_fseed, QDP_all);
  QDP_S_eq_seed_i_I(rs, sseed, fseed, QDP_all);

  QDP_M_eq_gaussian_S(mat, rs, QDP_all);

  QDP_destroy_I(fseed);
  QDP_destroy_S(rs);
}

void
set_fseed(QLA_Int *seed, int coords[])
{
  int i,t;

  t = coords[3];
  for(i=2; i>=0; --i) {
    t = t*lattice_size[i] + coords[i];
  }
  *seed = t;
}