5
5
#define CONTROL
6
6
#include "ks_eig_includes.h" /* definitions files and prototypes */
7
7
8
- EXTERN gauge_header start_lat_hdr ; /* Input gauge field header */
8
+ #ifdef HAVE_QUDA
9
+ #include <quda_milc_interface.h>
10
+ #endif
11
+ #ifdef U1_FIELD
12
+ #include "../include/io_u1lat.h"
13
+ #endif
9
14
10
15
int main ( int argc , char * * argv ){
11
16
register site * s ;
12
17
int i ,si ;
13
18
int prompt ;
14
19
double dtime ;
15
20
// su3_vector **eigVec ;
16
- su3_vector * tmp ;
17
21
// double *eigVal ;
18
22
int total_R_iters ;
19
23
double * resid = NULL ;
@@ -26,7 +30,9 @@ int main( int argc, char **argv ){
26
30
27
31
g_sync ();
28
32
/* set up */
33
+ STARTTIME ;
29
34
prompt = setup ();
35
+ ENDTIME ("setup" );
30
36
31
37
/* loop over input sets */
32
38
while ( readin (prompt ) == 0 ){
@@ -38,36 +44,74 @@ int main( int argc, char **argv ){
38
44
#ifdef HYPISQ_SVD_COUNTER
39
45
hypisq_svd_counter = 0 ;
40
46
#endif
41
- restore_fermion_links_from_site (fn_links , MILC_PRECISION );
42
47
48
+ /* Fix the gauge */
49
+
50
+ if ( param .fixflag == COULOMB_GAUGE_FIX )
51
+ {
52
+ if (this_node == 0 )
53
+ printf ("Fixing to Coulomb gauge\n" );
54
+
55
+ rephase ( OFF );
56
+ STARTTIME ;
57
+ gaugefix (TUP ,(Real )1.8 ,500 ,GAUGE_FIX_TOL );
58
+ ENDTIME ("gauge fix" );
59
+ }
60
+ else
61
+ if (this_node == 0 )printf ("COULOMB GAUGE FIXING SKIPPED.\n" );
62
+
63
+ /* save lattice if requested */
64
+ if ( param .saveflag != FORGET ){
65
+ rephase ( OFF );
66
+ savelat_p = save_lattice ( param .saveflag , param .savefile ,
67
+ param .stringLFN );
68
+ rephase ( ON );
69
+ }
70
+
71
+ #ifdef U1_FIELD
72
+ if ( param .save_u1flag != FORGET ){
73
+ save_u1_lattice ( param .save_u1flag , param .save_u1file );
74
+ }
75
+
76
+ /* Insert U(1) phases into the gauge links in the site structure */
77
+ u1phase_on (param .charge , u1_A );
78
+ invalidate_fermion_links (fn_links );
79
+ #endif
80
+
81
+ /* Recompute FN links from the links in the site structure */
82
+ restore_fermion_links_from_site (fn_links , MILC_PRECISION );
83
+
84
+ /* Eigenpair calculation */
85
+ STARTTIME ;
86
+
43
87
eigVal = (double * )malloc (param .eigen_param .Nvecs * sizeof (double ));
44
88
eigVec = (su3_vector * * )malloc (param .eigen_param .Nvecs * sizeof (su3_vector * ));
45
89
for (i = 0 ;i < param .eigen_param .Nvecs ;i ++ )
46
- eigVec [i ]=
47
- (su3_vector * )malloc (sites_on_node * sizeof (su3_vector ));
48
-
49
- /* call fermion_variable measuring routines */
50
- /* results are printed in output file */
51
- f_meas_imp ( 1 , MILC_PRECISION , F_OFFSET (phi ), F_OFFSET (xxx ), mass ,
52
- 0 , fn_links );
53
-
90
+ eigVec [i ]= (su3_vector * )malloc (sites_on_node * sizeof (su3_vector ));
91
+
54
92
fn = get_fm_links (fn_links );
93
+ /* Calculate eigenpairs on even sites */
55
94
param .eigen_param .parity = EVEN ;
56
95
total_R_iters = ks_eigensolve (eigVec , eigVal , & param .eigen_param , 1 );
57
- construct_eigen_odd ( eigVec , eigVal , & param . eigen_param , fn [ 0 ] );
96
+ fflush ( stdout );
58
97
98
+ /* Construct eigenpairs on odd sites */
99
+ construct_eigen_odd (eigVec , eigVal , & param .eigen_param , fn [0 ]);
100
+
59
101
/* Calculate and print the residues and norms of the eigenvectors */
60
102
resid = (double * )malloc (param .eigen_param .Nvecs * sizeof (double ));
61
103
node0_printf ("Even site residuals\n" );
62
104
check_eigres ( resid , eigVec , eigVal , param .eigen_param .Nvecs , EVEN , fn [0 ] );
63
105
node0_printf ("Odd site residuals\n" );
64
106
check_eigres ( resid , eigVec , eigVal , param .eigen_param .Nvecs , ODD , fn [0 ] );
65
-
107
+ fflush (stdout );
108
+
66
109
/* save eigenvectors if requested */
67
110
int status = save_ks_eigen (param .ks_eigen_saveflag , param .ks_eigen_savefile ,
68
- param .eigen_param .Nvecs , eigVal , eigVec , resid , 1 );
111
+ param .eigen_param .Nvecs , eigVal , eigVec , resid , 1 );
112
+
69
113
if (status != 0 ){
70
- node0_printf ("ERROR writing eigenvectors\n" );
114
+ node0_printf ("ERROR writing eigenvectors\n" ); fflush ( stdout );
71
115
}
72
116
73
117
/* print eigenvalues of iDslash */
@@ -80,92 +124,64 @@ int main( int argc, char **argv ){
80
124
node0_printf ("eigenval(%i): %10g\n" ,i ,chirality ) ;
81
125
}
82
126
83
- tmp = (su3_vector * )malloc (sites_on_node * sizeof (su3_vector ));
84
- for (i = 0 ;i < param .eigen_param .Nvecs ;i ++ )
85
- {
86
- // if ( eigVal[i] > 10.0*eigenval_tol )
87
- // {
88
- /* Construct to odd part of the vector. *
89
- * Note that the true odd part of the eigenvector is *
90
- * i/sqrt(eigVal) Dslash Psi. But since I only compute *
91
- * the chirality the i factor is irrelevant (-i)*i=1!! */
92
- dslash_field (eigVec [i ], tmp , ODD , fn [0 ]) ;
93
- FORSOMEPARITY (si ,s ,ODD ){
94
- scalar_mult_su3_vector ( & (tmp [si ]),
95
- 1.0 /sqrt (eigVal [i ]),
96
- & (eigVec [i ][si ]) ) ;
97
- }
98
-
99
- // measure_chirality(eigVec[i], &chirality, EVENANDODD);
100
- /* Here I divide by 2 since the EVEN vector is normalized to
101
- * 1. The EVENANDODD vector is normalized to 2. I could have
102
- * normalized the EVENANDODD vector to 1 and then not devide
103
- * by to. The measure_chirality routine assumes vectors
104
- * normalized to 1. */
105
- // node0_printf("Chirality(%i): %g\n",i,chirality/2) ;
106
- measure_chirality (eigVec [i ], & chir_ev , EVEN );
107
- measure_chirality (eigVec [i ], & chir_od , ODD );
108
- chirality = (chir_ev + chir_od ) / 2.0 ;
109
- node0_printf ("Chirality(%i) -- even, odd, total: %10g, %10g, %10g\n" ,
110
- i ,chir_ev ,chir_od ,chirality ) ;
111
- // }
112
- // else
113
- // {
114
- /* This is considered a "zero mode", and treated as such */
115
- // measure_chirality(eigVec[i], &chirality, EVEN);
116
- // node0_printf("Chirality(%i): %g\n",i,chirality) ;
117
- // }
118
- }
119
- #ifdef EO
120
- cleanup_dslash_temps ();
121
- #endif
122
- free (tmp );
123
- /**
124
- for(i=0;i<Nvecs;i++)
125
- {
126
- sprintf(label,"DENSITY(%i)",i) ;
127
- print_densities(eigVec[i], label, ny/2,nz/2,nt/2, EVEN) ;
128
- }
129
- **/
130
- for (i = 0 ;i < param .eigen_param .Nvecs ;i ++ )
131
- free (eigVec [i ]) ;
132
- free (eigVec ) ;
133
- free (eigVal ) ;
127
+ #ifdef U1_FIELD
128
+ /* Unapply the U(1) field phases */
129
+ u1phase_off ();
134
130
invalidate_fermion_links (fn_links );
135
- fflush (stdout );
131
+ #endif
132
+ ENDTIME ("calculate Dirac eigenpairs" );
136
133
137
- node0_printf ("RUNNING COMPLETED\n" ); fflush (stdout );
138
- dtime += dclock ();
139
- if (this_node == 0 ){
140
- printf ("Time = %e seconds\n" ,dtime );
141
- printf ("total_iters = %d\n" ,total_iters );
142
- printf ("total Rayleigh iters = %d\n" ,total_R_iters );
134
+ #ifdef CHIRALITY
135
+ measure_chirality (eigVec [i ], & chir_ev , EVEN );
136
+ measure_chirality (eigVec [i ], & chir_od , ODD );
137
+ chirality = (chir_ev + chir_od ) / 2.0 ;
138
+ node0_printf ("Chirality(%i) -- even, odd, total: %10g, %10g, %10g\n" ,
139
+ i ,chir_ev ,chir_od ,chirality ) ;
140
+ #endif
141
+ }
142
+ #ifdef EO
143
+ cleanup_dslash_temps ();
144
+ #endif
145
+ /**
146
+ for(i=0;i<Nvecs;i++)
147
+ {
148
+ sprintf(label,"DENSITY(%i)",i) ;
149
+ print_densities(eigVec[i], label, ny/2,nz/2,nt/2, EVEN) ;
150
+ }
151
+ **/
152
+
153
+ /* Clean up eigen storage */
154
+ for (i = 0 ; i < param .eigen_param .Nvecs ; i ++ ) free (eigVec [i ]);
155
+ free (eigVal ); free (eigVec ); free (resid );
156
+ invalidate_fermion_links (fn_links );
157
+
158
+ node0_printf ("RUNNING COMPLETED\n" ); fflush (stdout );
159
+ dtime += dclock ();
160
+ if (this_node == 0 ){
161
+ printf ("Time = %e seconds\n" ,dtime );
162
+ printf ("total Rayleigh iters = %d\n" ,total_R_iters );
143
163
#ifdef HISQ_SVD_COUNTER
144
- printf ("hisq_svd_counter = %d\n" ,hisq_svd_counter );
164
+ printf ("hisq_svd_counter = %d\n" ,hisq_svd_counter );
145
165
#endif
146
166
#ifdef HYPISQ_SVD_COUNTER
147
- printf ("hypisq_svd_counter = %d\n" ,hypisq_svd_counter );
167
+ printf ("hypisq_svd_counter = %d\n" ,hypisq_svd_counter );
148
168
#endif
149
- }
150
- fflush (stdout );
151
-
152
- destroy_ape_links_4D (ape_links );
153
-
154
- /* Destroy fermion links (created in readin() */
155
-
169
+ }
170
+ fflush (stdout );
171
+
172
+ /* Destroy fermion links (created in readin() */
173
+
156
174
#if FERM_ACTION == HISQ
157
- destroy_fermion_links_hisq (fn_links );
175
+ destroy_fermion_links_hisq (fn_links );
158
176
#else
159
- destroy_fermion_links (fn_links );
177
+ destroy_fermion_links (fn_links );
160
178
#endif
161
- fn_links = NULL ;
162
-
163
- } /* readin(prompt) */
164
-
179
+ fn_links = NULL ;
180
+
165
181
#ifdef HAVE_QUDA
166
182
qudaFinalize ();
167
183
#endif
168
-
184
+
169
185
normal_exit (0 );
170
186
return 0 ;
171
187
}
0 commit comments