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

QLA_Real
plaquette(QDP_ColorMatrix *link[])
{
  int mu, nu;
  QLA_Real plaq, tplaq;
  QDP_ColorMatrix *temp1, *temp2, *temp3, *temp4;

  temp1 = QDP_create_M();
  temp2 = QDP_create_M();
  temp3 = QDP_create_M();
  temp4 = QDP_create_M();

  plaq = 0;

  for(mu=0; mu<3; ++mu) {
    for(nu=mu+1; nu<4; ++nu) {

      QDP_M_eq_sM(temp1, link[nu], QDP_neighbor[mu], QDP_forward, QDP_all);
      QDP_M_eq_sM(temp2, link[mu], QDP_neighbor[nu], QDP_forward, QDP_all);

      QDP_M_eq_Ma_times_M(temp3, link[nu], link[mu], QDP_all);

      QDP_M_eq_M_times_M(temp4, temp3, temp1, QDP_all);

      QDP_r_eq_re_M_dot_M(&tplaq, temp2, temp4, QDP_all);

      plaq += tplaq;

    }
  }

  QDP_destroy_M(temp1);
  QDP_destroy_M(temp2);
  QDP_destroy_M(temp3);
  QDP_destroy_M(temp4);

  return plaq;
}