Skip to content
Browse files

This commit was manufactured by cvs2svn to create tag 'phmc_rc0'.

  • Loading branch information...
1 parent a77f33e commit 31b53e85dc2478bf620da7590f1b131e97e11ba0 cvs2svn committed Jan 18, 2007
Showing with 16,373 additions and 2,337 deletions.
  1. +4 −0 EVS.data
  2. +1 −0 INOUT.data
  3. +16 −6 Makefile.in
  4. +217 −54 Nondegenerate_Matrix.c
  5. +8 −0 Nondegenerate_Matrix.h
  6. +96 −0 PHMC-Makefile.in-PHMC
  7. +399 −0 Ptilde_nd.c
  8. +16 −0 Ptilde_nd.h
  9. +7 −0 README-for-PHMC
  10. +97 −0 Square_root_BR_roots.dat
  11. +391 −0 chebyshev_polynomial_nd.c
  12. +16 −0 chebyshev_polynomial_nd.h
  13. +258 −0 clenshaw_coef.c
  14. +6 −0 clenshaw_coef.h
  15. +5 −1 default_input_values.h
  16. +5 −3 eigenvalues_bi.c
  17. +1 −1 eigenvalues_bi.h
  18. +42 −2 global.h
  19. +197 −0 global.h-CARST
  20. +197 −0 global.h-TOM
  21. +11 −7 hmc.input
  22. +106 −31 hybrid_nondegenerate_update.c
  23. +10 −3 hybrid_nondegenerate_update.h
  24. +1 −0 hybrid_update.h
  25. +68 −0 init_chi_copy.c
  26. +11 −0 init_chi_copy.h
  27. +74 −0 init_chi_spinor_field.c
  28. +11 −0 init_chi_spinor_field.h
  29. +6 −1 max_eigenvalues_bi.c
  30. +1 −1 max_eigenvalues_bi.h
  31. +1 −0 normierungLocal.dat
  32. +854 −0 phmc_tm.c
  33. +843 −0 phmc_tm.c-COMPLETE
  34. +667 −0 phmc_tm.c-FULL
  35. +674 −0 phmc_tm.c-FULL-FIXED-POLY
  36. +640 −0 phmc_tm.c-ONLY-EV
  37. +2,412 −2,225 read_input.c
  38. +45 −1 read_input.l
  39. +98 −0 reweighting_factor_nd.c
  40. +8 −0 reweighting_factor_nd.h
  41. +5 −1 solver/Makefile.in
  42. +101 −0 solver/PHMC-Makefile.in-SOLVER
  43. +5 −0 solver/jdher_bi.c
  44. +6 −0 solver/pjdher_bi.c
  45. +40 −0 start.c
  46. +618 −0 start.c-NEW
  47. +578 −0 start.c-SAVE
  48. +1 −0 start.h
  49. +20 −0 start.h-NEW
  50. +19 −0 start.h-SAVE
  51. +611 −0 update_tm_nd.c
  52. +645 −0 update_tm_nd.c-CHECK-H-FINAL
  53. +680 −0 update_tm_nd.c-FOR-CHECKS
  54. +701 −0 update_tm_nd.c-FOR-CHECKS-COMPLETE
  55. +611 −0 update_tm_nd.c-GOOD
  56. +602 −0 update_tm_nd.c-NO-CORR-W-Qpoly
  57. +604 −0 update_tm_nd.c-NO-CORR-WITH-B
  58. +762 −0 update_tm_nd.c-SOME-INITIAL-CHECKS
  59. +636 −0 update_tm_nd.c-WITH-CORR-W-Qpoly
  60. +599 −0 update_tm_nd.c-WITH-Qpoly
  61. +9 −0 update_tm_nd.h
