#include <stdio.h>
#include <qdp.h>

static QDP_ColorVector *temp0, *temp1[4], *temp2[4], *temp3[4];

/* allocate fields before calling dslash */
/* this is so we can reuse the fields & gathers */
void
prepare_dslash(void)
{
  int mu;

  temp0 = QDP_create_V();
  for (mu = 0; mu < 4; mu++) {
    temp1[mu] = QDP_create_V();
    temp2[mu] = QDP_create_V();
    temp3[mu] = QDP_create_V();
  }
}

/* free temporary fields when we're done using dslash */
void
cleanup_dslash(void)
{
  int mu;

  QDP_destroy_V(temp0);
  for (mu = 0; mu < 4; mu++) {
    QDP_destroy_V(temp1[mu]);
    QDP_destroy_V(temp2[mu]);
    QDP_destroy_V(temp3[mu]);
  }
}

/* staggered dslash */
void
dslash(QDP_ColorVector *result,
       QDP_ColorMatrix *link[],
       QDP_ColorVector *source)
{
  int mu;

  /* shift from forward neighbors into temp1 */
  for (mu = 0; mu < 4; mu++) {
    QDP_V_eq_sV(temp1[mu], source, QDP_neighbor[mu], QDP_forward, QDP_all);
  }

  /* multiply by backwards link and shift */
  for (mu = 0; mu < 4; mu++) {
    QDP_V_eq_Ma_times_V(temp2[mu], link[mu], source, QDP_all);
    QDP_V_eq_sV(temp3[mu], temp2[mu], QDP_neighbor[mu], QDP_backward, QDP_all);
  }

  /* multiply link by forward shift result and accumulate into result */
  QDP_V_eq_M_times_V(result, link[0], temp1[0], QDP_all);
  QDP_discard_V(temp1[0]);
  for (mu = 1; mu < 4; mu++) {
    QDP_V_peq_M_times_V(result, link[mu], temp1[mu], QDP_all);
    QDP_discard_V(temp1[mu]);
  }

  /* subtract backwards shift from result */
  for (mu = 0; mu < 4; mu++) {
    QDP_V_meq_V(result, temp3[mu], QDP_all);
    QDP_discard_V(temp3[mu]);
  }

}