/********************* updatehiggs.c ************************/

#include "lattice.h"

void get_adjstaple(int parity, adjoint_matrix *as);

void updatehiggs(int relax)
{
  int parity;
  double rtot;
  adjoint_matrix *adjstaple;

  adjstaple = tmp_latfield( adjoint_matrix );

  forbothparities(parity) {

    if (is_multicanonical) set_mc_update(parity);

    get_adjstaple(parity, adjstaple);
    
    /* Here Goverrelax or Xoverrelax */
    if (relax) rtot = Goverrelax(parity, ahiggs, adjstaple);
    else rtot = HBHiggs(parity, ahiggs, adjstaple);

    if (is_multicanonical) mc_acceptance(parity, rtot);
  }

  free_tmp( adjstaple );
}


/*************************************************************
 *                                                           *
 *     calculate the gauge + adjoint Higgs link              *
 *                                                           *
 ************************************************************/

void get_adjstaple(int parity, adjoint_matrix *adjstaple)
{
  /*     calculate the adjoint-gauge link action ('staple')
   *
   *     ab |-- ub --> X |-- uf --> af
   *
   *     2 beta tr (a [u an u'])
   *
   *     Note that [ ] is also 'adjoint!'
   */

  register int i,dir, odir, otherparity;
  msg_tag *tag0,*tag1;
  su3_matrix tmat1,tmat2;
  adjoint_matrix tadj;

  /* Loop over directions, computing force from links */

  otherparity = opp_parity( parity );
  foralldir(dir) {
    odir = opp_dir(dir);
    
    /* start gather of up-adjoint link */

    tag0 = start_get( ahiggs, dir, parity );

    /* multiply adjoint here with up-link, for opp-parity */
    forparity(i,otherparity) {
      prefetch_matrix( &U[dir][i+1] );
      prefetch_adjoint( &ahiggs[i+1] );
      uncompress_adjmat(&ahiggs[i],&tmat1);
      mult_su3_an( &U[dir][i], &tmat1, &tmat2 );
      mult_su3_nn( &tmat2, &U[dir][i], &tmat1 );
      compress_adjmat( &tmat1, &adjstaple[i] );
    }

    tag1 = start_get( adjstaple, odir, parity );

    wait_get(tag0);

    /* multiply link with up-adjoint */
    if (dir == XUP) forparity(i,parity) {
      prefetch_matrix( &U[dir][i+1] );
      prefetch_adjoint( &ahiggs[nb(dir,i+1)] );
      uncompress_adjmat( &ahiggs[nb(dir,i)], &tmat1);
      mult_su3_nn( &U[dir][i], &tmat1, &tmat2 );
      mult_su3_na( &tmat2, &U[dir][i], &tmat1 );
      compress_adjmat( &tmat1, &adjstaple[i] );
    } else forparity(i,parity) {
      prefetch_matrix( &U[dir][i+1] );
      prefetch_adjoint( &ahiggs[nb(dir,i+1)] );
      prefetch_adjoint( &adjstaple[i+1] );
      uncompress_adjmat( &ahiggs[nb(dir,i)], &tmat1);
      mult_su3_nn( &U[dir][i], &tmat1, &tmat2 );
      mult_su3_na( &tmat2, &U[dir][i], &tmat1 );
      compress_adjmat( &tmat1, &tadj );
      add_adjmat( &adjstaple[i], &tadj, &adjstaple[i] );
    }

    wait_get(tag1);
    forparity(i,parity) {
      prefetch_adjoint( &adjstaple[i+1] );
      prefetch_adjoint( &adjstaple[nb(odir,i+1)] );
      add_adjmat( &adjstaple[i], &adjstaple[nb(odir,i)] , &adjstaple[i] );
    }
  }  
}


/******************************************************
 * needed for adjoint acceptance
 * note-first su3_matrix contains the 'local' adjoint
 * matrix uncompressed
 *
 * acc/rej with
 *
 *     a |-- u --> an
 *
 *     -2 beta tr (a u an u')
 *
 *****************************************************/

double act_gauge_adj(su3_matrix *a, su3_matrix *u,adjoint_matrix *b)
{
  su3_matrix am,m1,m2;

  uncompress_adjmat(b,&am);
  mult_su3_nn(a,u,&m1);
  mult_su3_na(u,&am,&m2);
  return((-2.0)*betaA*realtrace_su3(&m1,&m2));  /* m2' * m1 */
}

