CHROMA
polar_dec.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Decompose a complex matrix as C = exp(i\alpha) V P
3  */
4 
5 #include "chromabase.h"
6 #include "meas/gfix/polar_dec.h"
7 #include "util/gauge/reunit.h"
8 
9 namespace Chroma {
10 //! Decompose a complex matrix as C = exp(i\alpha) V P
11 /*!
12  * \ingroup gfix
13  *
14  * Decompose a complex matrix as C = exp(i\alpha) V P
15  * with V SU(Nc) and P = (C^\dagger C)^{1/2} positive hermitian
16  *
17  * \param c complex Nc x Nc matrix ( Modify )
18  * on exit it contains the hermitian matrix P
19  * \param v the projected SU(Nc) Matrix ( Write )
20  * \param alpha the phase ( Write )
21  * \param JacAccu accuracy in the Jacobi iteration ( Read )
22  * \param JacMax maximum number of Jacobi iterations ( Read )
23  */
24 
25 void polar_dec(LatticeColorMatrix& c, LatticeColorMatrix& v,
26  LatticeReal& alpha, const Real& JacAccu, int JacMax)
27 {
28  LatticeColorMatrix u_tmp;
29  LatticeColorMatrix w_tmp;
30  LatticeColorMatrix v_tmp;
31 
32  multi2d<LatticeComplex> mat1(Nc, Nc);
33  multi2d<LatticeComplex> mat2(Nc, Nc);
34 
35  LatticeComplex lc1;
36  LatticeComplex lc2;
37  LatticeComplex v12;
38  LatticeComplex v21;
39  multi1d<LatticeReal> diag(Nc);
40  LatticeReal rescale;
41  LatticeReal hundred;
42  LatticeReal v11;
43  LatticeReal dd;
44  LatticeReal diff_diag;
45  LatticeReal tt;
46  LatticeReal cc;
47  LatticeReal ss;
48  LatticeReal al1;
49  LatticeReal al2;
50  LatticeReal lr1;
51  LatticeReal lr2;
52  LatticeBoolean l_rot;
53  LatticeBoolean lb1;
54  LatticeBoolean lb2;
55  LatticeBoolean lb3;
56 
57  Double off_d;
58  Double diff_sq;
59  Double det_diff;
60  Double unit_test;
61  int iter;
62  int i_rot;
63  int n_rot;
64  int numbad;
65 
66  START_CODE();
67 
68  Real acc_sq = JacAccu * JacAccu / Real(100);
69 
70  u_tmp = adj(c) * c;
71 
72  /* Bring diagonal elements to be O(1) */
73  rescale = real(trace(u_tmp));
74  lr1 = 1;
75  lr1 = lr1 / rescale;
76  u_tmp = u_tmp * lr1;
77 
78  // Extract initial components
79  for(int i=0; i < Nc; ++i)
80  for(int j=0; j < Nc; ++j)
81  mat1[i][j] = peekColor(u_tmp, i, j);
82 
83  for(int i=0; i<Nc; i++)
84  {
85  for(int j=0; j<i; j++)
86  {
87  mat2[i][j] = 0;
88  mat2[j][i] = 0;
89  }
90  mat2[i][i] = 1;
91  diag[i] = real(mat1[i][i]);
92  }
93 
94  /* Diagonalize mat1 by Jacobi iteration and keep "rotation" in mat2 */
95 
96  hundred = 100;
97  iter = 0;
98  i_rot = Layout::vol();
99 
100  while ( iter < JacMax && i_rot > 0 )
101  {
102  iter++;
103  i_rot = 0;
104 
105  for(int j = 1; j < Nc; j++)
106  for(int i = 0; i < j; i++)
107  {
108  dd = real(mat1[i][j] * mat1[j][i]);
109  lr1 = fabs(diag[i] * diag[j] * acc_sq);
110 
111  l_rot = dd > lr1;
112  n_rot = toInt(sum(where(l_rot, LatticeInteger(1), LatticeInteger(0))));
113 
114  if( n_rot > 0 )
115  {
116  i_rot += n_rot;
117  dd = sqrt(dd);
118  lr1 = dd * hundred;
119  diff_diag = diag[j] - diag[i];
120 
121  lr2 = fabs(diff_diag);
122  lr1 += lr2;
123  lb1 = lr2 == lr1;
124  lb2 = l_rot & lb1;
125  lr1 = dd / diff_diag;
126  copymask(tt, lb2, lr1);
127 
128  lb1 = !lb1;
129  lb2 = l_rot & lb1;
130  lr1 = diff_diag / dd;
131  LatticeReal theta = 0.5 * lr1;
132  lr1 = fabs(theta) + sqrt(1 + theta * theta);
133  lr2 = 1 / lr1;
134  copymask(tt, lb2, lr2);
135 
136  lr1 = 0;
137  lb1 = theta < lr1;
138  lb3 = lb2 & lb1;
139  copymask(tt, lb3, LatticeReal(-tt));
140 
141  lr1 = 1 / sqrt(1 + tt*tt);
142  lr2 = tt * lr1;
143  copymask(cc, l_rot, lr1);
144  copymask(ss, l_rot, lr2);
145 
146  lr1 = cc * cc;
147  al1 = lr1 * diag[i];
148  al2 = lr1 * diag[j];
149  lr1 = ss * ss;
150  al1 += lr1 * diag[j];
151  al2 += lr1 * diag[i];
152  lr1 = cc * ss * dd;
153 
154  lr1 *= 2;
155  al1 -= lr1;
156  al2 += lr1;
157  copymask(diag[i], l_rot, al1);
158  copymask(diag[j], l_rot, al2);
159 
160  v11 = cc;
161  lr1 = ss / dd;
162  v12 = mat1[j][i] * lr1;
163  v21 = -adj(v12);
164  lc1 = 0;
165  copymask(mat1[j][i], l_rot, lc1);
166  copymask(mat1[i][j], l_rot, lc1);
167 
168  /* Rotate the remaining matrix elements */
169  for(int k = 0; k < Nc; k++)
170  if( k != i && k != j )
171  {
172  lc1 = mat1[k][i] * v11 - mat1[k][j] * v12;
173  lc2 = mat1[k][j] * v11 - mat1[k][i] * v21;
174  copymask(mat1[k][i], l_rot, lc1);
175  copymask(mat1[k][j], l_rot, lc2);
176  copymask(mat1[i][k], l_rot, LatticeComplex(adj(lc1)));
177  copymask(mat1[j][k], l_rot, LatticeComplex(adj(lc2)));
178  }
179 
180  /* Rotate the eigenstd::vector matrix mat2 */
181  for(int k = 0; k < Nc; k++)
182  {
183  lc1 = mat2[i][k] * v11 + mat2[j][k] * v21;
184  lc2 = mat2[j][k] * v11 + mat2[i][k] * v12;
185  copymask(mat2[i][k], l_rot, lc1);
186  copymask(mat2[j][k], l_rot, lc2);
187  }
188  }
189  }
190  }
191 
192 
193  for(int i=0; i < Nc; ++i)
194  for(int j=0; j < Nc; ++j)
195  pokeColor(w_tmp, mat2[i][j], i, j);
196 
197  v_tmp = w_tmp;
198 
199  /* Check unitarity */
200  reunit(v_tmp, lb1, numbad, REUNITARIZE_LABEL);
201  if (numbad > 0)
202  printf("Warning: %d rotation matrices not unitary\n", numbad);
203 /* */
204 
205  /* Check diagonalization */
206  lr1 = 0;
207  off_d = 0;
208  numbad = 0;
209  for(int i=0; i<Nc; i++)
210  {
211  for(int j=0; j<i; j++)
212  {
213  off_d += norm2(mat1[j][i]);
214  mat1[j][i] = 0;
215  mat1[i][j] = 0;
216  }
217  numbad += toInt(sum(where(diag[i] <= lr1, LatticeInteger(1), LatticeInteger(0))));
218  mat1[i][i] = cmplx(diag[i],lr1);
219  }
220 
221  for(int i=0; i < Nc; ++i)
222  for(int j=0; j < Nc; ++j)
223  pokeColor(v, mat1[i][j], i, j);
224 
225  v_tmp = v * w_tmp ;
226  v = adj(w_tmp) * v_tmp - u_tmp;
227  diff_sq = norm2(v) / Double(Layout::vol());
228  off_d /= Double(Layout::vol());
229  if (numbad > 0)
230  {
231  QDPIO::cout << "Warning: " << numbad << " C matrices have zero determiant" << std::endl;
232 #if 0
233  push(nml_out,"Bad_C_matrices");
234  write(nml_out, "numbad", numbad);
235  pop(nml_out);
236 #endif
237  }
238 
239  /* Rescale eigenvalues and construct P^{-1} in v and c*P^{-1} in u_tmp */
240  for(int i=0; i<Nc; i++)
241  {
242  diag[i] = sqrt(diag[i] * rescale);
243  lr2 = 1 / diag[i];
244  mat1[i][i] = cmplx(lr2,lr1);
245  }
246 
247  for(int i=0; i < Nc; ++i)
248  for(int j=0; j < Nc; ++j)
249  pokeColor(v, mat1[i][j], i, j);
250 
251  u_tmp = v * w_tmp;
252  v = adj(w_tmp) * u_tmp;
253  u_tmp = c * v;
254 
255  /* Constuct P in c */
256  for(int i=0; i<Nc; i++)
257  {
258  mat1[i][i] = cmplx(diag[i],lr1);
259  }
260 
261  for(int i=0; i < Nc; ++i)
262  for(int j=0; j < Nc; ++j)
263  pokeColor(c, mat1[i][j], i, j);
264 
265  v_tmp = c * w_tmp;
266  c = adj(w_tmp) * v_tmp;
267 
268  /* Check that that u_tmp^dagger * u_tmp = 1 */
269  v_tmp = 1;
270  v_tmp -= adj(u_tmp) * u_tmp;
271  unit_test = norm2(v_tmp);
272  unit_test /= Double(Layout::vol());
273 
274 
275  /* Now we just need to remove the phase alpha from u_tmp to make it v */
276  /* We compute the determinant by LU decomposition */
277  for(int i=0; i < Nc; ++i)
278  for(int j=0; j < Nc; ++j)
279  mat1[i][j] = peekColor(u_tmp, i, j);
280 
281  lr2 = 1;
282  for(int j = 0; j < Nc; j++)
283  {
284  for(int i = 0; i <= j; i++)
285  {
286  lc1 = mat1[j][i];
287  for(int k = 0; k < i; k++)
288  lc1 -= mat1[k][i] * mat1[j][k];
289 
290  mat1[j][i] = lc1;
291  }
292 
293  for(int i = (j+1); i < Nc; i++)
294  {
295  lc1 = mat1[j][i];
296  for(int k = 0; k < j; k++)
297  lc1 -= mat1[k][i] * mat1[j][k];
298 
299  lc2 = adj(mat1[j][j]) * mat1[j][j];
300  lr1 = real(lc2);
301  lr1 = lr2 / lr1;
302  lc2 = lc1 * lr1;
303  lc1 = adj(mat1[j][j]) * lc2;
304  mat1[j][i] = lc1;
305  }
306  }
307 
308  /* The determinant */
309  lc1 = mat1[0][0] * mat1[1][1];
310  for(int k = 2; k < Nc; k++)
311  {
312  lc2 = mat1[k][k] * lc1;
313  lc1 = lc2;
314  }
315  lr1 = real(lc1);
316  lr2 = imag(lc1);
317  alpha = atan2(lr2, lr1);
318 
319  lr1 = 1;
320  lr2 = real(trace(adj(lc1) * lc1));
321  lr1 -= lr2;
322  det_diff = norm2(lr1) / Double(Layout::vol());
323 
324  /* Make u_tmp unitary in v */
325  lr1 = Real(Nc);
326  alpha = alpha / lr1;
327  lr2 = log(lr2);
328  lr2 = lr2 / lr1;
329  lr2 = -lr2;
330  lr2 = exp(lr2);
331  cc = cos(alpha);
332  ss = sin(alpha);
333  lr1 = cc * lr2;
334  lr2 = -(ss * lr2);
335  lc1 = cmplx(lr1,lr2);
336  v = u_tmp * lc1;
337 
338  /* Make v unitary */
339  reunit(v, lb1, numbad, REUNITARIZE_LABEL);
340 
341 #if 0
342  XMLBufferWriter xml_out;
343  push(xml_out,"Diagonalization_test");
344  write(xml_out, "iter",iter);
345  write(xml_out, "off_d",off_d);
346  write(xml_out, "diff_sq", diff_sq);
347  write(xml_out, "unit_test", unit_test);
348  write(xml_out, "det_diff", det_diff);
349  write(xml_out, "numbad", numbad);
350  pop(xml_out);
351 #endif
352 
353 #if 0
354  // DEBUG DEBUG DEBUG
355  QDPIO::cout << "polar_dec::iter " << iter << "\n" ;
356  QDPIO::cout << "polar_dec::off_d " << off_d << "\n" ;
357  QDPIO::cout << "polar_dec::diff_sq " << diff_sq << "\n" ;
358  QDPIO::cout << "polar_dec::unit_test " << unit_test << "\n" ;
359  QDPIO::cout << "polar_dec::det_diff " << det_diff<< "\n" ;
360  QDPIO::cout << "polar_dec::numbad " << numbad<< "\n" ;
361 #endif
362 
363  END_CODE();
364 }
365 
366 } // End namespace Chroma
chromabase.h
Primary include file for CHROMA library code.
numbad
int numbad
Definition: cool.cc:28
Chroma::write
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
Definition: asqtad_fermact_params_s.cc:40
Chroma::polar_dec
void polar_dec(LatticeColorMatrix &c, LatticeColorMatrix &v, LatticeReal &alpha, const Real &JacAccu, int JacMax)
Decompose a complex matrix as C = exp(i\alpha) V P.
Definition: polar_dec.cc:25
j
unsigned j
Definition: ldumul_w.cc:35
Chroma::StagPhases::alpha
static const LatticeInteger & alpha(const int dim)
Definition: stag_phases_s.h:43
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::push
push(xml_out,"Condensates")
Chroma::c
Double c
Definition: invbicg.cc:108
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::reunit
void reunit(LatticeColorMatrixF3 &xa)
Definition: reunit.cc:467
Chroma::REUNITARIZE_LABEL
@ REUNITARIZE_LABEL
Definition: reunit.h:29
Chroma::pop
pop(xml_out)
Chroma::START_CODE
START_CODE()
Chroma::k
int k
Definition: invbicg.cc:119
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
polar_dec.h
Decompose a complex matrix as .
copymask
copymask(lcoord, lbit, ltmp_1, REPLACE)
sum
Double sum
Definition: qtopcor.cc:37
reunit.h
Reunitarize in place a color matrix to SU(N)
Generated by doxygen 1.9.1