Skip to content

Commit

Permalink
update_momenta.c: typo corrected
Browse files Browse the repository at this point in the history 
test/check_xchange.c broken, xchange_deri() needs an argument

clover.c: removed VOLUME parameter to init_sw_fields

clover_leaf.c: replaced uggly macros with proper inline functions

hybrid_update.c: renamed ini_momenta -> init_momenta

test/check_xchange.c: corrected missing argument to xchange_deri

cloverdet_monomial.c: removed non-even/odd branch

new xchange_deri version checked for 1-dim case
  • Loading branch information
@urbach
urbach committed Feb 10, 2012
1 parent 2c2aae0 commit d635b67
Show file tree
Hide file tree
Showing 13 changed files with 137 additions and 171 deletions.
9 changes: 5 additions & 4 deletions clover.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,7 +313,8 @@ void clover(const int ieo,
su3 ** sw1, ** sw_inv1;
su3 * _sw, *_sw_inv;

void init_sw_fields(const int V) {
void init_sw_fields() {
int V = VOLUME;
su3 * tmp;
static int sw_init = 0;

Expand All @@ -338,7 +339,7 @@ void init_sw_fields(const int V) {
}
sw[0] = sw1;
sw_inv[0] = sw_inv1;
for(int i = 1; i < VOLUME; i++) {
for(int i = 1; i < V; i++) {
sw[i] = sw[i-1]+3;
sw_inv[i] = sw_inv[i-1]+3;
}
Expand All @@ -350,15 +351,15 @@ void init_sw_fields(const int V) {
sw_inv[0][0] = _sw_inv;
# endif
tmp = sw[0][0];
for(int i = 0; i < VOLUME; i++) {
for(int i = 0; i < V; i++) {
for(int j = 0; j < 3; j++) {
sw[i][j] = tmp;
tmp = tmp+2;
}
}

tmp = sw_inv[0][0];
for(int i = 0; i < VOLUME; i++) {
for(int i = 0; i < V; i++) {
for(int j = 0; j < 3; j++) {
sw_inv[i][j] = tmp;
tmp = tmp+2;
Expand Down
2 changes: 1 addition & 1 deletion clover.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,6 @@ void Qsw_psi(spinor * const l, spinor * const k);
void Qsw_sq_psi(spinor * const l, spinor * const k);
void Msw_psi(spinor * const l, spinor * const k);
void H_eo_sw_inv_psi(spinor * const l, spinor * const k, const int ieo);
void init_sw_fields(const int V);
void init_sw_fields();

#endif
87 changes: 40 additions & 47 deletions clover_leaf.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -343,32 +343,35 @@ double six_det(complex a[6][6]) {
}

/*definitions needed for the functions sw_trace(int ieo) and sw_trace(int ieo)*/
#define _a_C(A, B, C) \
a[0+(A)][0+(B)]=(C).c00; \
a[0+(A)][1+(B)]=(C).c01; \
a[0+(A)][2+(B)]=(C).c02; \
a[1+(A)][0+(B)]=(C).c10; \
a[1+(A)][1+(B)]=(C).c11; \
a[1+(A)][2+(B)]=(C).c12; \
a[2+(A)][0+(B)]=(C).c20; \
a[2+(A)][1+(B)]=(C).c21; \
a[2+(A)][2+(B)]=(C).c22;

#define _C_a(A, B, C) \
(C).c00=a[0+(A)][0+(B)]; \
(C).c01=a[0+(A)][1+(B)]; \
(C).c02=a[0+(A)][2+(B)]; \
(C).c10=a[1+(A)][0+(B)]; \
(C).c11=a[1+(A)][1+(B)]; \
(C).c12=a[1+(A)][2+(B)]; \
(C).c20=a[2+(A)][0+(B)]; \
(C).c21=a[2+(A)][1+(B)]; \
(C).c22=a[2+(A)][2+(B)];
inline void populate_6x6_matrix(complex a[6][6], su3 * C, const int row, const int col) {
a[0+row][0+col] = C->c00;
a[0+row][1+col] = C->c01;
a[0+row][2+col] = C->c02;
a[1+row][0+col] = C->c10;
a[1+row][1+col] = C->c11;
a[1+row][2+col] = C->c12;
a[2+row][0+col] = C->c20;
a[2+row][1+col] = C->c21;
a[2+row][2+col] = C->c22;
return;
}

inline void get_3x3_block_matrix(su3 * C, complex a[6][6], const int row, const int col) {
C->c00 = a[0+row][0+col];
C->c01 = a[0+row][1+col];
C->c02 = a[0+row][2+col];
C->c10 = a[1+row][0+col];
C->c11 = a[1+row][1+col];
C->c12 = a[1+row][2+col];
C->c20 = a[2+row][0+col];
C->c21 = a[2+row][1+col];
C->c22 = a[2+row][2+col];
return;
}

double sw_trace(const int ieo) {
int i,x,icx,ioff;
su3 *w;
static su3 v2;
static su3 v;
static complex a[6][6];
static double tra;
static double ks,kc,tr,ts,tt;
Expand All @@ -385,14 +388,11 @@ double sw_trace(const int ieo) {
for(icx = ioff; icx < (VOLUME/2+ioff); icx++) {
x = g_eo2lexic[icx];
for(i=0;i<2;i++) {
w=&sw[x][0][i];
_a_C(0,0,*w);
w=&sw[x][1][i];
_a_C(0,3,*w);
_su3_dagger(v2,*w);
_a_C(3,0,v2);
w=&sw[x][2][i];
_a_C(3,3,*w);
populate_6x6_matrix(a, &sw[x][0][i], 0, 0);
populate_6x6_matrix(a, &sw[x][1][i], 0, 3);
_su3_dagger(v, sw[x][1][i]);
populate_6x6_matrix(a, &v, 3, 0);
populate_6x6_matrix(a, &sw[x][2][i], 3, 3);
tra = log(six_det(a));

tr=tra+kc;
Expand All @@ -417,8 +417,7 @@ double sw_trace(const int ieo) {
void sw_invert(const int ieo) {
int icx,ioff, err=0;
int i,x;
su3 *w;
static su3 v2;
static su3 v;
static complex a[6][6];
if(ieo==0) {
ioff=0;
Expand All @@ -431,14 +430,11 @@ void sw_invert(const int ieo) {
x = g_eo2lexic[icx];

for(i = 0; i < 2; i++) {
w=&sw[x][0][i];
_a_C(0,0,*w);
w=&sw[x][1][i];
_a_C(0,3,*w);
_su3_dagger(v2,*w);
_a_C(3,0,v2);
w=&sw[x][2][i];
_a_C(3,3,*w);
populate_6x6_matrix(a, &sw[x][0][i], 0, 0);
populate_6x6_matrix(a, &sw[x][1][i], 0, 3);
_su3_dagger(v, sw[x][1][i]);
populate_6x6_matrix(a, &v, 3, 0);
populate_6x6_matrix(a, &sw[x][2][i], 3, 3);

err = six_invert(a);
/* here we need to catch the error! */
Expand All @@ -448,12 +444,9 @@ void sw_invert(const int ieo) {
}

/* copy "a" back to sw_inv */
w=&sw_inv[x][0][i];
_C_a(0,0,*w);
w=&sw_inv[x][1][i];
_C_a(0,3,*w);
w=&sw_inv[x][2][i];
_C_a(3,3,*w);
get_3x3_block_matrix(&sw_inv[x][0][i], a, 0, 0);
get_3x3_block_matrix(&sw_inv[x][1][i], a, 0, 3);
get_3x3_block_matrix(&sw_inv[x][2][i], a, 3, 3);
}
}
return;
Expand Down
2 changes: 1 addition & 1 deletion clover_trlog_monomial.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ void clover_trlog_heatbath(const int id, hamiltonian_field_t * const hf) {
monomial * mnl = &monomial_list[id];
mnl->energy0 = 0.;

init_sw_fields(VOLUME);
init_sw_fields();
sw_term(hf->gaugefield, mnl->kappa, mnl->c_sw);
/*compute the contribution from the clover trlog term */
mnl->energy0 = -2.*sw_trace(EO);
Expand Down
171 changes: 71 additions & 100 deletions cloverdet_monomial.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,71 +57,58 @@ void cloverdet_derivative(const int id, hamiltonian_field_t * const hf) {
/* This factor 2 a missing factor 2 in trace_lambda */
(*mnl).forcefactor = 2.;

if(mnl->even_odd_flag) {
/*********************************************************************
*
* even/odd version
*
* This a term is det(\hat Q^2(\mu))
*
*********************************************************************/

g_mu = 0.;
boundary(mnl->kappa);

sw_term(hf->gaugefield, mnl->kappa, mnl->c_sw);
sw_invert(EO);

if(mnl->solver != CG && g_proc_id == 0) {
fprintf(stderr, "Bicgstab currently not implemented, using CG instead! (cloverdet_monomial.c)\n");
}

/* Invert Q_{+} Q_{-} */
/* X_o -> DUM_DERI+1 */
chrono_guess(g_spinor_field[DUM_DERI+1], mnl->pf, mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, mnl->csg_n, VOLUME/2, mnl->Qsq);
mnl->iter1 += cg_her(g_spinor_field[DUM_DERI+1], mnl->pf, mnl->maxiter, mnl->forceprec,
g_relative_precision_flag, VOLUME/2, mnl->Qsq);
chrono_add_solution(g_spinor_field[DUM_DERI+1], mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, &mnl->csg_n, VOLUME/2);

/* Y_o -> DUM_DERI */
mnl->Qm(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1]);

/* apply Hopping Matrix M_{eo} */
/* to get the even sites of X_e */
H_eo_sw_inv_psi(g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+1], EO);
/* \delta Q sandwitched by Y_o^\dagger and X_e */
deriv_Sb(OE, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+2], hf);

/* to get the even sites of Y_e */
H_eo_sw_inv_psi(g_spinor_field[DUM_DERI+3], g_spinor_field[DUM_DERI], EO);
/* \delta Q sandwitched by Y_e^\dagger and X_o */
deriv_Sb(EO, g_spinor_field[DUM_DERI+3], g_spinor_field[DUM_DERI+1], hf);

/* here comes the clover term... */
gamma5(g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+2], VOLUME/2);
sw_spinor(EO, g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+3]);

/* compute the contribution for the det-part */
gamma5(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI], VOLUME/2);
sw_spinor(OE, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1]);

sw_deriv(EO);
sw_all(hf, mnl->kappa, mnl->c_sw);
}
else {
/*********************************************************************
* non even/odd version not existent
*
* This term is det(Q^2 + \mu_1^2)
*
*********************************************************************/
if(g_proc_id == 0) {
fprintf(stderr, "only even/odd implemented right now! (cloverdet_monomial.c)\n");
}

/*********************************************************************
*
* even/odd version
*
* This a term is det(\hat Q^2(\mu))
*
*********************************************************************/

g_mu = 0.;
boundary(mnl->kappa);

sw_term(hf->gaugefield, mnl->kappa, mnl->c_sw);
sw_invert(EO);

if(mnl->solver != CG && g_proc_id == 0) {
fprintf(stderr, "Bicgstab currently not implemented, using CG instead! (cloverdet_monomial.c)\n");
}

/* Invert Q_{+} Q_{-} */
/* X_o -> DUM_DERI+1 */
chrono_guess(g_spinor_field[DUM_DERI+1], mnl->pf, mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, mnl->csg_n, VOLUME/2, mnl->Qsq);
mnl->iter1 += cg_her(g_spinor_field[DUM_DERI+1], mnl->pf, mnl->maxiter, mnl->forceprec,
g_relative_precision_flag, VOLUME/2, mnl->Qsq);
chrono_add_solution(g_spinor_field[DUM_DERI+1], mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, &mnl->csg_n, VOLUME/2);

/* Y_o -> DUM_DERI */
mnl->Qm(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1]);

/* apply Hopping Matrix M_{eo} */
/* to get the even sites of X_e */
H_eo_sw_inv_psi(g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+1], EO);
/* \delta Q sandwitched by Y_o^\dagger and X_e */
deriv_Sb(OE, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+2], hf);

/* to get the even sites of Y_e */
H_eo_sw_inv_psi(g_spinor_field[DUM_DERI+3], g_spinor_field[DUM_DERI], EO);
/* \delta Q sandwitched by Y_e^\dagger and X_o */
deriv_Sb(EO, g_spinor_field[DUM_DERI+3], g_spinor_field[DUM_DERI+1], hf);

/* here comes the clover term... */
gamma5(g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+2], VOLUME/2);
sw_spinor(EO, g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+3]);

/* compute the contribution for the det-part */
gamma5(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI], VOLUME/2);
sw_spinor(OE, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1]);

sw_deriv(EO);
sw_all(hf, mnl->kappa, mnl->c_sw);

g_mu = g_mu1;
boundary(g_kappa);

Expand All @@ -139,22 +126,17 @@ void cloverdet_heatbath(const int id, hamiltonian_field_t * const hf) {
mnl->iter0 = 0;
mnl->iter1 = 0;

init_sw_fields(VOLUME);
init_sw_fields();
sw_term(hf->gaugefield, mnl->kappa, mnl->c_sw);
sw_invert(EO);
if(mnl->even_odd_flag) {
random_spinor_field(g_spinor_field[2], VOLUME/2, mnl->rngrepro);
mnl->energy0 = square_norm(g_spinor_field[2], VOLUME/2, 1);

mnl->Qp(mnl->pf, g_spinor_field[2]);
chrono_add_solution(mnl->pf, mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, &mnl->csg_n, VOLUME/2);
}
else {
if(g_proc_id == 0) {
fprintf(stderr, "only even/odd implemented right now! (cloverdet_monomial.c)\n");
}
}
random_spinor_field(g_spinor_field[2], VOLUME/2, mnl->rngrepro);
mnl->energy0 = square_norm(g_spinor_field[2], VOLUME/2, 1);

mnl->Qp(mnl->pf, g_spinor_field[2]);
chrono_add_solution(mnl->pf, mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, &mnl->csg_n, VOLUME/2);

g_mu = g_mu1;
boundary(g_kappa);
if(g_proc_id == 0 && g_debug_level > 3) {
Expand All @@ -166,41 +148,30 @@ void cloverdet_heatbath(const int id, hamiltonian_field_t * const hf) {

double cloverdet_acc(const int id, hamiltonian_field_t * const hf) {
monomial * mnl = &monomial_list[id];
int save_iter = ITER_MAX_BCG;
int save_sloppy = g_sloppy_precision_flag;

sw_term(hf->gaugefield, mnl->kappa, mnl->c_sw);
sw_invert(EO);

g_mu = 0.;
boundary(mnl->kappa);
if(mnl->even_odd_flag) {

if(mnl->solver == CG) {
ITER_MAX_BCG = 0;
}
chrono_guess(g_spinor_field[2], mnl->pf, mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, mnl->csg_n, VOLUME/2, mnl->Qsq);
g_sloppy_precision_flag = 0;
mnl->iter0 = cg_her(g_spinor_field[2], mnl->pf, mnl->maxiter, mnl->accprec,
g_relative_precision_flag, VOLUME/2, mnl->Qsq);
mnl->Qm(g_spinor_field[2], g_spinor_field[2]);

g_sloppy_precision_flag = save_sloppy;
/* Compute the energy contr. from first field */
mnl->energy1 = square_norm(g_spinor_field[2], VOLUME/2, 1);
}
else {
if(g_proc_id == 0) {
fprintf(stderr, "only even/odd implemented right now! (cloverdet_monomial.c)\n");
}
}

chrono_guess(g_spinor_field[2], mnl->pf, mnl->csg_field, mnl->csg_index_array,
mnl->csg_N, mnl->csg_n, VOLUME/2, mnl->Qsq);
g_sloppy_precision_flag = 0;
mnl->iter0 = cg_her(g_spinor_field[2], mnl->pf, mnl->maxiter, mnl->accprec,
g_relative_precision_flag, VOLUME/2, mnl->Qsq);
mnl->Qm(g_spinor_field[2], g_spinor_field[2]);

g_sloppy_precision_flag = save_sloppy;
/* Compute the energy contr. from first field */
mnl->energy1 = square_norm(g_spinor_field[2], VOLUME/2, 1);

g_mu = g_mu1;
boundary(g_kappa);
if(g_proc_id == 0 && g_debug_level > 3) {
printf("called cloverdet_acc for id %d %d dH = %1.4e\n",
id, mnl->even_odd_flag, mnl->energy1 - mnl->energy0);
}
ITER_MAX_BCG = save_iter;
return(mnl->energy1 - mnl->energy0);
}
2 changes: 1 addition & 1 deletion hybrid_update.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ double moment_energy(su3adj ** const momenta) {
* with the gaussian distribution
*
**************************************/
double ini_momenta(const int repro, su3adj ** const momenta) {
double init_momenta(const int repro, su3adj ** const momenta) {

su3adj *xm;
int i, mu;
Expand Down
Loading

0 comments on commit d635b67

Please sign in to comment.