Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mat_SU_N.h
Go to the documentation of this file.
1 
15 #ifndef MAT_SU_N_INCLUDED
16 #define MAT_SU_N_INCLUDED
17 
18 #include <valarray>
19 //#include <complex>
20 
21 class RandomNumbers; // H.Matsufuru added [15 Mar 2012].
22 
23 namespace SU_N {
24  class Mat_SU_N {
25  private:
26  int m_Nc;
27  std::valarray<double> va;
28  public:
29  Mat_SU_N(int Nc, double r = 0.0) : m_Nc(Nc) { va.resize(2 * m_Nc * m_Nc, r); }
30 
31  Mat_SU_N& dag();
32  Mat_SU_N& ht();
33  Mat_SU_N& ah();
34  Mat_SU_N& at();
35  Mat_SU_N& unit();
36  Mat_SU_N& zero();
37  Mat_SU_N& xI();
38  Mat_SU_N& reunit();
39 
40  Mat_SU_N& set_random(RandomNumbers *); // H.Matsufuru added [15 Mar 2012].
41 
42  int nc() const { return m_Nc; }
43 
44  Mat_SU_N& operator-();
45 
46  Mat_SU_N& operator=(const double&);
47 
48  // Mat_SU_N& operator=(const std::complex<double>&);
49 
50  Mat_SU_N& operator+=(const Mat_SU_N&);
51  Mat_SU_N& operator+=(const double&);
52 
53  // Mat_SU_N& operator+=(const std::complex<double>&);
54 
55  Mat_SU_N& operator-=(const Mat_SU_N&);
56  Mat_SU_N& operator-=(const double&);
57 
58  // Mat_SU_N& operator-=(const std::complex<double>&);
59 
60  Mat_SU_N& operator*=(const Mat_SU_N&);
61  Mat_SU_N& operator*=(const double&);
62 
63  // Mat_SU_N& operator*=(const std::complex<double>&);
64 
65  Mat_SU_N& operator/=(const double&);
66 
67  // Mat_SU_N& operator/=(const std::complex<double>&);
68 
69  int size() const { return va.size(); }
70  double r(int c) const { return va[2 * c]; }
71  double i(int c) const { return va[2 * c + 1]; }
72 
73  double r(int c1, int c2) const { return r(m_Nc * c1 + c2); }
74  double i(int c1, int c2) const { return i(m_Nc * c1 + c2); }
75 
76  void setr(int c, const double& re) { va[2 * c] = re; }
77  void seti(int c, const double& im) { va[2 * c + 1] = im; }
78 
79  void setr(int c1, int c2, const double& re)
80  {
81  setr(m_Nc * c1 + c2, re);
82  }
83 
84  void seti(int c1, int c2, const double& im)
85  {
86  seti(m_Nc * c1 + c2, im);
87  }
88 
89  void set(int c, double re, const double& im)
90  {
91  va[2 * c] = re;
92  va[2 * c + 1] = im;
93  }
94 
95  void set(int c1, int c2, const double& re, const double& im)
96  {
97  set(m_Nc * c1 + c2, re, im);
98  }
99 
100  void add(int c, const double& re, const double& im)
101  {
102  va[2 * c] += re;
103  va[2 * c + 1] += im;
104  }
105 
106  void add(int c1, int c2, const double& re, const double& im)
107  {
108  add(m_Nc * c1 + c2, re, im);
109  }
110 
111  double norm2() // H.Matsufuru added [15 Mar 2012]
112  {
113  double t = 0.0;
114 
115  for (int i = 0; i < 2 * m_Nc * m_Nc; ++i) { t += va[i] * va[i]; }
116  return t;
117  }
118 
119  void mult_nn(const Mat_SU_N& u1, const Mat_SU_N& u2)
120  {
121  for (int a = 0; a < m_Nc; ++a) {
122  for (int b = 0; b < m_Nc; ++b) {
123  va[2 * (b + m_Nc * a)] = 0.0;
124  va[2 * (b + m_Nc * a) + 1] = 0.0;
125  for (int c = 0; c < m_Nc; ++c) {
126  va[2 * (b + m_Nc * a)] += u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (b + m_Nc * c)]
127  - u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (b + m_Nc * c) + 1];
128  va[2 * (b + m_Nc * a) + 1] += u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (b + m_Nc * c)]
129  + u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (b + m_Nc * c) + 1];
130  }
131  }
132  }
133  }
134 
135  void multadd_nn(const Mat_SU_N& u1, const Mat_SU_N& u2)
136  {
137  for (int a = 0; a < m_Nc; ++a) {
138  for (int b = 0; b < m_Nc; ++b) {
139  for (int c = 0; c < m_Nc; ++c) {
140  va[2 * (b + m_Nc * a)] += u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (b + m_Nc * c)]
141  - u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (b + m_Nc * c) + 1];
142  va[2 * (b + m_Nc * a) + 1] += u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (b + m_Nc * c)]
143  + u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (b + m_Nc * c) + 1];
144  }
145  }
146  }
147  }
148 
149  void mult_nd(const Mat_SU_N& u1, const Mat_SU_N& u2)
150  {
151  for (int a = 0; a < m_Nc; ++a) {
152  for (int b = 0; b < m_Nc; ++b) {
153  va[2 * (b + m_Nc * a)] = 0.0;
154  va[2 * (b + m_Nc * a) + 1] = 0.0;
155  for (int c = 0; c < m_Nc; ++c) {
156  va[2 * (b + m_Nc * a)] += u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (c + m_Nc * b)]
157  + u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (c + m_Nc * b) + 1];
158  va[2 * (b + m_Nc * a) + 1] += u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (c + m_Nc * b)]
159  - u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (c + m_Nc * b) + 1];
160  }
161  }
162  }
163  }
164 
165  void multadd_nd(const Mat_SU_N& u1, const Mat_SU_N& u2)
166  {
167  for (int a = 0; a < m_Nc; ++a) {
168  for (int b = 0; b < m_Nc; ++b) {
169  for (int c = 0; c < m_Nc; ++c) {
170  va[2 * (b + m_Nc * a)] += u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (c + m_Nc * b)]
171  + u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (c + m_Nc * b) + 1];
172  va[2 * (b + m_Nc * a) + 1] += u1.va[2 * (c + m_Nc * a) + 1] * u2.va[2 * (c + m_Nc * b)]
173  - u1.va[2 * (c + m_Nc * a)] * u2.va[2 * (c + m_Nc * b) + 1];
174  }
175  }
176  }
177  }
178 
179  void mult_dn(const Mat_SU_N& u1, const Mat_SU_N& u2)
180  {
181  for (int a = 0; a < m_Nc; ++a) {
182  for (int b = 0; b < m_Nc; ++b) {
183  va[2 * (b + m_Nc * a)] = 0.0;
184  va[2 * (b + m_Nc * a) + 1] = 0.0;
185  for (int c = 0; c < m_Nc; ++c) {
186  va[2 * (b + m_Nc * a)] += u1.va[2 * (a + m_Nc * c)] * u2.va[2 * (b + m_Nc * c)]
187  + u1.va[2 * (a + m_Nc * c) + 1] * u2.va[2 * (b + m_Nc * c) + 1];
188  va[2 * (b + m_Nc * a) + 1] -= u1.va[2 * (a + m_Nc * c) + 1] * u2.va[2 * (b + m_Nc * c)]
189  - u1.va[2 * (a + m_Nc * c)] * u2.va[2 * (b + m_Nc * c) + 1];
190  }
191  }
192  }
193  }
194 
195  void multadd_dn(const Mat_SU_N& u1, const Mat_SU_N& u2)
196  {
197  for (int a = 0; a < m_Nc; ++a) {
198  for (int b = 0; b < m_Nc; ++b) {
199  for (int c = 0; c < m_Nc; ++c) {
200  va[2 * (b + m_Nc * a)] += u1.va[2 * (a + m_Nc * c)] * u2.va[2 * (b + m_Nc * c)]
201  + u1.va[2 * (a + m_Nc * c) + 1] * u2.va[2 * (b + m_Nc * c) + 1];
202  va[2 * (b + m_Nc * a) + 1] -= u1.va[2 * (a + m_Nc * c) + 1] * u2.va[2 * (b + m_Nc * c)]
203  - u1.va[2 * (a + m_Nc * c)] * u2.va[2 * (b + m_Nc * c) + 1];
204  }
205  }
206  }
207  }
208 
209  void zcopy(double re, double im, const Mat_SU_N& v)
210  {
211  for (int cc = 0; cc < m_Nc * m_Nc; ++cc) {
212  va[2 * cc] = re * v.va[2 * cc] - im * v.va[2 * cc + 1];
213  va[2 * cc + 1] = re * v.va[2 * cc + 1] + im * v.va[2 * cc];
214  }
215  }
216 
217  void zaxpy(double re, double im, const Mat_SU_N& v)
218  {
219  for (int cc = 0; cc < m_Nc * m_Nc; ++cc) {
220  va[2 * cc] += re * v.va[2 * cc] - im * v.va[2 * cc + 1];
221  va[2 * cc + 1] += re * v.va[2 * cc + 1] + im * v.va[2 * cc];
222  }
223  }
224  };
225 
226  inline Mat_SU_N& Mat_SU_N::dag()
227  {
228  for (int a = 0; a < m_Nc; ++a) {
229  for (int b = a; b < m_Nc; ++b) {
230  double re = va[2 * (m_Nc * a + b)];
231  double im = va[2 * (m_Nc * a + b) + 1];
232 
233  va[2 * (m_Nc * a + b)] = va[2 * (m_Nc * b + a)];
234  va[2 * (m_Nc * a + b) + 1] = -va[2 * (m_Nc * b + a) + 1];
235 
236  va[2 * (m_Nc * b + a)] = re;
237  va[2 * (m_Nc * b + a) + 1] = -im;
238  }
239  }
240  return *this;
241  }
242 
243 
244  inline Mat_SU_N& Mat_SU_N::ht() // hermitian traceless
245  {
246  for (int a = 0; a < m_Nc; ++a) {
247  for (int b = a + 1; b < m_Nc; ++b) {
248  double re = va[2 * (m_Nc * a + b)] + va[2 * (m_Nc * b + a)];
249  double im = va[2 * (m_Nc * a + b) + 1] - va[2 * (m_Nc * b + a) + 1];
250 
251  va[2 * (m_Nc * a + b)] = 0.5 * re;
252  va[2 * (m_Nc * a + b) + 1] = 0.5 * im;
253 
254  va[2 * (m_Nc * b + a)] = 0.5 * re;
255  va[2 * (m_Nc * b + a) + 1] = -0.5 * im;
256  }
257  }
258  double tr = 0.0;
259  for (int cc = 0; cc < m_Nc; ++cc) {
260  tr += va[2 * (m_Nc * cc + cc)];
261  }
262  tr = tr / m_Nc;
263  for (int cc = 0; cc < m_Nc; ++cc) {
264  va[2 * (m_Nc * cc + cc)] -= tr;
265  va[2 * (m_Nc * cc + cc) + 1] = 0.0;
266  }
267  return *this;
268  }
269 
270 
272  inline Mat_SU_N& Mat_SU_N::at()
273  {
274  for (int a = 0; a < m_Nc; ++a) {
275  for (int b = a + 1; b < m_Nc; ++b) {
276  double re = va[2 * (m_Nc * a + b)] - va[2 * (m_Nc * b + a)];
277  double im = va[2 * (m_Nc * a + b) + 1] + va[2 * (m_Nc * b + a) + 1];
278 
279  va[2 * (m_Nc * a + b)] = 0.5 * re;
280  va[2 * (m_Nc * a + b) + 1] = 0.5 * im;
281 
282  va[2 * (m_Nc * b + a)] = -0.5 * re;
283  va[2 * (m_Nc * b + a) + 1] = 0.5 * im;
284  }
285  }
286  double tr = 0.0;
287  for (int cc = 0; cc < m_Nc; ++cc) {
288  tr += va[2 * (m_Nc * cc + cc) + 1];
289  }
290  tr = tr / m_Nc;
291  for (int cc = 0; cc < m_Nc; ++cc) {
292  va[2 * (m_Nc * cc + cc)] = 0.0;
293  va[2 * (m_Nc * cc + cc) + 1] -= tr;
294  }
295  return *this;
296  }
297 
298 
300  inline Mat_SU_N& Mat_SU_N::ah()
301  {
302  for (int a = 0; a < m_Nc; ++a) {
303  for (int b = a; b < m_Nc; ++b) {
304  double re = va[2 * (m_Nc * a + b)] - va[2 * (m_Nc * b + a)];
305  double im = va[2 * (m_Nc * a + b) + 1] + va[2 * (m_Nc * b + a) + 1];
306  va[2 * (m_Nc * a + b)] = 0.5 * re;
307  va[2 * (m_Nc * a + b) + 1] = 0.5 * im;
308  va[2 * (m_Nc * b + a)] = -0.5 * re;
309  va[2 * (m_Nc * b + a) + 1] = 0.5 * im;
310  }
311  }
312  return *this;
313  }
314 
315 
316  inline Mat_SU_N& Mat_SU_N::unit()
317  {
318  va = 0.0;
319  for (int c = 0; c < m_Nc; ++c) {
320  va[2 * (m_Nc + 1) * c] = 1.0;
321  }
322  return *this;
323  }
324 
325 
326  inline Mat_SU_N& Mat_SU_N::zero()
327  {
328  va = 0.0;
329  return *this;
330  }
331 
332 
333  inline Mat_SU_N& Mat_SU_N::xI()
334  {
335  for (int c = 0; c < va.size() / 2; ++c) {
336  double tmp = va[2 * c];
337  va[2 * c] = -va[2 * c + 1];
338  va[2 * c + 1] = tmp;
339  }
340  return *this;
341  }
342 
343 
344  inline Mat_SU_N& Mat_SU_N::operator-()
345  {
346  va = -va;
347  return *this;
348  }
349 
350 
351  inline Mat_SU_N& Mat_SU_N::operator+=(const Mat_SU_N& rhs)
352  {
353  va += rhs.va;
354  return *this;
355  }
356 
357 
358  inline Mat_SU_N& Mat_SU_N::operator-=(const Mat_SU_N& rhs)
359  {
360  va -= rhs.va;
361  return *this;
362  }
363 
364 
365  inline Mat_SU_N& Mat_SU_N::operator=(const double& rhs)
366  {
367  va = rhs;
368  return *this;
369  }
370 
371 
372  inline Mat_SU_N& Mat_SU_N::operator+=(const double& rhs)
373  {
374  va += rhs;
375  return *this;
376  }
377 
378 
379  inline Mat_SU_N& Mat_SU_N::operator-=(const double& rhs)
380  {
381  va -= rhs;
382  return *this;
383  }
384 
385 
386  inline Mat_SU_N& Mat_SU_N::operator*=(const Mat_SU_N& rhs)
387  {
388  std::valarray<double> tmp(2 * m_Nc * m_Nc);
389 
390  for (int a = 0; a < m_Nc; ++a) {
391  for (int b = 0; b < m_Nc; ++b) {
392  tmp[2 * (m_Nc * a + b)] = 0.0;
393  tmp[2 * (m_Nc * a + b) + 1] = 0.0;
394  for (int c = 0; c < m_Nc; ++c) {
395  tmp[2 * (m_Nc * a + b)] += va[2 * (m_Nc * a + c)] * rhs.va[2 * (m_Nc * c + b)]
396  - va[2 * (m_Nc * a + c) + 1] * rhs.va[2 * (m_Nc * c + b) + 1];
397  tmp[2 * (m_Nc * a + b) + 1] += va[2 * (m_Nc * a + c) + 1] * rhs.va[2 * (m_Nc * c + b)]
398  + va[2 * (m_Nc * a + c)] * rhs.va[2 * (m_Nc * c + b) + 1];
399  }
400  }
401  }
402  va = tmp;
403  return *this;
404  }
405 
406 
407  inline Mat_SU_N& Mat_SU_N::operator*=(const double& rhs)
408  {
409  va *= rhs;
410  return *this;
411  }
412 
413 
414 /*
415 inline Mat_SU_N& Mat_SU_N::operator*=(const std::complex<double>& rhs){
416  std::valarray<double> tmp = va;
417  for(int c = 0; c < va.size()/2; ++c){
418  va[2*c ] = (tmp[2*c]*real(rhs) -tmp[2*c+1]*imag(rhs));
419  va[2*c+1] = (tmp[2*c]*imag(rhs) +tmp[2*c+1]*real(rhs));
420  }
421  return *this;
422 }
423 */
424  inline Mat_SU_N& Mat_SU_N::operator/=(const double& rhs)
425  {
426  va /= rhs;
427  return *this;
428  }
429 
430 
431  inline double ReTr(const Mat_SU_N& m)
432  {
433  int Nc = m.nc();
434  double tr = 0.0;
435 
436  for (int c = 0; c < Nc; ++c) {
437  tr += m.r(c, c);
438  }
439  return tr;
440  }
441 
442 
443  inline double ImTr(const Mat_SU_N& m)
444  {
445  int Nc = m.nc();
446  double tr = 0.0;
447 
448  for (int c = 0; c < Nc; ++c) {
449  tr += m.i(c, c);
450  }
451  return tr;
452  }
453 
454 
455  inline Mat_SU_N dag(const Mat_SU_N& u)
456  {
457  int Nc = u.nc();
458  Mat_SU_N tmp(Nc);
459 
460  for (int a = 0; a < Nc; a++) {
461  for (int b = 0; b < Nc; b++) {
462  tmp.set(a, b, u.r(b, a), -u.i(b, a));
463  }
464  }
465  return tmp;
466  }
467 
468 
469  inline Mat_SU_N xI(const Mat_SU_N& u)
470  {
471  int Nc = u.nc();
472  Mat_SU_N tmp(Nc);
473 
474  for (int c = 0; c < u.size() / 2; ++c) {
475  tmp.set(c, -u.i(c), u.r(c));
476  }
477  return tmp;
478  }
479 
480 
481  inline Mat_SU_N operator+(const Mat_SU_N& m1, const Mat_SU_N& m2)
482  {
483  return Mat_SU_N(m1) += m2;
484  }
485 
486 
487  inline Mat_SU_N operator-(const Mat_SU_N& m1, const Mat_SU_N& m2)
488  {
489  return Mat_SU_N(m1) -= m2;
490  }
491 
492 
493  inline Mat_SU_N operator*(const Mat_SU_N& m1, const Mat_SU_N& m2)
494  {
495  return Mat_SU_N(m1) *= m2;
496  }
497 
498 
499  inline Mat_SU_N reunit(const Mat_SU_N& m)
500  {
501  Mat_SU_N tmp = m;
502 
503  tmp.reunit();
504  return tmp;
505  }
506 } // namespace SU_N
507 #endif