View
4 EVS.data
@@ -0,0 +1,4 @@
+ EV_min Ev_max n s_low s_max normalisation
+ 0.218281 3.401855 23 0.064000 1.000000 0.494939
+ 0.138196 2.587464 23 0.053000 1.000000 0.567508
+ 0.213994 2.519209 19 0.084000 1.000000 0.575145
View
1 INOUT.data
@@ -0,0 +1 @@
+ 0.031000 1.000000 0.494939
View
22 Makefile.in
@@ -27,27 +27,35 @@ COMPILE = ${CC} $(DEFS) $(INCLUDES) -o $@ ${CFLAGS} ${OPTARGS}
MODULES = init_gauge_field init_geometry_indices init_spinor_field \
init_bispinor_field init_moment_field init_gauge_tmp \
+ init_chi_spinor_field init_chi_copy \
read_input gamma linsolve xchange_field xchange_gauge \
xchange_deri expo hybrid_update observables start \
+ hybrid_nondegenerate_update \
geometry_eo linalg_eo ranlxd io Hopping_Matrix \
boundary deriv_Sb tm_operators getopt sighandler \
- mpi_init update_tm Hopping_Matrix_nocom \
+ mpi_init Hopping_Matrix_nocom \
+ update_tm_nd \
test/check_geometry test/check_xchange invert_eo \
measure_rectangles get_rectangle_staples \
derivative_psf ext_integrator polyakov_loop \
2mn_integrator get_staples update_backward_gauge \
- reweight_kappac
+ reweight_kappac reweighting_factor_nd \
+ chebyshev_polynomial_nd Ptilde_nd Nondegenerate_Matrix \
+ eigenvalues_bi max_eigenvalues_bi \
+### update_tm \
### eigenvalues_bi max_eigenvalues_bi Nondegenerate_Matrix \
-PROGRAMS = hmc_tm thermal_cycle_tm benchmark invert gwc2ildg ildg2gwc single2double double2single
+PROGRAMS = phmc_tm
+### hmc_tm
### test/test_eigenvalues
ALLOBJ = ${MODULES} ${PROGRAMS}
SUBDIRS = linalg solver
.SUFFIXES:
-all: Makefile all-recursive dep ${LINKLIBS} hmc_tm invert
+all: Makefile all-recursive dep ${LINKLIBS} phmc_tm
+### hmc_tm
### test/test_eigenvalues
.NOTPARALLEL:
@@ -74,10 +82,12 @@ compile-clean: compile-clean-recursive Makefile
rm -f *.o *.d
clean: clean-recursive Makefile
- rm -f hmc_tm hybrid *.o *.d
+ rm -f phmc_tm *.o *.d
+### rm -f hmc_tm hybrid phmc_tm *.o *.d
distclean: distclean-recursive Makefile
- rm -f hmc_tm hybrid *.o *.d *~ Makefile config.log config.status fixed_volume.h
+ rm -f phmc_tm *.o *.d *~ Makefile config.log config.status fixed_volume.h
+### rm -f hmc_tm hybrid phmc_tm *.o *.d *~ Makefile config.log config.status fixed_volume.h
.PHONY: all clean compile-clean distclean dep \
flex_read_input ${PROGRAMS} all-recursive \
View
271 Nondegenerate_Matrix.c
@@ -8,6 +8,9 @@
* see notes of R. Frezzoti and T. Chiarappa for more details *
* Author: Karl Jansen *
* Karl.Jansen@desy.de *
+ * *
+ * Adapted by Thomas Chiarappa <Thomas.Chiarappa@mib.infn.it> *
+ * *
**************************************************************/
#ifdef HAVE_CONFIG_H
@@ -22,12 +25,10 @@
#include "gamma.h"
/* in piu` */
#include "linsolve.h"
-/* fine in piu` */
+
#include "linalg_eo.h"
#include "Nondegenerate_Matrix.h"
-/* internal */
-
void mul_one_minus_imubar(spinor * const l, spinor * const k);
/******************************************
@@ -63,7 +64,6 @@ void mul_one_plus_imubar(spinor * const l, spinor * const k);
void QNon_degenerate(spinor * const l_strange, spinor * const l_charm,
spinor * const k_strange, spinor * const k_charm){
- /* double invmaxev=1./20.;*/
double nrm = 1./(1.+g_mubar*g_mubar-g_epsbar*g_epsbar);
/* Here the M_oe Mee^-1 M_eo implementation */
@@ -99,10 +99,8 @@ void QNon_degenerate(spinor * const l_strange, spinor * const l_charm,
gamma5(l_charm, l_charm, VOLUME/2);
/* At the end, the normalisation by the max. eigenvalue */
- /*
mul_r(l_strange, invmaxev, l_strange, VOLUME/2);
mul_r(l_charm, invmaxev, l_charm, VOLUME/2);
- */
}
/******************************************
@@ -125,7 +123,6 @@ void QNon_degenerate(spinor * const l_strange, spinor * const l_charm,
void QdaggerNon_degenerate(spinor * const l_strange, spinor * const l_charm,
spinor * const k_strange, spinor * const k_charm){
- /* double invmaxev=1./20.;*/
double nrm = 1./(1.+g_mubar*g_mubar-g_epsbar*g_epsbar);
/* Here the M_oe Mee^-1 M_eo implementation */
@@ -161,10 +158,8 @@ void QdaggerNon_degenerate(spinor * const l_strange, spinor * const l_charm,
gamma5(l_strange, l_strange, VOLUME/2);
/* At the end, the normalisation by the max. eigenvalue */
- /*
mul_r(l_charm, invmaxev, l_charm, VOLUME/2);
mul_r(l_strange, invmaxev, l_strange, VOLUME/2);
- */
}
@@ -188,10 +183,8 @@ void QdaggerNon_degenerate(spinor * const l_strange, spinor * const l_charm,
void Q_Qdagger_ND(spinor * const l_strange, spinor * const l_charm,
spinor * const k_strange, spinor * const k_charm){
- /* double invmaxev=1./20.;*/
double nrm = 1./(1.+g_mubar*g_mubar-g_epsbar*g_epsbar);
-
/* FIRST THE Qhat(2x2)^dagger PART*/
/* Here the M_oe Mee^-1 M_eo implementation */
@@ -206,8 +199,6 @@ void Q_Qdagger_ND(spinor * const l_strange, spinor * const l_charm,
mul_r(g_spinor_field[DUM_MATRIX+2], nrm, g_spinor_field[DUM_MATRIX+2], VOLUME/2);
mul_r(g_spinor_field[DUM_MATRIX+3], nrm, g_spinor_field[DUM_MATRIX+3], VOLUME/2);
- /* where nrm (= 1/(1+mu^2 -eps^2)) has been defined at the beginning of
- the subroutine */
Hopping_Matrix(OE, g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+2]);
Hopping_Matrix(OE, g_spinor_field[DUM_MATRIX+1], g_spinor_field[DUM_MATRIX+3]);
@@ -226,18 +217,23 @@ void Q_Qdagger_ND(spinor * const l_strange, spinor * const l_charm,
gamma5(g_spinor_field[DUM_MATRIX+2], g_spinor_field[DUM_MATRIX+4], VOLUME/2);
gamma5(g_spinor_field[DUM_MATRIX+3], g_spinor_field[DUM_MATRIX+5], VOLUME/2);
- /* At the end, the normalisation by the max. eigenvalue */
- /*
- mul_r(k_charm, invmaxev, g_spinor_field[DUM_MATRIX+2], VOLUME/2);
- mul_r(k_strange, invmaxev, g_spinor_field[DUM_MATRIX+3], VOLUME/2);
- */
+ /* The normalisation by the max. eigenvalue is done twice at the end */
+
+
+ /* We have to reassigin as follows to avoid overwriting */
+ /* Recall in fact that Q^hat = tau_1 Q tau_1 , hence */
+
+ /* ABOVE: dum_matrix+2 is l_charm goes to dum_matrix+6 :BELOW */
+ /* ABOVE: dum_matrix+3 is l_strange goes to dum_matrix+7 :BELOW */
+ assign(g_spinor_field[DUM_MATRIX+6], g_spinor_field[DUM_MATRIX+2], VOLUME/2);
+ assign(g_spinor_field[DUM_MATRIX+7], g_spinor_field[DUM_MATRIX+3], VOLUME/2);
/* AND THEN THE Qhat(2x2) PART */
/* Here the M_oe Mee^-1 M_eo implementation */
- Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX], k_strange);
- Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX+1], k_charm);
+ Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+7]);
+ Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX+1], g_spinor_field[DUM_MATRIX+6]);
mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+2], g_spinor_field[DUM_MATRIX]);
mul_one_plus_imubar(g_spinor_field[DUM_MATRIX+3], g_spinor_field[DUM_MATRIX+1]);
@@ -247,18 +243,16 @@ void Q_Qdagger_ND(spinor * const l_strange, spinor * const l_charm,
mul_r(g_spinor_field[DUM_MATRIX+2], nrm, g_spinor_field[DUM_MATRIX+2], VOLUME/2);
mul_r(g_spinor_field[DUM_MATRIX+3], nrm, g_spinor_field[DUM_MATRIX+3], VOLUME/2);
- /* where nrm (= 1/(1+mu^2 -eps^2)) has been defined at the beginning of
- the subroutine */
-
+
Hopping_Matrix(OE, l_strange, g_spinor_field[DUM_MATRIX+2]);
Hopping_Matrix(OE, l_charm, g_spinor_field[DUM_MATRIX+3]);
/* Here the M_oo implementation */
- mul_one_plus_imubar(g_spinor_field[DUM_MATRIX], k_strange);
- mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+1], k_charm);
+ mul_one_plus_imubar(g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+7]);
+ mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+1], g_spinor_field[DUM_MATRIX+6]);
- assign_add_mul_r(g_spinor_field[DUM_MATRIX], k_charm, -g_epsbar, VOLUME/2);
- assign_add_mul_r(g_spinor_field[DUM_MATRIX+1], k_strange, -g_epsbar, VOLUME/2);
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+6], -g_epsbar, VOLUME/2);
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX+1], g_spinor_field[DUM_MATRIX+7], -g_epsbar, VOLUME/2);
diff(l_strange, g_spinor_field[DUM_MATRIX], l_strange, VOLUME/2);
diff(l_charm, g_spinor_field[DUM_MATRIX+1], l_charm, VOLUME/2);
@@ -267,15 +261,147 @@ void Q_Qdagger_ND(spinor * const l_strange, spinor * const l_charm,
gamma5(l_strange, l_strange, VOLUME/2);
gamma5(l_charm, l_charm, VOLUME/2);
+
+ /* At the end, the normalisation by the max. eigenvalue */
+ /* Twice invmaxev since we consider here D Ddag !!! */
+ mul_r(l_charm, invmaxev*invmaxev, l_charm, VOLUME/2);
+ mul_r(l_strange, invmaxev*invmaxev, l_strange, VOLUME/2);
+
+}
+
+
+/******************************************
+ *
+ * This is the implementation of
+ *
+ * L_POLY_MIN_CCONST = M - z_k
+ *
+ * with M = Qhat(2x2) tau_1 and z_k \in Complex
+ *
+ *
+ * needed in the evaluation of the forces when
+ * the Polynomial approximation is used
+ *
+ *
+ * For details, see documentation and comments of the
+ * above mentioned routines
+ *
+ * k_charm and k_strange are the input fields
+ * l_* the output fields
+ *
+ * it acts only on the odd part or only
+ * on a half spinor
+ ******************************************/
+void L_POLY_MIN_CCONST(spinor * const l_strange, spinor * const l_charm,
+ spinor * const k_strange, spinor * const k_charm, const complex z){
+
+
+ int ix;
+ spinor *r, *s;
+ static su3_vector phi1;
+
+ double nrm = 1./(1.+g_mubar*g_mubar-g_epsbar*g_epsbar);
+
+
+ /* tau_1 inverts the k_charm <-> k_strange spinors */
+ /* Apply first Qhat(2x2) and finally substract the constant */
+
+ /* Here the M_oe Mee^-1 M_eo implementation */
+
+ Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX], k_charm);
+ Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX+1], k_strange);
+
+
+ mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+4], g_spinor_field[DUM_MATRIX]);
+ mul_one_plus_imubar(g_spinor_field[DUM_MATRIX+3], g_spinor_field[DUM_MATRIX+1]);
+
+
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX+4], g_spinor_field[DUM_MATRIX+1], g_epsbar, VOLUME/2);
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX+3], g_spinor_field[DUM_MATRIX], g_epsbar, VOLUME/2);
+
+ mul_r(g_spinor_field[DUM_MATRIX+4], nrm, g_spinor_field[DUM_MATRIX+4], VOLUME/2);
+ mul_r(g_spinor_field[DUM_MATRIX+3], nrm, g_spinor_field[DUM_MATRIX+3], VOLUME/2);
+ /* where nrm (= 1/(1+mu^2 -eps^2)) has been defined at the beginning of
+ the subroutine */
+
+
+ Hopping_Matrix(OE, l_strange, g_spinor_field[DUM_MATRIX+4]);
+ Hopping_Matrix(OE, l_charm, g_spinor_field[DUM_MATRIX+3]);
+
+
+ /* Here the M_oo implementation */
+ mul_one_plus_imubar(g_spinor_field[DUM_MATRIX], k_charm);
+ mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+1], k_strange);
+
+
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX], k_strange, -g_epsbar, VOLUME/2);
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX+1], k_charm, -g_epsbar, VOLUME/2);
+
+
+ diff(l_strange, g_spinor_field[DUM_MATRIX], l_strange, VOLUME/2);
+ diff(l_charm, g_spinor_field[DUM_MATRIX+1], l_charm, VOLUME/2);
+
+
+ /* and finally the gamma_5 multiplication */
+ gamma5(l_strange, l_strange, VOLUME/2);
+ gamma5(l_charm, l_charm, VOLUME/2);
+
/* At the end, the normalisation by the max. eigenvalue */
- /*
- mul_r(l_charm, invmaxev, l_charm, VOLUME/2);
mul_r(l_strange, invmaxev, l_strange, VOLUME/2);
+ mul_r(l_charm, invmaxev, l_charm, VOLUME/2);
+
+ /*
+ printf(" IN UP: %f %f \n", l_strange[0].s2.c1.re, l_strange[0].s2.c1.im);
+ printf(" IN DN: %f %f \n", l_charm[0].s2.c1.re, l_charm[0].s2.c1.im);
*/
+ /* AND FINALLY WE SUBSTRACT THE C-CONSTANT */
+
+
+ /************ loop over all lattice sites ************/
+ for(ix = 0; ix < (VOLUME/2); ix++){
+
+ r=l_strange + ix;
+ s=k_strange + ix;
+
+ _complex_times_vector(phi1, z, (*s).s0);
+ _vector_sub_assign((*r).s0, phi1);
+ _complex_times_vector(phi1, z, (*s).s1);
+ _vector_sub_assign((*r).s1, phi1);
+ _complex_times_vector(phi1, z, (*s).s2);
+ _vector_sub_assign((*r).s2, phi1);
+ _complex_times_vector(phi1, z, (*s).s3);
+ _vector_sub_assign((*r).s3, phi1);
+
+ r=l_charm + ix;
+ s=k_charm + ix;
+
+ _complex_times_vector(phi1, z, (*s).s0);
+ _vector_sub_assign((*r).s0, phi1);
+ _complex_times_vector(phi1, z, (*s).s1);
+ _vector_sub_assign((*r).s1, phi1);
+ _complex_times_vector(phi1, z, (*s).s2);
+ _vector_sub_assign((*r).s2, phi1);
+ _complex_times_vector(phi1, z, (*s).s3);
+ _vector_sub_assign((*r).s3, phi1);
+ }
+
+ /*
+ printf(" IN 2 UP: %f %f \n", l_strange[0].s2.c1.re, l_strange[0].s2.c1.im);
+ printf(" IN 2 DN: %f %f \n", l_charm[0].s2.c1.re, l_charm[0].s2.c1.im);
+ */
+
+ /* Finally, we multiply by the constant Cpol */
+ /* which renders the polynomial in monomials */
+ /* identical to the polynomial a la clenshaw */;
+ mul_r(l_strange, Cpol, l_strange, VOLUME/2);
+ mul_r(l_charm, Cpol, l_charm, VOLUME/2);
+
}
+
+
/******************************************
*
* This is the same implementation as above of
@@ -296,7 +422,6 @@ void Q_Qdagger_ND(spinor * const l_strange, spinor * const l_charm,
******************************************/
void Q_Qdagger_ND_BI(bispinor * const bisp_l, bispinor * const bisp_k){
- /* double invmaxev=1./20.;*/
double nrm = 1./(1.+g_mubar*g_mubar-g_epsbar*g_epsbar);
@@ -324,7 +449,6 @@ void Q_Qdagger_ND_BI(bispinor * const bisp_l, bispinor * const bisp_k){
/* CREATE 2 SPINORS OUT OF 1 (INPUT) BISPINOR */
decompact(&k_strange[0], &k_charm[0], &bisp_k[0]);
-
/* FIRST THE Qhat(2x2)^dagger PART*/
/* Here the M_oe Mee^-1 M_eo implementation */
@@ -359,23 +483,23 @@ void Q_Qdagger_ND_BI(bispinor * const bisp_l, bispinor * const bisp_k){
gamma5(g_spinor_field[DUM_MATRIX+2], g_spinor_field[DUM_MATRIX+4], VOLUME/2);
gamma5(g_spinor_field[DUM_MATRIX+3], g_spinor_field[DUM_MATRIX+5], VOLUME/2);
- /* At the end, the normalisation by the max. eigenvalue */
- /*
- mul_r(k_charm, invmaxev, g_spinor_field[DUM_MATRIX+2], VOLUME/2);
- mul_r(k_strange, invmaxev, g_spinor_field[DUM_MATRIX+3], VOLUME/2);
- */
-
+ /* The normalisation by the max. eigenvalue is done twice at the end */
- /* !!! I HAVE REPLACED THE ABOVE LINES BY ASSIGNMENTS !!! */
- assign(k_charm, g_spinor_field[DUM_MATRIX+2], VOLUME/2);
- assign(k_strange, g_spinor_field[DUM_MATRIX+3], VOLUME/2);
+
+ /* We have to reassigin as follows to avoid overwriting */
+ /* Recall in fact that Q^hat = tau_1 Q tau_1 , hence */
+
+ /* ABOVE: dum_matrix+2 is l_charm goes to dum_matrix+6 :BELOW */
+ /* ABOVE: dum_matrix+3 is l_strange goes to dum_matrix+7 :BELOW */
+ assign(g_spinor_field[DUM_MATRIX+6], g_spinor_field[DUM_MATRIX+2], VOLUME/2);
+ assign(g_spinor_field[DUM_MATRIX+7], g_spinor_field[DUM_MATRIX+3], VOLUME/2);
/* AND THEN THE Qhat(2x2) PART */
/* Here the M_oe Mee^-1 M_eo implementation */
- Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX], k_strange);
- Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX+1], k_charm);
+ Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+7]);
+ Hopping_Matrix(EO, g_spinor_field[DUM_MATRIX+1], g_spinor_field[DUM_MATRIX+6]);
mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+2], g_spinor_field[DUM_MATRIX]);
mul_one_plus_imubar(g_spinor_field[DUM_MATRIX+3], g_spinor_field[DUM_MATRIX+1]);
@@ -392,11 +516,11 @@ void Q_Qdagger_ND_BI(bispinor * const bisp_l, bispinor * const bisp_k){
Hopping_Matrix(OE, l_charm, g_spinor_field[DUM_MATRIX+3]);
/* Here the M_oo implementation */
- mul_one_plus_imubar(g_spinor_field[DUM_MATRIX], k_strange);
- mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+1], k_charm);
+ mul_one_plus_imubar(g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+7]);
+ mul_one_minus_imubar(g_spinor_field[DUM_MATRIX+1], g_spinor_field[DUM_MATRIX+6]);
- assign_add_mul_r(g_spinor_field[DUM_MATRIX], k_charm, -g_epsbar, VOLUME/2);
- assign_add_mul_r(g_spinor_field[DUM_MATRIX+1], k_strange, -g_epsbar, VOLUME/2);
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+6], -g_epsbar, VOLUME/2);
+ assign_add_mul_r(g_spinor_field[DUM_MATRIX+1], g_spinor_field[DUM_MATRIX+7], -g_epsbar, VOLUME/2);
diff(l_strange, g_spinor_field[DUM_MATRIX], l_strange, VOLUME/2);
diff(l_charm, g_spinor_field[DUM_MATRIX+1], l_charm, VOLUME/2);
@@ -406,10 +530,10 @@ void Q_Qdagger_ND_BI(bispinor * const bisp_l, bispinor * const bisp_k){
gamma5(l_charm, l_charm, VOLUME/2);
/* At the end, the normalisation by the max. eigenvalue */
- /*
- mul_r(l_charm, invmaxev, l_charm, VOLUME/2);
- mul_r(l_strange, invmaxev, l_strange, VOLUME/2);
- */
+ /* Twice invmaxev since we consider here D Ddag !!! */
+ mul_r(l_charm, invmaxev*invmaxev, l_charm, VOLUME/2);
+ mul_r(l_strange, invmaxev*invmaxev, l_strange, VOLUME/2);
+
/* !!! I HAVE REPLACED THE ABOVE LINES BY ASSIGNMENTS !!! */
@@ -430,14 +554,51 @@ void Q_Qdagger_ND_BI(bispinor * const bisp_l, bispinor * const bisp_k){
free(l_charm);
#endif
+}
+
+
+/******************************************
+ *
+ * This is the implementation of
+ *
+ * (M_{ee}^\pm)^{-1}M_{eo}
+ *
+ * see documentation for details
+ * k is the number of the input field
+ * l is the number of the output field
+ *
+ * it acts only on the odd part or only
+ * on a half spinor
+ ******************************************/
+void H_eo_ND(spinor * const l_strange, spinor * const l_charm,
+ spinor * const k_strange, spinor * const k_charm,
+ const int ieo){
+
+ double nrm = 1./(1.+g_mubar*g_mubar-g_epsbar*g_epsbar);
+
+
+ /* recall: strange <-> up while charm <-> dn */
+ Hopping_Matrix(ieo, g_spinor_field[DUM_MATRIX], k_strange);
+ Hopping_Matrix(ieo, g_spinor_field[DUM_MATRIX+1], k_charm);
+
+ mul_one_minus_imubar(l_strange, g_spinor_field[DUM_MATRIX+1]);
+ mul_one_plus_imubar(l_charm, g_spinor_field[DUM_MATRIX]);
+
+ assign_add_mul_r(l_strange, g_spinor_field[DUM_MATRIX], g_epsbar, VOLUME/2);
+ assign_add_mul_r(l_charm, g_spinor_field[DUM_MATRIX+1], g_epsbar, VOLUME/2);
+
+ mul_r(l_strange, nrm, l_strange, VOLUME/2);
+ mul_r(l_charm, nrm, l_charm, VOLUME/2);
}
+
+
+
void Q_test_epsilon(spinor * const l_strange, spinor * const l_charm,
spinor * const k_strange, spinor * const k_charm){
- /* double invmaxev=1./20.; */
double nrm = 1./(1.+g_mubar*g_mubar-g_epsbar*g_epsbar);
/* Here the M_oe Mee^-1 M_eo implementation */
@@ -458,10 +619,8 @@ void Q_test_epsilon(spinor * const l_strange, spinor * const l_charm,
gamma5(l_charm, l_charm, VOLUME/2);
/* At the end, the normalisation by the max. eigenvalue */
- /*
mul_r(l_charm, invmaxev, l_charm, VOLUME/2);
mul_r(l_strange, invmaxev, l_strange, VOLUME/2);
- */
}
@@ -522,5 +681,9 @@ void mul_one_plus_imubar(spinor * const l, spinor * const k){
}
-
static char const rcsid[] = "$Id$";
+
+
+
+
+
View
8 Nondegenerate_Matrix.h
@@ -13,6 +13,14 @@ void Q_Qdagger_ND(spinor * const l_strange, spinor * const l_charm,
void Q_Qdagger_ND_BI(bispinor * const bisp_l, bispinor * const bisp_k);
+void L_POLY_MIN_CCONST(spinor * const l_strange, spinor * const l_charm,
+ spinor * const k_strange, spinor * const k_charm,
+ const complex z);
+
+void H_eo_ND(spinor * const l_strange, spinor * const l_charm,
+ spinor * const k_strange, spinor * const k_charm,
+ const int ieo);
+
void Q_test_epsilon(spinor * const l_strange, spinor * const l_charm,
spinor * const k_strange, spinor * const k_charm);
View
96 PHMC-Makefile.in-PHMC
@@ -0,0 +1,96 @@
+srcdir = @srcdir@
+top_srcdir = @top_srcdir@
+top_builddir = .
+abs_top_builddir = @abs_top_builddir@
+builddir = @builddir@
+subdir = .
+
+
+CC = @CC@
+CCDEP = @CCDEP@
+CFLAGS = @CFLAGS@
+LDFLAGS = @LDFLAGS@
+DEPFLAGS = @DEPFLAGS@
+CCLD = $(CC)
+LEX = @LEX@
+AUTOCONF = @AUTOCONF@
+LIBS = @LIBS@
+SHELL = @SHELL@
+OPTARGS = @OPTARGS@
+DEFS = @DEFS@
+
+INCLUDES = @INCLUDES@
+LINK = $(CCLD) -o $@ ${LDFLAGS} ${OPTARGS}
+LINKLIBS = ${top_builddir}/linalg/liblinalg.a \
+ ${top_builddir}/solver/libsolver.a
+COMPILE = ${CC} $(DEFS) $(INCLUDES) -o $@ ${CFLAGS} ${OPTARGS}
+
+MODULES = init_gauge_field init_geometry_indices init_spinor_field \
+ init_bispinor_field init_moment_field init_gauge_tmp \
+ init_chi_spinor_field init_chi_copy \
+ read_input gamma linsolve xchange_field xchange_gauge \
+ xchange_deri expo hybrid_update observables start \
+ hybrid_nondegenerate_update \
+ geometry_eo linalg_eo ranlxd io Hopping_Matrix \
+ boundary deriv_Sb tm_operators getopt sighandler \
+ mpi_init Hopping_Matrix_nocom \
+ update_tm_nd \
+ test/check_geometry test/check_xchange invert_eo \
+ measure_rectangles get_rectangle_staples \
+ derivative_psf ext_integrator polyakov_loop \
+ 2mn_integrator get_staples update_backward_gauge \
+ reweight_kappac \
+ chebyshev_polynomial_nd Ptilde_nd Nondegenerate_Matrix \
+ eigenvalues_bi max_eigenvalues_bi \
+
+### update_tm \
+### eigenvalues_bi max_eigenvalues_bi Nondegenerate_Matrix \
+
+PROGRAMS = phmc_tm
+### hmc_tm
+### test/test_eigenvalues
+ALLOBJ = ${MODULES} ${PROGRAMS}
+SUBDIRS = linalg solver
+
+.SUFFIXES:
+
+all: Makefile all-recursive dep ${LINKLIBS} phmc_tm
+### hmc_tm
+### test/test_eigenvalues
+
+.NOTPARALLEL:
+
+-include $(addsuffix .d,$(ALLOBJ))
+
+include ${top_srcdir}/Makefile.global
+
+flex_read_input: read_input.l
+ ${LEX} -i -t $< > read_input.c
+
+${addsuffix .o, ${MODULES}}: %.o: ${srcdir}/%.c %.d Makefile
+ ${COMPILE} -c $<
+
+${addsuffix .o, ${PROGRAMS}}: %.o: ${srcdir}/%.c %.d Makefile
+ ${COMPILE} -c $<
+
+${PROGRAMS}: %: %.o ${addsuffix .o, ${MODULES}} ${LINKLIBS}
+ ${LINK} ${addsuffix .o, ${MODULES}} $@.o $(LIBS)
+
+dep: $(addsuffix .d,$(ALLOBJ))
+
+compile-clean: compile-clean-recursive Makefile
+ rm -f *.o *.d
+
+clean: clean-recursive Makefile
+ rm -f phmc_tm *.o *.d
+### rm -f hmc_tm hybrid phmc_tm *.o *.d
+
+distclean: distclean-recursive Makefile
+ rm -f phmc_tm *.o *.d *~ Makefile config.log config.status fixed_volume.h
+### rm -f hmc_tm hybrid phmc_tm *.o *.d *~ Makefile config.log config.status fixed_volume.h
+
+.PHONY: all clean compile-clean distclean dep \
+ flex_read_input ${PROGRAMS} all-recursive \
+ all-debug-recursive all-profile-recursive \
+ clean-recursive distclean-recursive \
+ compile-clean-recursive
View
399 Ptilde_nd.c
@@ -0,0 +1,399 @@
+/* $Id$ */
+
+#ifdef HAVE_CONFIG_H
+# include<config.h>
+#endif
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include "global.h"
+#include "linsolve.h"
+#include "linalg_eo.h"
+#include "start.h"
+#include "tm_operators.h"
+#include "Nondegenerate_Matrix.h"
+#include "chebyshev_polynomial_nd.h"
+#include "Ptilde_nd.h"
+
+
+#define PI 3.141592653589793
+
+
+double func_tilde(double u, double exponent){
+
+ double ff=0.0;
+ double d=0,ddd=0, sv, z, z2;
+ int j;
+ double res=0.0;
+
+ z = (2.0*u - cheb_evmin - cheb_evmax)/(double)(cheb_evmax - cheb_evmin);
+ z2 = 2.0*z;
+
+ for(j=dop_n_cheby-1; j>=1; j--){
+ sv = d;
+ d = z2*d - ddd + dop_cheby_coef[j];
+ ddd = sv;
+ }
+
+ res = z*d - ddd + 0.5*dop_cheby_coef[0];
+
+ ff = (double)(res * sqrt(u));
+
+ return(pow(ff,exponent));
+}
+
+void Ptilde_cheb_coefs(double aa, double bb, double dd[], int n, double exponent){
+ int k,j;
+ double fac,bpa,bma,*f;
+ double inv_n;
+
+ inv_n=1./(double)n;
+ f=calloc(n,sizeof(double));/*vector(0,n-1);*/
+ if(g_proc_id == g_stdio_proc){
+ printf("\n hello in PTILDE-chebyshev_polynomial\n");
+ printf("n= %d inv_n=%e \n",n,inv_n);
+ printf("allocation !!!\n");
+ }
+ fflush(stdout);
+ bma=0.5*(bb-aa);
+ bpa=0.5*(bb+aa);
+ for (k=0;k<n;k++) {
+ double y=cos(PI*(k+0.5)*inv_n);
+ f[k]=func_tilde(y*bma+bpa,exponent);
+ }
+ fac=2.0*inv_n;
+ for (j=0;j<n;j++) {
+ double sum=0.0;
+ for (k=0;k<n;k++)
+ sum += f[k]*cos(PI*j*(k+0.5)*inv_n);
+ dd[j]=fac*sum;
+ }
+ free(f);
+}
+#undef PI
+
+
+/****************************************************************************
+ *
+ * computation of the second poplynomial, Ptilde, on a vector
+ * by using the chebyshev approximation for the function ()^1/4
+ * subtraction of low-lying eigenvalues is not yet implemented for this
+ *
+ **************************************************************************/
+
+void Poly_tilde_ND(spinor *R_s, spinor *R_c, double *dd, int n,
+ spinor *S_s, spinor *S_c){
+
+ int j;
+ double fact1, fact2, temp, temp1, temp2, temp3, temp4;
+
+ spinor *svs_=NULL, *svs=NULL, *ds_=NULL, *ds=NULL, *dds_=NULL, *dds=NULL,
+ *auxs_=NULL, *auxs=NULL, *aux2s_=NULL, *aux2s=NULL, *aux3s_=NULL,
+ *aux3s=NULL;
+ spinor *svc_=NULL, *svc=NULL, *dc_=NULL, *dc=NULL, *ddc_=NULL,
+ *ddc=NULL, *auxc_=NULL, *auxc=NULL, *aux2c_=NULL, *aux2c=NULL,
+ *aux3c_=NULL, *aux3c=NULL;
+
+
+#if ( defined SSE || defined SSE2 )
+ svs_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ svs = (spinor *)(((unsigned int)(svs_)+ALIGN_BASE)&~ALIGN_BASE);
+ ds_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ ds = (spinor *)(((unsigned int)(ds_)+ALIGN_BASE)&~ALIGN_BASE);
+ dds_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ dds = (spinor *)(((unsigned int)(dds_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxs_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ auxs = (spinor *)(((unsigned int)(auxs_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2s_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux2s = (spinor *)(((unsigned int)(aux2s_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux3s_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux3s = (spinor *)(((unsigned int)(aux3s_)+ALIGN_BASE)&~ALIGN_BASE);
+ svc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ svc = (spinor *)(((unsigned int)(svc_)+ALIGN_BASE)&~ALIGN_BASE);
+ dc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ dc = (spinor *)(((unsigned int)(dc_)+ALIGN_BASE)&~ALIGN_BASE);
+ ddc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ ddc = (spinor *)(((unsigned int)(ddc_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ auxc = (spinor *)(((unsigned int)(auxc_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2c_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux2c = (spinor *)(((unsigned int)(aux2c_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux3c_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux3c = (spinor *)(((unsigned int)(aux3c_)+ALIGN_BASE)&~ALIGN_BASE);
+#else
+ svs_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ svs = svs_;
+ ds_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ ds = ds_;
+ dds_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ dds = dds_;
+ auxs_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ auxs = auxs_;
+ aux2s_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux2s = aux2s_;
+ aux3s_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux3s = aux3s_;
+ svc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ svc = svc_;
+ dc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ dc = dc_;
+ ddc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ ddc = ddc_;
+ auxc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ auxc = auxc_;
+ aux2c_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux2c = aux2c_;
+ aux3c_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux3c = aux3c_;
+#endif
+
+
+ fact1=4/(cheb_evmax-cheb_evmin);
+ fact2=-2*(cheb_evmax+cheb_evmin)/(cheb_evmax-cheb_evmin);
+
+ zero_spinor_field(&ds[0],VOLUME/2);
+ zero_spinor_field(&dds[0],VOLUME/2);
+ zero_spinor_field(&dc[0],VOLUME/2);
+ zero_spinor_field(&ddc[0],VOLUME/2);
+
+ /* sub_low_ev(&aux3[0], &S[0]); */
+ assign(&aux3s[0], &S_s[0],VOLUME/2);
+ assign(&aux3c[0], &S_c[0],VOLUME/2);
+
+ /* Use the Clenshaw's recursion for the Chebysheff polynomial */
+ for (j=n-1; j>=1; j--) {
+ assign(&svs[0],&ds[0],VOLUME/2);
+ assign(&svc[0],&dc[0],VOLUME/2);
+
+ /*
+ if ( (j%10) == 0 ) {
+ sub_low_ev(&aux[0], &d[0]);
+ }
+ else { */
+ assign(&auxs[0], &ds[0], VOLUME/2);
+ assign(&auxc[0], &dc[0], VOLUME/2);
+ /* } */
+
+
+ Q_Qdagger_ND(&R_s[0], &R_c[0], &auxs[0], &auxc[0]);
+
+ temp1=-1.0;
+ temp2=dd[j];
+ assign_mul_add_mul_add_mul_add_mul_r(&ds[0] , &R_s[0], &dds[0], &aux3s[0], fact2, fact1, temp1, temp2,VOLUME/2);
+ assign_mul_add_mul_add_mul_add_mul_r(&dc[0] , &R_c[0], &ddc[0], &aux3c[0], fact2, fact1, temp1, temp2,VOLUME/2);
+ assign(&dds[0], &svs[0],VOLUME/2);
+ assign(&ddc[0], &svc[0],VOLUME/2);
+ }
+
+ /* sub_low_ev(&R[0],&d[0]); */
+ assign(&R_s[0], &ds[0],VOLUME/2);
+ assign(&R_c[0], &dc[0],VOLUME/2);
+
+
+ Q_Qdagger_ND(&auxs[0], &auxc[0], &R_s[0], &R_c[0]);
+
+ temp1=-1.0;
+ temp2=dd[0]/2;
+ temp3=fact1/2;
+ temp4=fact2/2;
+ assign_mul_add_mul_add_mul_add_mul_r(&auxs[0], &ds[0], &dds[0], &aux3s[0], temp3, temp4, temp1, temp2,VOLUME/2);
+ assign_mul_add_mul_add_mul_add_mul_r(&auxc[0], &dc[0], &ddc[0], &aux3c[0], temp3, temp4, temp1, temp2,VOLUME/2);
+ assign(&R_s[0], &auxs[0],VOLUME/2);
+ assign(&R_c[0], &auxc[0],VOLUME/2);
+
+ /* addproj_q_invsqrt(&R[0], &S[0]); */
+
+ /*
+#ifndef _SOLVER_OUTPUT
+ if(g_proc_id == g_stdio_proc){
+ printf("Order of Chebysheff approximation = %d\n",j);
+ fflush( stdout);};
+#endif
+ */
+
+
+ free(svs_);
+ free(ds_);
+ free(dds_);
+ free(auxs_);
+ free(aux2s_);
+ free(aux3s_);
+ free(svc_);
+ free(dc_);
+ free(ddc_);
+ free(auxc_);
+ free(aux2c_);
+ free(aux3c_);
+
+}
+
+
+double chebtilde_eval(int M, double *dd, double s){
+
+ double d=0,ddd=0, sv, z, z2, res;
+ int j;
+
+ z = (2.0*s - cheb_evmin - cheb_evmax)/(double)(cheb_evmax - cheb_evmin);
+ z2 = 2.0*z;
+
+ for(j=M-1; j>=1; j--){
+ sv = d;
+ d = z2*d - ddd + dd[j];
+ ddd = sv;
+ }
+
+ res = z*d - ddd + 0.5*dd[0];
+
+ return(res);
+}
+
+/**************************************************************************
+ *
+ * The externally accessible function is
+ *
+ * void degree_of_Ptilde(void)
+ * Computation of (QdaggerQ)^1/4
+ * by using the chebyshev approximation for the function ()^1/4
+ *
+ * Author: Thomas Chiarappa <Thomas.Chiarappa@mib.infn.it> Mai 2006
+ *
+*****************************************************************************/
+
+
+
+void degree_of_Ptilde(){
+ int i, j;
+ double temp, temp2;
+ static int ini=0;
+
+ double sum=0.0;
+
+ spinor *ss=NULL, *ss_=NULL, *sc=NULL, *sc_=NULL;
+ spinor *auxs=NULL, *auxs_=NULL, *auxc=NULL, *auxc_=NULL;
+ spinor *aux2s=NULL, *aux2s_=NULL, *aux2c=NULL, *aux2c_=NULL;
+
+ if(ini==0){
+ ptilde_cheby_coef = calloc(NTILDE_CHEBYMAX,sizeof(double));
+ ini=1;
+ }
+
+#if ( defined SSE || defined SSE2 || defined SSE3)
+ ss_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ auxs_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ aux2s_= calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ sc_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ auxc_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ aux2c_= calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+
+ ss = (spinor *)(((unsigned int)(ss_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxs = (spinor *)(((unsigned int)(auxs_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2s = (spinor *)(((unsigned int)(aux2s_)+ALIGN_BASE)&~ALIGN_BASE);
+ sc = (spinor *)(((unsigned int)(sc_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxc = (spinor *)(((unsigned int)(auxc_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2c = (spinor *)(((unsigned int)(aux2c_)+ALIGN_BASE)&~ALIGN_BASE);
+
+#else
+ ss =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ auxs =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ aux2s=calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ sc =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ auxc =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ aux2c=calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+#endif
+
+
+ Ptilde_cheb_coefs(cheb_evmin, cheb_evmax, ptilde_cheby_coef, NTILDE_CHEBYMAX, -1.0);
+
+ random_spinor_field(ss,VOLUME/2);
+ random_spinor_field(sc,VOLUME/2);
+
+ if(g_proc_id == g_stdio_proc){
+ printf(" \n In Ptilde: EVmin = %f EVmax = %f\n", cheb_evmin, cheb_evmax);
+ printf("\n determine the degree of the polynomial : Stop=%e \n", g_acc_Ptilde);
+ fflush(stdout);
+ }
+
+ ptilde_n_cheby = dop_n_cheby;
+
+ for(i = 0;i < 100 ; i++){
+
+ if (ptilde_n_cheby > NTILDE_CHEBYMAX) {
+ if(g_proc_id == g_stdio_proc){
+ printf("Error: n_cheby=%d > NTILDE_CHEBYMAX=%d\n",ptilde_n_cheby,NTILDE_CHEBYMAX);
+ printf("Increase n_chebymax\n");
+ }
+ errorhandler(35,"degree_of_polynomial");
+ }
+
+
+ /* Ptilde P S P Ptilde X - X */
+ Poly_tilde_ND(&auxs[0], &auxc[0], ptilde_cheby_coef, ptilde_n_cheby, &ss[0], &sc[0]);
+ QdaggerQ_poly(&aux2s[0], &aux2c[0], dop_cheby_coef, dop_n_cheby, &auxs[0], &auxc[0]);
+ Q_Qdagger_ND(&auxs[0], &auxc[0], &aux2s[0], &aux2c[0]);
+ QdaggerQ_poly(&aux2s[0], &aux2c[0], dop_cheby_coef, dop_n_cheby, &auxs[0], &auxc[0]);
+ Poly_tilde_ND(&auxs[0], &auxc[0], ptilde_cheby_coef, ptilde_n_cheby, &aux2s[0], &aux2c[0]);
+
+
+ diff(&aux2s[0],&auxs[0],&ss[0],VOLUME/2);
+ temp=square_norm(&aux2s[0],VOLUME/2)/square_norm(&ss[0],VOLUME/2)/4.0;
+
+ diff(&aux2c[0],&auxc[0],&sc[0],VOLUME/2);
+ temp2=square_norm(&aux2c[0],VOLUME/2)/square_norm(&sc[0],VOLUME/2)/4.0;
+
+ if(g_epsbar == 0){
+ temp2 = 0.0;
+ }
+ if(g_proc_id == g_stdio_proc) {
+ printf("At n=%d || differences ||^2 : UP=%e DN=%e \n", ptilde_n_cheby, temp, temp2);
+ }
+
+
+ sum=0;
+ for(j=ptilde_n_cheby; j<NTILDE_CHEBYMAX; j++){
+ sum += fabs(ptilde_cheby_coef[j]);
+ }
+ if(g_proc_id == g_stdio_proc) printf(" Sum remaining | d_n | = %e\n", sum);
+
+ /* if(temp < g_acc_ptilde && temp2 < g_acc_ptilde){ */ /* break; */
+
+ if(sum < g_acc_Ptilde){
+ if(g_proc_id == g_stdio_proc){
+ printf("\n Achieved Accuracies for Ptilde : Stop=%e \n", g_acc_Ptilde);
+ printf(" Uniform: Sum |d_n|=%e \n", sum);
+ printf(" RND: || (Ptilde P S P Ptilde - 1)X ||^2 / || 2X ||^2 : UP=%e DN=%e \n",temp, temp2);
+ }
+
+ temp = chebtilde_eval(ptilde_n_cheby, ptilde_cheby_coef, cheb_evmin);
+ temp *= cheb_eval(dop_n_cheby, dop_cheby_coef, cheb_evmin);
+ temp *= cheb_evmin;
+ temp *= cheb_eval(dop_n_cheby, dop_cheby_coef, cheb_evmin);
+ temp *= chebtilde_eval(ptilde_n_cheby, ptilde_cheby_coef, cheb_evmin);
+ temp = 0.5*fabs(temp - 1);
+ if(g_proc_id == g_stdio_proc){
+ printf(" Delta_IR at s=%f: | Ptilde P s_low P Ptilde - 1 |/2 = %e \n", cheb_evmin, temp);
+ printf("\n Latest (TILDE) polynomial degree = %d\n\n", ptilde_n_cheby);
+ }
+ break;
+ }
+
+ ptilde_n_cheby+=1;
+ }
+
+
+#if ( defined SSE || defined SSE2 || defined SSE3)
+ free(ss_);
+ free(auxs_);
+ free(aux2s_);
+ free(sc_);
+ free(auxc_);
+ free(aux2c_);
+#else
+ free(ss);
+ free(auxs);
+ free(aux2s);
+ free(sc);
+ free(auxc);
+ free(aux2c);
+#endif
+
+}
View
16 Ptilde_nd.h
@@ -0,0 +1,16 @@
+/* $Id$ */
+#ifndef _PTILDE_ND_H
+#define _PTILDE_ND_H
+
+
+double func_tilde(double u, double exponent);
+
+void Ptilde_cheb_coefs(double a, double b, double dd[], int n, double exponent);
+
+void Poly_tilde_ND(spinor *R_s, spinor *R_c, double *dd, int n, spinor *S_s, spinor *S_c);
+
+double chebtilde_eval(int M, double *dd, double s);
+
+void degree_of_Ptilde();
+
+#endif
View
7 README-for-PHMC
@@ -0,0 +1,7 @@
+
+ Untill today, May 17. 2006
+
+ If one want to compute EV, remember to configure with
+
+ --disable-gaugecopy
+
View
97 Square_root_BR_roots.dat
@@ -0,0 +1,97 @@
+ Nr. Re Im
+ 1 -0.891594033208813918633950379444 0.0478046728267078127605493875762
+ 2 -0.0721945631074747268263180899339 0.0644452739791291501214587356117
+ 3 0.860015776358629802089694749156 0.0533783936903994629674663485730
+ 4 1.00495403041123276821622312127 0.00335325450715577479815499017946
+ 5 -0.549587877226238696870552757900 0.0822111314524478636878157544743
+ 6 0.494741269037981323819508361339 0.0844369235580245125483145329781
+ 7 -0.980288435724586904029820288997 0.0232218630408449651369995336836
+ 8 -0.319568757981259221612191367967 0.0869996126360912946928749533981
+ 9 0.700028613987436654220175569208 0.0723786395349235306406754375530
+ 10 0.963958600452450742324117527460 0.0296422177609976349177411947267
+ 11 -0.744753913455955096623029021430 0.0681892943582186744633233388413
+ 12 0.258609123021910025652658760009 0.0859606320884811481164433644153
+ 13 -0.943673956035893768401479064778 0.0359025371040272087141964618695
+ 14 -0.196950493280740906998360628677 0.0833009974950428139228719714993
+ 15 0.786441564269944715093174636422 0.0636038546266462501144189900515
+ 16 0.992595561608037302292473214038 0.0166765814867637640617203942384
+ 17 -0.652455003235876684897220911807 0.0761372932210211544878930567393
+ 18 0.379454821692793431431312001223 0.0869441517611857922487317296145
+ 19 -1.00082882884021739933189110161 0.0100418149048726642463824632046
+ 20 -0.437942079193747113574630702715 0.0860437819939531001178778524263
+ 21 0.602235950014153864628951851046 0.0794284416786994185155634795592
+ 22 0.919518902615733768612926724018 0.0419681119227433077623068413686
+ 23 -0.824915996552996411139702104265 0.0586558423558887903626057891415
+ 24 0.134968186623928010359207974034 0.0776669202577349310212539990061
+ 25 -0.919518902615733768612926724018 0.0419681119227433077623068413686
+ 26 -0.134968186623928010359207974034 0.0776669202577349310212539990061
+ 27 0.824915996552996411139702104265 0.0586558423558887903626057891415
+ 28 1.00082882884021739933189110161 0.0100418149048726642463824632046
+ 29 -0.602235950014153864628951851046 0.0794284416786994185155634795592
+ 30 0.437942079193747113574630702715 0.0860437819939531001178778524263
+ 31 -0.992595561608037302292473214038 0.0166765814867637640617203942384
+ 32 -0.379454821692793431431312001223 0.0869441517611857922487317296145
+ 33 0.652455003235876684897220911807 0.0761372932210211544878930567393
+ 34 0.943673956035893768401479064778 0.0359025371040272087141964618695
+ 35 -0.786441564269944715093174636422 0.0636038546266462501144189900515
+ 36 0.196950493280740906998360628677 0.0833009974950428139228719714993
+ 37 -0.963958600452450742324117527460 0.0296422177609976349177411947267
+ 38 -0.258609123021910025652658760009 0.0859606320884811481164433644153
+ 39 0.744753913455955096623029021430 0.0681892943582186744633233388413
+ 40 0.980288435724586904029820288997 0.0232218630408449651369995336836
+ 41 -0.700028613987436654220175569208 0.0723786395349235306406754375530
+ 42 0.319568757981259221612191367967 0.0869996126360912946928749533981
+ 43 -1.00495403041123276821622312127 0.00335325450715577479815499017946
+ 44 -0.494741269037981323819508361339 0.0844369235580245125483145329781
+ 45 0.549587877226238696870552757900 0.0822111314524478636878157544743
+ 46 0.891594033208813918633950379444 0.0478046728267078127605493875762
+ 47 -0.860015776358629802089694749156 0.0533783936903994629674663485730
+ 48 0.0721945631074747268263180899339 0.0644452739791291501214587356117
+ 49 0.0721945631074747268263180899339 -0.0644452739791291501214587356117
+ 50 -0.860015776358629802089694749156 -0.0533783936903994629674663485730
+ 51 0.891594033208813918633950379444 -0.0478046728267078127605493875762
+ 52 0.549587877226238696870552757900 -0.0822111314524478636878157544743
+ 53 -0.494741269037981323819508361339 -0.0844369235580245125483145329781
+ 54 -1.00495403041123276821622312127 -0.00335325450715577479815499017946
+ 55 0.319568757981259221612191367967 -0.0869996126360912946928749533981
+ 56 -0.700028613987436654220175569208 -0.0723786395349235306406754375530
+ 57 0.980288435724586904029820288997 -0.0232218630408449651369995336836
+ 58 0.744753913455955096623029021430 -0.0681892943582186744633233388413
+ 59 -0.258609123021910025652658760009 -0.0859606320884811481164433644153
+ 60 -0.963958600452450742324117527460 -0.0296422177609976349177411947267
+ 61 0.196950493280740906998360628677 -0.0833009974950428139228719714993
+ 62 -0.786441564269944715093174636422 -0.0636038546266462501144189900515
+ 63 0.943673956035893768401479064778 -0.0359025371040272087141964618695
+ 64 0.652455003235876684897220911807 -0.0761372932210211544878930567393
+ 65 -0.379454821692793431431312001223 -0.0869441517611857922487317296145
+ 66 -0.992595561608037302292473214038 -0.0166765814867637640617203942384
+ 67 0.437942079193747113574630702715 -0.0860437819939531001178778524263
+ 68 -0.602235950014153864628951851046 -0.0794284416786994185155634795592
+ 69 1.00082882884021739933189110161 -0.0100418149048726642463824632046
+ 70 0.824915996552996411139702104265 -0.0586558423558887903626057891415
+ 71 -0.134968186623928010359207974034 -0.0776669202577349310212539990061
+ 72 -0.919518902615733768612926724018 -0.0419681119227433077623068413686
+ 73 0.134968186623928010359207974034 -0.0776669202577349310212539990061
+ 74 -0.824915996552996411139702104265 -0.0586558423558887903626057891415
+ 75 0.919518902615733768612926724018 -0.0419681119227433077623068413686
+ 76 0.602235950014153864628951851046 -0.0794284416786994185155634795592
+ 77 -0.437942079193747113574630702715 -0.0860437819939531001178778524263
+ 78 -1.00082882884021739933189110161 -0.0100418149048726642463824632046
+ 79 0.379454821692793431431312001223 -0.0869441517611857922487317296145
+ 80 -0.652455003235876684897220911807 -0.0761372932210211544878930567393
+ 81 0.992595561608037302292473214038 -0.0166765814867637640617203942384
+ 82 0.786441564269944715093174636422 -0.0636038546266462501144189900515
+ 83 -0.196950493280740906998360628677 -0.0833009974950428139228719714993
+ 84 -0.943673956035893768401479064778 -0.0359025371040272087141964618695
+ 85 0.258609123021910025652658760009 -0.0859606320884811481164433644153
+ 86 -0.744753913455955096623029021430 -0.0681892943582186744633233388413
+ 87 0.963958600452450742324117527460 -0.0296422177609976349177411947267
+ 88 0.700028613987436654220175569208 -0.0723786395349235306406754375530
+ 89 -0.319568757981259221612191367967 -0.0869996126360912946928749533981
+ 90 -0.980288435724586904029820288997 -0.0232218630408449651369995336836
+ 91 0.494741269037981323819508361339 -0.0844369235580245125483145329781
+ 92 -0.549587877226238696870552757900 -0.0822111314524478636878157544743
+ 93 1.00495403041123276821622312127 -0.00335325450715577479815499017946
+ 94 0.860015776358629802089694749156 -0.0533783936903994629674663485730
+ 95 -0.0721945631074747268263180899339 -0.0644452739791291501214587356117
+ 96 -0.891594033208813918633950379444 -0.0478046728267078127605493875762
View
391 chebyshev_polynomial_nd.c
@@ -0,0 +1,391 @@
+/* $Id$ */
+
+#ifdef HAVE_CONFIG_H
+# include<config.h>
+#endif
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include "global.h"
+#include "linsolve.h"
+#include "linalg_eo.h"
+#include "start.h"
+#include "tm_operators.h"
+#include "chebyshev_polynomial_nd.h"
+#include "Nondegenerate_Matrix.h"
+
+
+#define PI 3.141592653589793
+
+double func(double u, double exponent){
+ return pow(u,exponent);
+}
+
+
+void chebyshev_coefs(double aa, double bb, double c[], int n, double exponent){
+ int k,j;
+ double fac,bpa,bma,*f;
+ double inv_n;
+
+
+ inv_n=1./(double)n;
+ f=calloc(n,sizeof(double));/*vector(0,n-1);*/
+ if(g_proc_id == g_stdio_proc){
+ printf("\n hello in chebyshev_polynomial\n");
+ printf("n= %d inv_n=%e \n",n,inv_n);
+ printf("allocation !!!\n");
+ }
+ fflush(stdout);
+ bma=0.5*(bb-aa);
+ bpa=0.5*(bb+aa);
+ for (k=0;k<n;k++) {
+ double y=cos(PI*(k+0.5)*inv_n);
+ f[k]=func(y*bma+bpa,exponent);
+ }
+ fac=2.0*inv_n;
+ for (j=0;j<n;j++) {
+ double sum=0.0;
+ for (k=0;k<n;k++)
+ sum += f[k]*cos(PI*j*(k+0.5)*inv_n);
+ c[j]=fac*sum;
+ }
+ free(f);
+
+
+}
+#undef PI
+
+
+/****************************************************************************
+ *
+ * computation of, despite of the name, (Q Q^dagger) on a vector
+ * by using the chebyshev approximation for the function ()^1/4
+ * subtraction of low-lying eigenvalues is not yet implemented for this
+ *
+ **************************************************************************/
+
+
+void QdaggerQ_poly(spinor *R_s, spinor *R_c, double *c, int n,
+ spinor *S_s, spinor *S_c){
+
+ int j;
+ double fact1, fact2, temp, temp1, temp2, temp3, temp4;
+
+ spinor *svs_=NULL, *svs=NULL, *ds_=NULL, *ds=NULL, *dds_=NULL, *dds=NULL,
+ *auxs_=NULL, *auxs=NULL, *aux2s_=NULL, *aux2s=NULL, *aux3s_=NULL,
+ *aux3s=NULL;
+ spinor *svc_=NULL, *svc=NULL, *dc_=NULL, *dc=NULL, *ddc_=NULL,
+ *ddc=NULL, *auxc_=NULL, *auxc=NULL, *aux2c_=NULL, *aux2c=NULL,
+ *aux3c_=NULL, *aux3c=NULL;
+
+
+#if ( defined SSE || defined SSE2 )
+ svs_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ svs = (spinor *)(((unsigned int)(svs_)+ALIGN_BASE)&~ALIGN_BASE);
+ ds_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ ds = (spinor *)(((unsigned int)(ds_)+ALIGN_BASE)&~ALIGN_BASE);
+ dds_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ dds = (spinor *)(((unsigned int)(dds_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxs_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ auxs = (spinor *)(((unsigned int)(auxs_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2s_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux2s = (spinor *)(((unsigned int)(aux2s_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux3s_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux3s = (spinor *)(((unsigned int)(aux3s_)+ALIGN_BASE)&~ALIGN_BASE);
+ svc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ svc = (spinor *)(((unsigned int)(svc_)+ALIGN_BASE)&~ALIGN_BASE);
+ dc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ dc = (spinor *)(((unsigned int)(dc_)+ALIGN_BASE)&~ALIGN_BASE);
+ ddc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ ddc = (spinor *)(((unsigned int)(ddc_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxc_ = calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ auxc = (spinor *)(((unsigned int)(auxc_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2c_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux2c = (spinor *)(((unsigned int)(aux2c_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux3c_= calloc(VOLUMEPLUSRAND+1, sizeof(spinor));
+ aux3c = (spinor *)(((unsigned int)(aux3c_)+ALIGN_BASE)&~ALIGN_BASE);
+#else
+ svs_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ svs = svs_;
+ ds_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ ds = ds_;
+ dds_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ dds = dds_;
+ auxs_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ auxs = auxs_;
+ aux2s_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux2s = aux2s_;
+ aux3s_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux3s = aux3s_;
+ svc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ svc = svc_;
+ dc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ dc = dc_;
+ ddc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ ddc = ddc_;
+ auxc_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ auxc = auxc_;
+ aux2c_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux2c = aux2c_;
+ aux3c_=calloc(VOLUMEPLUSRAND, sizeof(spinor));
+ aux3c = aux3c_;
+#endif
+
+
+ fact1=4/(cheb_evmax-cheb_evmin);
+ fact2=-2*(cheb_evmax+cheb_evmin)/(cheb_evmax-cheb_evmin);
+
+ zero_spinor_field(&ds[0],VOLUME/2);
+ zero_spinor_field(&dds[0],VOLUME/2);
+ zero_spinor_field(&dc[0],VOLUME/2);
+ zero_spinor_field(&ddc[0],VOLUME/2);
+
+
+ /* sub_low_ev(&aux3[0], &S[0]); */
+ assign(&aux3s[0], &S_s[0],VOLUME/2);
+ assign(&aux3c[0], &S_c[0],VOLUME/2);
+
+ /* Use the Clenshaw's recursion for the Chebysheff polynomial */
+ for (j=n-1; j>=1; j--) {
+ assign(&svs[0],&ds[0],VOLUME/2);
+ assign(&svc[0],&dc[0],VOLUME/2);
+
+ /*
+ if ( (j%10) == 0 ) {
+ sub_low_ev(&aux[0], &d[0]);
+ }
+ else { */
+ assign(&auxs[0], &ds[0], VOLUME/2);
+ assign(&auxc[0], &dc[0], VOLUME/2);
+ /* } */
+
+
+ Q_Qdagger_ND(&R_s[0], &R_c[0], &auxs[0], &auxc[0]);
+
+ temp1=-1.0;
+ temp2=c[j];
+ assign_mul_add_mul_add_mul_add_mul_r(&ds[0] , &R_s[0], &dds[0], &aux3s[0], fact2, fact1, temp1, temp2,VOLUME/2);
+ assign_mul_add_mul_add_mul_add_mul_r(&dc[0] , &R_c[0], &ddc[0], &aux3c[0], fact2, fact1, temp1, temp2,VOLUME/2);
+ assign(&dds[0], &svs[0],VOLUME/2);
+ assign(&ddc[0], &svc[0],VOLUME/2);
+
+ }
+
+ /* sub_low_ev(&R[0],&d[0]); */
+ assign(&R_s[0], &ds[0],VOLUME/2);
+ assign(&R_c[0], &dc[0],VOLUME/2);
+
+
+ Q_Qdagger_ND(&auxs[0], &auxc[0], &R_s[0], &R_c[0]);
+
+ temp1=-1.0;
+ temp2=c[0]/2;
+ temp3=fact1/2;
+ temp4=fact2/2;
+ assign_mul_add_mul_add_mul_add_mul_r(&auxs[0], &ds[0], &dds[0], &aux3s[0], temp3, temp4, temp1, temp2,VOLUME/2);
+ assign_mul_add_mul_add_mul_add_mul_r(&auxc[0], &dc[0], &ddc[0], &aux3c[0], temp3, temp4, temp1, temp2,VOLUME/2);
+ assign(&R_s[0], &auxs[0],VOLUME/2);
+ assign(&R_c[0], &auxc[0],VOLUME/2);
+
+ /* addproj_q_invsqrt(&R[0], &S[0]); */
+
+ /*
+#ifndef _SOLVER_OUTPUT
+ if(g_proc_id == g_stdio_proc){
+ printf("Order of Chebysheff approximation = %d\n",j);
+ fflush( stdout);};
+#endif
+ */
+
+
+ free(svs_);
+ free(ds_);
+ free(dds_);
+ free(auxs_);
+ free(aux2s_);
+ free(aux3s_);
+ free(svc_);
+ free(dc_);
+ free(ddc_);
+ free(auxc_);
+ free(aux2c_);
+ free(aux3c_);
+
+}
+
+
+
+double cheb_eval(int M, double *c, double s){
+
+ double d=0,dd=0, sv, z, z2, res;
+ int j;
+
+ z = (2.0*s - cheb_evmin - cheb_evmax)/(double)(cheb_evmax - cheb_evmin);
+ z2 = 2.0*z;
+
+ for(j=M-1; j>=1; j--){
+ sv = d;
+ d = z2*d - dd + c[j];
+ dd = sv;
+ }
+
+ res = z*d - dd + 0.5*c[0];
+
+ return(res);
+}
+
+/**************************************************************************
+ *
+ * The externally accessible function is
+ *
+ * void degree_of_polynomial_nd(void)
+ * Computation of (QdaggerQ)^1/4
+ * by using the chebyshev approximation for the function ()^1/4
+ *
+ * Author: Mauro Papinutto <papinutt@mail.desy.de> Apr 2003
+ * adapted by Ines Wetzorke <Ines.Wetzorke@desy.de> May 2003
+ * adapted by Karl Jansen <Karl.Jansen@desy.de> June 2005
+ * adapted Thomas Chiarappa <Thomas.Chiarappa@mib.infn.it> Mai 2006
+ *
+*****************************************************************************/
+
+
+void degree_of_polynomial_nd(){
+ int i, j;
+ double temp, temp2;
+ static int ini=0;
+
+ double sum=0.0;
+
+ spinor *ss=NULL, *ss_=NULL, *sc=NULL, *sc_=NULL;
+ spinor *auxs=NULL, *auxs_=NULL, *auxc=NULL, *auxc_=NULL;
+ spinor *aux2s=NULL, *aux2s_=NULL, *aux2c=NULL, *aux2c_=NULL;
+
+
+ if(ini==0){
+ dop_cheby_coef = calloc(N_CHEBYMAX,sizeof(double));
+ ini=1;
+ }
+
+
+#if ( defined SSE || defined SSE2 || defined SSE3)
+ ss_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ auxs_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ aux2s_= calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ sc_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ auxc_ = calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+ aux2c_= calloc(VOLUMEPLUSRAND/2+1, sizeof(spinor));
+
+ ss = (spinor *)(((unsigned int)(ss_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxs = (spinor *)(((unsigned int)(auxs_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2s = (spinor *)(((unsigned int)(aux2s_)+ALIGN_BASE)&~ALIGN_BASE);
+ sc = (spinor *)(((unsigned int)(sc_)+ALIGN_BASE)&~ALIGN_BASE);
+ auxc = (spinor *)(((unsigned int)(auxc_)+ALIGN_BASE)&~ALIGN_BASE);
+ aux2c = (spinor *)(((unsigned int)(aux2c_)+ALIGN_BASE)&~ALIGN_BASE);
+
+#else
+ ss =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ auxs =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ aux2s=calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ sc =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ auxc =calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+ aux2c=calloc(VOLUMEPLUSRAND/2, sizeof(spinor));
+#endif
+
+
+ chebyshev_coefs(cheb_evmin, cheb_evmax, dop_cheby_coef, N_CHEBYMAX, -0.5);
+ /*
+ printf(" \n NchebyMAX = %d \n ", N_CHEBYMAX);
+ for(j=0; j<49; j++){
+ printf(" At %d Coef=%20.18f \n", j, dop_cheby_coef[j]);
+ }
+ */
+
+ random_spinor_field(ss,VOLUME/2);
+ random_spinor_field(sc,VOLUME/2);
+
+ if(g_proc_id == g_stdio_proc){
+ printf(" \n In P: EVmin = %f EVmax = %f \n", cheb_evmin, cheb_evmax);
+ printf("\n determine the degree of the polynomial : Stop=%e \n", g_acc_Pfirst);
+ fflush(stdout);
+ }
+
+ dop_n_cheby=49;
+ for(i = 0;i < 100 ; i++){
+
+ if (dop_n_cheby > N_CHEBYMAX) {
+ if(g_proc_id == g_stdio_proc){
+ printf("Error: n_cheby=%d > N_CHEBYMAX=%d\n",dop_n_cheby,N_CHEBYMAX);
+ printf("Increase n_chebymax\n");
+ }
+ errorhandler(35,"degree_of_polynomial");
+ }
+
+
+ QdaggerQ_poly(&auxs[0], &auxc[0], dop_cheby_coef, dop_n_cheby, &ss[0], &sc[0]);
+
+ Q_Qdagger_ND(&aux2s[0], &aux2c[0], &auxs[0], &auxc[0]);
+
+ QdaggerQ_poly(&auxs[0], &auxc[0], dop_cheby_coef, dop_n_cheby, &aux2s[0], &aux2c[0]);
+
+
+ diff(&aux2s[0],&auxs[0],&ss[0],VOLUME/2);
+ temp=square_norm(&aux2s[0],VOLUME/2)/square_norm(&ss[0],VOLUME/2)/4.0;
+
+ diff(&aux2c[0],&auxc[0],&sc[0],VOLUME/2);
+ temp2=square_norm(&aux2c[0],VOLUME/2)/square_norm(&sc[0],VOLUME/2)/4.0;
+
+ if(g_epsbar == 0){
+ temp2 = 0.0;
+ }
+ if(g_proc_id == g_stdio_proc) {
+ printf("At n=%d || differences ||^2 : UP=%e DN=%e \n",dop_n_cheby, temp, temp2);
+ }
+
+
+ sum=0;
+ for(j=dop_n_cheby; j<N_CHEBYMAX; j++){
+ sum += fabs(dop_cheby_coef[j]);
+ }
+ if(g_proc_id == g_stdio_proc) printf(" Sum remaining | c_n |=%e \n", sum);
+
+ if(sum < g_acc_Pfirst){
+ if(g_proc_id == g_stdio_proc){
+ printf("\n Achieved Accuracies for P : Stop=%e \n", g_acc_Pfirst);
+ printf(" Uniform: Sum |c_n|=%e \n", sum);
+ printf(" RND: || (P S P - 1)X ||^2 /|| 2X ||^2 : UP=%e DN=%e \n",temp, temp2);
+ }
+
+ temp = cheb_eval(dop_n_cheby, dop_cheby_coef, cheb_evmin);
+ temp *= cheb_evmin;
+ temp *= cheb_eval(dop_n_cheby, dop_cheby_coef, cheb_evmin);
+ temp = 0.5*fabs(temp - 1);
+ if(g_proc_id == g_stdio_proc){
+ printf(" Delta_IR at s=%f: | P s_low P - 1 |/2 = %e \n", cheb_evmin, temp);
+ printf("\n Latest (FIRST) polynomial degree = %d \n \n", dop_n_cheby);
+ }
+ break;
+ }
+
+ /* RECALL THAT WE NEED AN EVEN DEGREE !!!! */
+ dop_n_cheby+=2;
+ }
+
+
+#if ( defined SSE || defined SSE2 || defined SSE3)
+ free(ss_);
+ free(auxs_);
+ free(aux2s_);
+ free(sc_);
+ free(auxc_);
+ free(aux2c_);
+#else
+ free(ss);
+ free(auxs);
+ free(aux2s);
+ free(sc);
+ free(auxc);
+ free(aux2c);
+#endif
+
+}
View
16 chebyshev_polynomial_nd.h
@@ -0,0 +1,16 @@
+/* $Id$ */
+#ifndef _CHEBYSHEV_POLYNOMIAL_ND_H
+#define _CHEBYSHEV_POLYNOMIAL_ND_H
+
+
+double func(double u, double exponent);
+
+void chebyshev_coefs(double a, double b, double c[], int n, double exponent);
+
+void QdaggerQ_poly(spinor *R_s, spinor *R_c, double *c, int n, spinor *S_s, spinor *S_c);
+
+double cheb_eval(int M, double *c, double s);
+
+void degree_of_polynomial_nd();
+
+#endif
View
258 clenshaw_coef.c
@@ -0,0 +1,258 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include "clenshaw_coef.h"
+#define Pi 3.141592653589793
+
+extern long double c[3000];
+
+extern double a, b;
+
+long double D[300];
+
+void clenscoef(int M){
+
+ long long int j, jmax, k, N2, N, ij, i, imax;
+
+ long double s[500][500];
+ long double l[500][500];
+ long double sgn;
+ long double sum;
+ long double snew, sold, lnew, lold;
+
+ long double jj;
+
+ long long int j2;
+ long double A, B, A2, B2;
+
+
+ FILE *coeroot;
+ char *filename_stub7 = "Cheby_coeff_for_roots_";
+ char *filename7;
+ char buf7[100];
+
+
+ FILE *factors;
+ char *filename_stub9 = "Pre-factors_";
+ char *filename9;
+ char buf9[100];
+
+
+ filename7=buf7;
+ sprintf(filename7,"%s%d.dat", filename_stub7,M);
+
+ filename9=buf9;
+ sprintf(filename9,"%s%d.dat", filename_stub9,M);
+
+ coeroot = fopen(filename7,"w");
+ fprintf(coeroot,"### Chebishev coeff. in ascending order (pow. coef.) \n");
+ /* fprintf(coeroot,"### power j coeff. \n"); */
+ fclose(coeroot);
+
+
+ N = (long long int)(M - 1);
+ N2 = (long long int)(N/2);
+
+ A = (long double)(2./(long double)(b-a));
+ B = (long double)((b+a)/(long double)(b-a));
+ A2 = (long double)(2*A);
+ B2 = (long double)(2*B);
+
+ /* Initialisation */
+ for(k=0; k<M; k++){
+ for(j=0; j<M; j++){
+ s[k][j] = 0.0;
+ l[k][j] = 0.0;
+ }
+ D[k] = 0.0;
+ }
+
+ factors = fopen(filename9,"w");
+ fprintf(factors," Val. of s[k][j] and l[k][j] \n");
+ fprintf(factors," At k = 1 \n");
+ fclose(factors);
+
+ factors = fopen(filename9,"a");
+ /* Coefficient sequences */
+
+ /* First the k = 1 case */
+ for(j=1; j<M; j++){
+ jj = (long double)(j);
+ s[1][j] = (long double)(2*jj - 1);
+ l[1][j] = (long double)(jj);
+ /* printf(" At j = %d s = %d \n", j, s[1][j]); */
+
+ fprintf(factors," %20.18lle %20.18lle \n", s[1][j], l[1][j]);
+
+ }
+
+ /* then the remaining k cases */
+
+ for(k=2; k<M; k++){
+ fprintf(factors," At k = %d \n", k);
+
+ sold = 0.0;
+ lold = 0.0;
+
+ for(j=1; j<M; j++){
+ snew = (long double)(sold + s[k-1][j]);
+ sold = (long double)(snew);
+ s[k][j] = (long double)(snew);
+
+ lnew = (long double)(lold + l[k-1][j]);
+ lold = (long double)(lnew);
+ l[k][j] = (long double)(lnew);
+
+ fprintf(factors," %20.18lle %20.18lle \n", s[k][j], l[k][j]);
+
+ }
+ }
+
+ fclose(factors);
+
+ for(j=N2; j>=1; j--){
+
+ sgn = -1.0;
+
+ j2=2*j;
+
+ ij = (long long int)((j+1)/2) - (long long int)(j/2);
+ /*
+ printf(" ij=%lld \n", ij);
+ printf(" j=%lld j2=%lld \n", j, j2);
+ */
+ if(ij == 0) sgn = -sgn;
+ /*
+ printf(" sgn=%llf \n", sgn);
+ printf(" C=%llf \n", c[j2]);
+ */
+
+ D[0]+= (long double)(c[j2]*sgn);
+ /*
+ printf(" D=%llf \n", D[0]);
+ */
+ }
+
+
+ D[0] = (long double)(D[0] + 0.5*c[0]);
+ /*
+ printf(" Pre final D=%llf \n", D[0]);
+ */
+
+ /*
+ printf(" D0 = %llf \n", D[0]);
+ */
+
+ /* Evaluate first the coefficient of x^0 */
+ for(i=1; i<M; i++){
+
+ sgn = -1.0;
+
+ ij = (long long int)(i/2) - (long long int)((i-1)/2);
+ if(ij == 0) sgn = -sgn;
+ sum = 0.0;
+ jmax = N2 -(long long int)((i-1)/2);
+ /*
+ printf(" ij=%lld jmax=%lld \n", ij, jmax);
+ */
+ for(j=1; j<=jmax; j++){
+
+ j2 = 2*j + i - 2;
+ sgn = -sgn;
+
+ sum += (long double)(c[j2]*sgn*s[i][j]);
+ /*
+ printf(" j=%lld j2=%lld \n", j, j2);
+ printf(" sgn=%llf \n", sgn);
+ printf(" C=%llf \n", c[j2]);
+ printf(" Sum=%llf \n", sum);
+ */
+ }
+
+
+ sum = (long double)(sum*B);
+
+ if (i > 1) sum = (long double)(sum*powl(B2,(i-1)));
+
+
+ D[0] = (long double)(sum + D[0]);
+ /*
+ printf("At i=%lld Sum=%llf D=%llf D=%20.18lle\n", i, sum, D[0], D[0]);
+ */
+ }
+
+
+
+
+ /* Evaluate the Block of coefficients [1, N-1] */
+
+ for(k=1; k<N; k++){ /* LOOP over degrees */
+ /* for(k=1; k<2; k++){ */
+
+ imax = N - k + 1;
+
+ /* printf(" \n Degree %d Max loop imax=%d \n", k, imax); */
+ /* for i > 1 LOOP over inner loop */
+ for(i=1; i<=imax; i++){
+
+ sgn = 1.0;
+ ij = (long long int)(i/2) - (long long int)((i-1)/2);
+ if(ij == 0) sgn = -sgn;
+ sum = 0.0;
+ jmax = (long long int)((N-k+3-i)/2);
+
+ /* printf(" \n At i=%d ij=%d jmax=%d \n", i, ij, jmax); */
+ for(j=1; j<=jmax; j++){
+
+ j2 = k + 2*j + i - 3;
+ sgn = -sgn;
+ /*
+ printf("At k=%d i=%d jmax=%d j=%d j2=%d \n", k, i, jmax, j, j2);
+ */
+ sum += (long double)(c[j2]*sgn*s[k+i-1][j]);
+ /*
+ printf("s=%d sgn=%llf sum=%llf \n", s[k+i-1][j], sgn, sum);
+ */
+ }
+
+ /* printf(" At k=%d and i=%d Value is %d \n", k, i, l[k][i]); */
+ /* D[k] += sum * l[k][i]; */
+ /* printf(" At degree %d The value is %12.10e \n", k,D[k]); */
+
+ sum = (long double)(sum*l[k][i]*powl(B2,(i-1)));
+
+ D[k] = (long double)((sum + D[k]));
+ /*
+ printf(" At k=%d i=%d, l=%d sum=%llf D=%llf \n", k,i,l[k][i], sum, D[k]);
+ */
+ }
+ D[k] = (long double)(D[k]*powl(A2,k)/2);
+ }
+
+ /* And finally the highest degree coefficient k=N */
+
+ D[N] = (long double)(powl(A2,(N-1))*A*c[N]);
+
+ /* If normalisation is required */
+ /*
+ for(k=0; k<M; k++){
+
+ D[k] = (long double)(D[k]/D[N]);
+
+ }
+ */
+
+
+
+ /* Write all the Clenshaw coefficients in a file */
+ for(k=0; k<M; k++){
+ coeroot = fopen(filename7,"a");
+ fprintf(coeroot," %lld %20.18lle \n", k,D[k]);
+
+ fclose(coeroot);
+ }
+
+
+}
+
+#undef PI
View
6 clenshaw_coef.h
@@ -0,0 +1,6 @@
+#ifndef _CLENSHAW_COEF_H
+#define _CLENSHAW_COEF_H
+
+void clenscoef(int M);
+
+#endif
View
6 default_input_values.h
@@ -21,6 +21,10 @@
#define _default_N_PROC_Y 1
#define _default_N_PROC_Z 1
#define _default_g_kappa 0.125
+#define _default_g_acc_Pfirst 1.e-02
+#define _default_g_acc_Ptilde 1.e-06
+#define _default_g_acc_Hfin 1.e-04
+#define _default_g_rec_ev 1.e+05
#define _default_g_mubar 0.0
#define _default_g_epsbar 0.0
#define _default_g_mu 0.0
@@ -62,7 +66,7 @@
#define _default_max_solver_iterations 5000
#define _default_solver_precision 1.e-15
#define _default_mass_number 0
-#define _default_g_rgi_C1 0.
+#define _default_g_rgi_C1 -0.08333333
#define _default_g_eps_sq_force 1.0e-7
#define _default_g_eps_sq_acc 1.0e-16
#define _default_g_eps_sq_force1 -1.
View
8 eigenvalues_bi.c
@@ -64,7 +64,7 @@ bispinor *eigenvectors = NULL;
double * eigenvls = NULL;
int eigenvalues_for_cg_computed = 0;
-void eigenvalues(int * nr_of_eigenvalues, const int operator_flag,
+void eigenvalues_bi(int * nr_of_eigenvalues, const int operator_flag,
const int max_iterations, const double precision) {
@@ -103,7 +103,6 @@ void eigenvalues(int * nr_of_eigenvalues, const int operator_flag,
/* strcpy(filename,optarg);*/
-
if(g_proc_id == g_stdio_proc) printf("\nNumber of lowest eigenvalues to compute = %d\n\n",(*nr_of_eigenvalues));
eigenvalues_for_cg_computed = 1;
@@ -222,7 +221,10 @@ void eigenvalues(int * nr_of_eigenvalues, const int operator_flag,
(*nr_of_eigenvalues) = converged;
v0dim = converged;
-
+ /*
+ printf(" Lowest EV = %22.15e \n", eigenvls[0]);
+ */
+ cheb_evmin = eigenvls[0];
free(eigenvls);
}
View
2 eigenvalues_bi.h
@@ -5,6 +5,6 @@
extern bispinor * eigenvectors;
extern double * eigenvls;
extern int eigenvalues_for_cg_computed;
-void eigenvalues(int * nev, const int operator_flag, const int max_iterations, const double prec);
+void eigenvalues_bi(int * nev, const int operator_flag, const int max_iterations, const double prec);
#endif
View
44 global.h
@@ -29,7 +29,23 @@
#define DUM_MATRIX (DUM_SOLVER+6)
/* if you want to include bicgstabell */
/* #define DUM_MATRIX (DUM_SOLVER+11) */
-#define NO_OF_SPINORFIELDS (DUM_MATRIX+6)
+
+/* IF PHMC need dum_matrix + 8 instead of + 6 */
+#define NO_OF_SPINORFIELDS (DUM_MATRIX+8)
+/* END IF PHMC */
+
+
+/* IF PHMC:
+ Define here how many bispinors and chi`s one needs ... = Pol degree */
+#define DUM_BI_DERI 6
+#define DUM_BI_SOLVER (DUM_BI_DERI+7)
+#define DUM_BI_MATRIX (DUM_BI_SOLVER+6)
+/* if you want to include bicgstabell */
+/* #define DUM_BI_MATRIX DUM_BI_SOLVER+11 */
+#define NO_OF_BISPINORFIELDS (DUM_BI_MATRIX+6)
+
+/* End IF PHMC*/
+
/* Here you can define antiperiodic */
/* boundary conditions with e.g. */
@@ -43,7 +59,8 @@
#define EPS_SQ3 1.0e-3
#define tiny_t 1.0e-20
-#define N_CHEBYMAX 2000
+#define N_CHEBYMAX 49
+#define NTILDE_CHEBYMAX 2000
#if defined MAIN_PROGRAM
@@ -60,6 +77,19 @@
# define ALIGN
#endif
+/* IF PHMC */
+EXTERN double Cpol;
+EXTERN double cheb_evmin, cheb_evmax;
+EXTERN double invmaxev;
+EXTERN complex * roo;
+EXTERN int dop_n_cheby;
+EXTERN double * dop_cheby_coef;
+EXTERN int ptilde_n_cheby;
+EXTERN double * ptilde_cheby_coef;
+EXTERN double stilde_low, stilde_max;
+/* END PHMC */
+
+
EXTERN double g_eps_sq_force, g_eps_sq_acc;
EXTERN double g_eps_sq_force1, g_eps_sq_force2, g_eps_sq_force3;
EXTERN double g_eps_sq_acc1, g_eps_sq_acc2, g_eps_sq_acc3;
@@ -90,6 +120,13 @@ EXTERN spinor ** g_spinor_field;
EXTERN bispinor ** g_bispinor_field;
+/* IF PHMC */
+EXTERN spinor ** g_chi_up_spinor_field;
+EXTERN spinor ** g_chi_dn_spinor_field;
+EXTERN spinor * g_chi_up_copy;
+EXTERN spinor * g_chi_dn_copy;
+/* End IF PHMC */
+
EXTERN spinor ** g_csg_field[4];
EXTERN int * g_csg_index_array[4];
EXTERN int g_csg_N[8];
@@ -120,6 +157,9 @@ EXTERN double g_kappa, g_c_sw, g_ka_csw_8, g_beta;
EXTERN double g_rgi_C0, g_rgi_C1;
EXTERN double g_mu, g_mu1, g_mu2, g_mu3;
/* Parameters for non-degenrate case */
+EXTERN double g_acc_Pfirst, g_acc_Ptilde;
+EXTERN double g_acc_Hfin;
+EXTERN int g_rec_ev;
EXTERN double g_mubar, g_epsbar;
EXTERN int g_use_clover_flag, g_nr_of_psf;
View
197 global.h-CARST
@@ -0,0 +1,197 @@
+/* $Id$ */
+#ifndef _GLOBAL_H
+#define _GLOBAL_H
+/***************************************************************
+ *
+ * File global.h
+ *
+ * Global parameters and arrays
+ *
+ * Author: Martin Luescher <luscher@mail.desy.de>
+ * Date: 16.03.2001
+ *
+ * Adapted for the HMC-Program by M. Hasenbusch 2002
+ *
+ ***************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#ifdef MPI
+# include <mpi.h>
+#endif
+#ifdef FIXEDVOLUME
+# include "fixed_volume.h"
+#endif
+#include"su3.h"
+#include"su3adj.h"
+
+#define DUM_DERI 6
+#define DUM_SOLVER (DUM_DERI+7)
+#define DUM_MATRIX (DUM_SOLVER+6)
+/* if you want to include bicgstabell */
+/* #define DUM_MATRIX (DUM_SOLVER+11) */
+
+/* IF PHMC need dum_matrix + 8 instead of + 6 */
+#define NO_OF_SPINORFIELDS (DUM_MATRIX+8)
+/* END IF PHMC */
+
+
+/* IF PHMC:
+ Define here how many bispinors and chi`s one needs ... = Pol degree */
+#define DUM_BI_DERI 6
+#define DUM_BI_SOLVER (DUM_BI_DERI+7)
+#define DUM_BI_MATRIX (DUM_BI_SOLVER+6)
+/* if you want to include bicgstabell */
+/* #define DUM_BI_MATRIX DUM_BI_SOLVER+11 */
+#define NO_OF_BISPINORFIELDS (DUM_BI_MATRIX+6)
+
+
+#define NO_OF_CHI_UP_SPINORFIELDS (100+2)
+#define NO_OF_CHI_DN_SPINORFIELDS (100+2)
+/* End IF PHMC*/
+
+
+/* Here you can define antiperiodic */
+/* boundary conditions with e.g. */
+/* #define X1 1. (in x-direction) */
+#define X1 0.
+#define X2 0.
+#define X3 0.
+#define EPS_SQ0 1.0e-16
+#define EPS_SQ1 1.0e-7
+#define EPS_SQ2 1.0e-5
+#define EPS_SQ3 1.0e-3
+#define tiny_t 1.0e-20
+
+#define N_CHEBYMAX 29
+
+
+#if defined MAIN_PROGRAM
+# define EXTERN
+#else
+# define EXTERN extern
+#endif
+
+#if ((defined SSE)||(defined SSE2)||(defined SSE3))
+# include "sse.h"
+#elif defined BGL
+# include "bgl.h"
+#else
+# define ALIGN
+#endif
+
+/* IF PHMC */
+EXTERN double Cpol;
+EXTERN double cheb_evmin, cheb_evmax;
+EXTERN double invmaxev;
+EXTERN complex * roo;
+EXTERN int dop_n_cheby;
+EXTERN double * dop_cheby_coef;
+EXTERN int ptilde_n_cheby;
+EXTERN double * ptilde_cheby_coef;
+/* END PHMC */
+
+
+EXTERN double g_eps_sq_force, g_eps_sq_acc;
+EXTERN double g_eps_sq_force1, g_eps_sq_force2, g_eps_sq_force3;
+EXTERN double g_eps_sq_acc1, g_eps_sq_acc2, g_eps_sq_acc3;
+EXTERN int g_relative_precision_flag;
+EXTERN int g_debug_level;
+
+EXTERN int T_global;
+#ifndef FIXEDVOLUME
+EXTERN int T, L, LX, LY, LZ, VOLUME;
+EXTERN int N_PROC_T, N_PROC_X, N_PROC_Y, N_PROC_Z;
+EXTERN int RAND, EDGES, VOLUMEPLUSRAND;
+#endif
+
+/* translates from lexicographic order to even/odd order */
+EXTERN int * g_lexic2eo;
+/* translates from even/odd orderto lexicograhic order */
+EXTERN int * g_eo2lexic;
+EXTERN int * g_lexic2eosub;
+
+EXTERN int **** g_ipt;
+EXTERN int ** g_iup;
+EXTERN int ** g_idn;
+
+EXTERN int * g_field_z_ipt_even;
+EXTERN int * g_field_z_ipt_odd;
+
+EXTERN spinor ** g_spinor_field;
+
+EXTERN bispinor ** g_bispinor_field;
+
+/* IF PHMC */
+EXTERN spinor ** g_chi_up_spinor_field;
+EXTERN spinor ** g_chi_dn_spinor_field;
+EXTERN spinor * g_chi_up_copy;
+EXTERN spinor * g_chi_dn_copy;
+/* End IF PHMC */
+
+EXTERN spinor ** g_csg_field[4];
+EXTERN int * g_csg_index_array[4];
+EXTERN int g_csg_N[8];
+
+EXTERN su3 ** g_gauge_field;
+EXTERN su3 ** g_gauge_field_back;
+EXTERN su3 ** g_gauge_field_copy;
+/* This is dirty, but dow not allocate memory */
+/* if no clover is used. */
+EXTERN su3adj ** moment;
+EXTERN su3adj ** df0;
+EXTERN su3adj ** ddummy;
+#ifdef CLOVER
+EXTERN su3adj dclover[VOLUMEPLUSRAND][4] ALIGN;
+EXTERN su3 sw[VOLUME][3][2] ALIGN;
+EXTERN su3 sw_inv[VOLUME][3][2] ALIGN;
+EXTERN su3 swp[VOLUME][4] ALIGN;
+EXTERN su3 swm[VOLUME][4] ALIGN;
+#else
+EXTERN su3adj dclover[1][1] ALIGN;
+EXTERN su3 sw[1][1][1] ALIGN;
+EXTERN su3 sw_inv[1][1][1] ALIGN;
+EXTERN su3 swp[1][1] ALIGN;
+EXTERN su3 swm[1][1] ALIGN;
+#endif
+EXTERN int count00,count01,count10,count11,count20,count21;
+EXTERN double g_kappa, g_c_sw, g_ka_csw_8, g_beta;
+EXTERN double g_rgi_C0, g_rgi_C1;
+EXTERN double g_mu, g_mu1, g_mu2, g_mu3;
+/* Parameters for non-degenrate case */
+EXTERN double g_mubar, g_epsbar;
+EXTERN int g_use_clover_flag, g_nr_of_psf;
+
+/* MPI information */
+EXTERN int g_proc_id, g_nproc, g_stdio_proc, g_nproc_t, g_nproc_x, g_nproc_y, g_nproc_z, g_cart_id;
+EXTERN int g_proc_coords[4];
+EXTERN int g_dbw2rand;
+#ifdef MPI
+EXTERN MPI_Status status;
+EXTERN MPI_Request req1,req2,req3,req4;
+EXTERN MPI_Comm g_cart_grid;
+
+/* the next neighbours for MPI */
+EXTERN int g_nb_x_up, g_nb_x_dn;
+EXTERN int g_nb_y_up, g_nb_y_dn;
+EXTERN int g_nb_t_up, g_nb_t_dn;
+EXTERN int g_nb_z_up, g_nb_z_dn;
+
+#endif
+
+/* something to evaluate time elaps */
+EXTERN double DeltaTtot, DeltaTcd, DeltaTev;
+EXTERN int counter_Spsi;
+/* end of the something ... */
+
+EXTERN int ITER_MAX_BCG;
+EXTERN int ITER_MAX_CG;
+
+#undef EXTERN
+/* #undef ALIGN */
+
+#ifdef MAIN_PROGRAM
+static char const g_rcsid[] = "$Id$";
+#endif
+
+#endif
+
View
197 global.h-TOM
@@ -0,0 +1,197 @@
+/* $Id$ */
+#ifndef _GLOBAL_H
+#define _GLOBAL_H
+/***************************************************************
+ *
+ * File global.h
+ *
+ * Global parameters and arrays
+ *
+ * Author: Martin Luescher <luscher@mail.desy.de>
+ * Date: 16.03.2001
+ *
+ * Adapted for the HMC-Program by M. Hasenbusch 2002
+ *
+ ***************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#ifdef MPI
+# include <mpi.h>
+#endif
+#ifdef FIXEDVOLUME
+# include "fixed_volume.h"
+#endif
+#include"su3.h"
+#include"su3adj.h"
+
+#define DUM_DERI 6
+#define DUM_SOLVER (DUM_DERI+7)
+#define DUM_MATRIX (DUM_SOLVER+6)
+/* if you want to include bicgstabell */
+/* #define DUM_MATRIX (DUM_SOLVER+11) */
+
+/* IF PHMC need dum_matrix + 8 instead of + 6 */
+#define NO_OF_SPINORFIELDS (DUM_MATRIX+8)
+/* END IF PHMC */
+
+
+/* IF PHMC:
+ Define here how many bispinors and chi`s one needs ... = Pol degree */
+#define DUM_BI_DERI 6
+#define DUM_BI_SOLVER (DUM_BI_DERI+7)
+#define DUM_BI_MATRIX (DUM_BI_SOLVER+6)
+/* if you want to include bicgstabell */
+/* #define DUM_BI_MATRIX DUM_BI_SOLVER+11 */
+#define NO_OF_BISPINORFIELDS (DUM_BI_MATRIX+6)
+
+
+#define NO_OF_CHI_UP_SPINORFIELDS (100+2)
+#define NO_OF_CHI_DN_SPINORFIELDS (100+2)
+/* End IF PHMC*/
+
+
+/* Here you can define antiperiodic */
+/* boundary conditions with e.g. */
+/* #define X1 1. (in x-direction) */
+#define X1 0.
+#define X2 0.
+#define X3 0.
+#define EPS_SQ0 1.0e-16
+#define EPS_SQ1 1.0e-7
+#define EPS_SQ2 1.0e-5
+#define EPS_SQ3 1.0e-3
+#define tiny_t 1.0e-20
+
+#define N_CHEBYMAX 2000
+
+
+#if defined MAIN_PROGRAM
+# define EXTERN
+#else
+# define EXTERN extern
+#endif
+
+#if ((defined SSE)||(defined SSE2)||(defined SSE3))
+# include "sse.h"
+#elif defined BGL
+# include "bgl.h"
+#else
+# define ALIGN
+#endif
+
+/* IF PHMC */
+EXTERN double Cpol;
+EXTERN double cheb_evmin, cheb_evmax;
+EXTERN double invmaxev;
+EXTERN complex * roo;
+EXTERN int dop_n_cheby;
+EXTERN double * dop_cheby_coef;
+EXTERN int ptilde_n_cheby;
+EXTERN double * ptilde_cheby_coef;
+/* END PHMC */
+
+
+EXTERN double g_eps_sq_force, g_eps_sq_acc;
+EXTERN double g_eps_sq_force1, g_eps_sq_force2, g_eps_sq_force3;
+EXTERN double g_eps_sq_acc1, g_eps_sq_acc2, g_eps_sq_acc3;
+EXTERN int g_relative_precision_flag;
+EXTERN int g_debug_level;
+
+EXTERN int T_global;
+#ifndef FIXEDVOLUME
+EXTERN int T, L, LX, LY, LZ, VOLUME;
+EXTERN int N_PROC_T, N_PROC_X, N_PROC_Y, N_PROC_Z;
+EXTERN int RAND, EDGES, VOLUMEPLUSRAND;
+#endif
+
+/* translates from lexicographic order to even/odd order */
+EXTERN int * g_lexic2eo;
+/* translates from even/odd orderto lexicograhic order */
+EXTERN int * g_eo2lexic;
+EXTERN int * g_lexic2eosub;
+
+EXTERN int **** g_ipt;
+EXTERN int ** g_iup;
+EXTERN int ** g_idn;
+
+EXTERN int * g_field_z_ipt_even;
+EXTERN int * g_field_z_ipt_odd;
+
+EXTERN spinor ** g_spinor_field;
+
+EXTERN bispinor ** g_bispinor_field;
+
+/* IF PHMC */
+EXTERN spinor ** g_chi_up_spinor_field;
+EXTERN spinor ** g_chi_dn_spinor_field;
+EXTERN spinor * g_chi_up_copy;
+EXTERN spinor * g_chi_dn_copy;
+/* End IF PHMC */
+