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

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

  treal1 = QDP_create_R();
  treal2 = QDP_create_R();
  QDP_R_eq_zero(treal2, QDP_all);

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

  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(treal1, temp2, temp4, QDP_all);
      QDP_R_peq_R(treal2, treal1, QDP_all);
    }
  }
  QDP_r_eq_sum_R(&plaq, treal2, QDP_all);

  QDP_destroy_R(treal1);
  QDP_destroy_R(treal2);

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

  return plaq;
}