LatticeYangMills
su3.h
Go to the documentation of this file.
1 /******************************************************************************
2 *
3 * MIT License
4 *
5 * Copyright (c) 2018 Giovanni Pederiva
6 *
7 * Permission is hereby granted, free of charge, to any person obtaining a copy
8 * of this software and associated documentation files (the "Software"), to deal
9 * in the Software without restriction, including without limitation the rights
10 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11 * copies of the Software, and to permit persons to whom the Software is
12 * furnished to do so, subject to the following conditions:
13 *
14 * The above copyright notice and this permission notice shall be included in
15 * all copies or substantial portions of the Software.
16 *
17 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
23 * SOFTWARE.
24 ******************************************************************************/
25 
34 #pragma once
35 #include <array>
36 #include "operators.hpp"
37 
53 struct SU3
54  : tao::operators::commutative_addable< SU3 >,
55  tao::operators::subtractable< SU3 >,
56  tao::operators::multipliable< SU3 >,
57  tao::operators::addable_left< SU3, double >,
58  tao::operators::subtractable_left< SU3, double >,
59  tao::operators::addable< SU3, double >,
60  tao::operators::subtractable< SU3, double >,
61  tao::operators::multipliable< SU3, double >
62 {
63  std::array<double, 18> mat;
64 
65  // constructors
66  SU3 () noexcept;
67  SU3 (double value) noexcept;
68  SU3 (const SU3 &source) noexcept;
69  SU3 (SU3 && source) noexcept;
70 
71  // Operations between elements
72  inline SU3& operator=(const SU3& other) noexcept;
73  inline SU3& operator=(SU3&& other) noexcept;
74 
75  inline SU3& operator+=( const SU3& other ) noexcept;
76  inline SU3& operator+=( SU3&& other ) noexcept;
77 
78  inline SU3& operator-=( const SU3& other ) noexcept;
79  inline SU3& operator-=( SU3&& other ) noexcept;
80 
81  inline SU3& operator*=( const SU3& other ) noexcept;
82  inline SU3& operator*=( SU3&& other ) noexcept;
83 
84  // operations with a scalar
85  inline SU3& operator+=( const double scalar ) noexcept;
86  inline SU3& operator-=( const double scalar ) noexcept;
87  inline SU3& operator*=( const double scalar ) noexcept;
88 
89  // Setters
90  void setSU3Identity();
91  void setSU3Zero();
92  void setSU3Random();
93 
94  // Traces
95  double realTrace();
96  double imagTrace();
97 
98  // Get the exponential of the matrix
99  SU3 exp();
100 
101  // Print Function
102  void printSU3();
103 };
104 
105 // assignment
106 SU3& SU3::operator=(SU3&& other) noexcept{
107  mat = std::move(other.mat);
108  return *this;
109 }
110 
111 SU3& SU3::operator=(const SU3& other) noexcept{
112  mat = other.mat;
113  return *this;
114 }
115 
116 // sum
117 SU3& SU3::operator+=(SU3&& other) noexcept{
118  for(int i = 0; i < 18; i++)
119  mat[i] += other.mat[i];
120  return *this;
121 }
122 
123 SU3& SU3::operator+=(const SU3& other) noexcept{
124  for(int i = 0; i < 18; i++)
125  mat[i] += other.mat[i];
126  return *this;
127 }
128 
129 // subtraction
130 SU3& SU3::operator-=(SU3&& other) noexcept{
131  for(int i = 0; i < 18; i++)
132  mat[i] -= other.mat[i];
133  return *this;
134 }
135 
136 SU3& SU3::operator-=(const SU3& other) noexcept{
137  for(int i = 0; i < 18; i++)
138  mat[i] -= other.mat[i];
139  return *this;
140 }
141 
142 // multiplication
143 SU3& SU3::operator*=(const SU3& other) noexcept{
144  SU3 result;
145  //REAL
146  result.mat[0] = mat[0] * other.mat[0] - mat[1]*other.mat[1] + mat[2] * other.mat[6] - mat[3]*other.mat[7] + mat[4] * other.mat[12] - mat[5]*other.mat[13];
147  result.mat[2] = mat[0] * other.mat[2] - mat[1]*other.mat[3] + mat[2] * other.mat[8] - mat[3]*other.mat[9] + mat[4] * other.mat[14] - mat[5]*other.mat[15];
148  result.mat[4] = mat[0] * other.mat[4] - mat[1]*other.mat[5] + mat[2] * other.mat[10] - mat[3]*other.mat[11] + mat[4] * other.mat[16] - mat[5]*other.mat[17];
149  result.mat[6] = mat[6] * other.mat[0] - mat[7]*other.mat[1] + mat[8] * other.mat[6] - mat[9]*other.mat[7] + mat[10] * other.mat[12] - mat[11]*other.mat[13];
150  result.mat[8] = mat[6] * other.mat[2] - mat[7]*other.mat[3] + mat[8] * other.mat[8] - mat[9]*other.mat[9] + mat[10] * other.mat[14] - mat[11]*other.mat[15];
151  result.mat[10] = mat[6] * other.mat[4] - mat[7]*other.mat[5] + mat[8] * other.mat[10] - mat[9]*other.mat[11] + mat[10] * other.mat[16] - mat[11]*other.mat[17];
152  result.mat[12] = mat[12] * other.mat[0] - mat[13]*other.mat[1] + mat[14] * other.mat[6] - mat[15]*other.mat[7] + mat[16] * other.mat[12] - mat[17]*other.mat[13];
153  result.mat[14] = mat[12] * other.mat[2] - mat[13]*other.mat[3] + mat[14] * other.mat[8] - mat[15]*other.mat[9] + mat[16] * other.mat[14] - mat[17]*other.mat[15];
154  result.mat[16] = mat[12] * other.mat[4] - mat[13]*other.mat[5] + mat[14] * other.mat[10] - mat[15]*other.mat[11] + mat[16] * other.mat[16] - mat[17]*other.mat[17];
155 
156 
157  // COMPLEX
158  result.mat[1] = mat[0] * other.mat[1] + mat[1]*other.mat[0] + mat[2] * other.mat[7] + mat[3]*other.mat[6] + mat[4] * other.mat[13] + mat[5]*other.mat[12];
159  result.mat[3] = mat[0] * other.mat[3] + mat[1]*other.mat[2] + mat[2] * other.mat[9] + mat[3]*other.mat[8] + mat[4] * other.mat[15] + mat[5]*other.mat[14];
160  result.mat[5] = mat[0] * other.mat[5] + mat[1]*other.mat[4] + mat[2] * other.mat[11] + mat[3]*other.mat[10] + mat[4] * other.mat[17] + mat[5]*other.mat[16];
161  result.mat[7] = mat[6] * other.mat[1] + mat[7]*other.mat[0] + mat[8] * other.mat[7] + mat[9]*other.mat[6] + mat[10] * other.mat[13] + mat[11]*other.mat[12];
162  result.mat[9] = mat[6] * other.mat[3] + mat[7]*other.mat[2] + mat[8] * other.mat[9] + mat[9]*other.mat[8] + mat[10] * other.mat[15] + mat[11]*other.mat[14];
163  result.mat[11] = mat[6] * other.mat[5] + mat[7]*other.mat[4] + mat[8] * other.mat[11] + mat[9]*other.mat[10] + mat[10] * other.mat[17] + mat[11]*other.mat[16];
164  result.mat[13] = mat[12] * other.mat[1] + mat[13]*other.mat[0] + mat[14] * other.mat[7] + mat[15]*other.mat[6] + mat[16] * other.mat[13] + mat[17]*other.mat[12];
165  result.mat[15] = mat[12] * other.mat[3] + mat[13]*other.mat[2] + mat[14] * other.mat[9] + mat[15]*other.mat[8] + mat[16] * other.mat[15] + mat[17]*other.mat[14];
166  result.mat[17] = mat[12] * other.mat[5] + mat[13]*other.mat[4] + mat[14] * other.mat[11] + mat[15]*other.mat[10] + mat[16] * other.mat[17] + mat[17]*other.mat[16];
167 
168  /*
169  for(int i = 0; i < 3; i++){
170  for(int j = 0; j < 3; j++){
171  c.mat[2 * (3*i + j)] = 0;
172  c.mat[2 * (3*i + j) + 1] = 0;
173  for(int k = 0; k < 3; k++){
174  c.mat[2 * (3*i + j)] += a.mat[2 * (3*i + k)] * b.mat[2 * (3*k + j)]
175  -a.mat[2 * (3*i + k) + 1] * b.mat[2 * (3*k + j) + 1];
176  c.mat[2 * (3*i + j) + 1] += a.mat[2 * (3*i + k)] * b.mat[2 * (3*k + j) + 1]
177  +a.mat[2 * (3*i + k) + 1] * b.mat[2 * (3*k + j)];
178 
179  }
180  }
181  }*/
182  return *this = result;
183 }
184 
185 SU3& SU3::operator*=(SU3&& other) noexcept{
186  SU3 result;
187  //REAL
188  result.mat[0] = mat[0] * other.mat[0] - mat[1]*other.mat[1] + mat[2] * other.mat[6] - mat[3]*other.mat[7] + mat[4] * other.mat[12] - mat[5]*other.mat[13];
189  result.mat[2] = mat[0] * other.mat[2] - mat[1]*other.mat[3] + mat[2] * other.mat[8] - mat[3]*other.mat[9] + mat[4] * other.mat[14] - mat[5]*other.mat[15];
190  result.mat[4] = mat[0] * other.mat[4] - mat[1]*other.mat[5] + mat[2] * other.mat[10] - mat[3]*other.mat[11] + mat[4] * other.mat[16] - mat[5]*other.mat[17];
191  result.mat[6] = mat[6] * other.mat[0] - mat[7]*other.mat[1] + mat[8] * other.mat[6] - mat[9]*other.mat[7] + mat[10] * other.mat[12] - mat[11]*other.mat[13];
192  result.mat[8] = mat[6] * other.mat[2] - mat[7]*other.mat[3] + mat[8] * other.mat[8] - mat[9]*other.mat[9] + mat[10] * other.mat[14] - mat[11]*other.mat[15];
193  result.mat[10] = mat[6] * other.mat[4] - mat[7]*other.mat[5] + mat[8] * other.mat[10] - mat[9]*other.mat[11] + mat[10] * other.mat[16] - mat[11]*other.mat[17];
194  result.mat[12] = mat[12] * other.mat[0] - mat[13]*other.mat[1] + mat[14] * other.mat[6] - mat[15]*other.mat[7] + mat[16] * other.mat[12] - mat[17]*other.mat[13];
195  result.mat[14] = mat[12] * other.mat[2] - mat[13]*other.mat[3] + mat[14] * other.mat[8] - mat[15]*other.mat[9] + mat[16] * other.mat[14] - mat[17]*other.mat[15];
196  result.mat[16] = mat[12] * other.mat[4] - mat[13]*other.mat[5] + mat[14] * other.mat[10] - mat[15]*other.mat[11] + mat[16] * other.mat[16] - mat[17]*other.mat[17];
197 
198 
199  // COMPLEX
200  result.mat[1] = mat[0] * other.mat[1] + mat[1]*other.mat[0] + mat[2] * other.mat[7] + mat[3]*other.mat[6] + mat[4] * other.mat[13] + mat[5]*other.mat[12];
201  result.mat[3] = mat[0] * other.mat[3] + mat[1]*other.mat[2] + mat[2] * other.mat[9] + mat[3]*other.mat[8] + mat[4] * other.mat[15] + mat[5]*other.mat[14];
202  result.mat[5] = mat[0] * other.mat[5] + mat[1]*other.mat[4] + mat[2] * other.mat[11] + mat[3]*other.mat[10] + mat[4] * other.mat[17] + mat[5]*other.mat[16];
203  result.mat[7] = mat[6] * other.mat[1] + mat[7]*other.mat[0] + mat[8] * other.mat[7] + mat[9]*other.mat[6] + mat[10] * other.mat[13] + mat[11]*other.mat[12];
204  result.mat[9] = mat[6] * other.mat[3] + mat[7]*other.mat[2] + mat[8] * other.mat[9] + mat[9]*other.mat[8] + mat[10] * other.mat[15] + mat[11]*other.mat[14];
205  result.mat[11] = mat[6] * other.mat[5] + mat[7]*other.mat[4] + mat[8] * other.mat[11] + mat[9]*other.mat[10] + mat[10] * other.mat[17] + mat[11]*other.mat[16];
206  result.mat[13] = mat[12] * other.mat[1] + mat[13]*other.mat[0] + mat[14] * other.mat[7] + mat[15]*other.mat[6] + mat[16] * other.mat[13] + mat[17]*other.mat[12];
207  result.mat[15] = mat[12] * other.mat[3] + mat[13]*other.mat[2] + mat[14] * other.mat[9] + mat[15]*other.mat[8] + mat[16] * other.mat[15] + mat[17]*other.mat[14];
208  result.mat[17] = mat[12] * other.mat[5] + mat[13]*other.mat[4] + mat[14] * other.mat[11] + mat[15]*other.mat[10] + mat[16] * other.mat[17] + mat[17]*other.mat[16];
209 
210  /*
211  for(int i = 0; i < 3; i++){
212  for(int j = 0; j < 3; j++){
213  c.mat[2 * (3*i + j)] = 0;
214  c.mat[2 * (3*i + j) + 1] = 0;
215  for(int k = 0; k < 3; k++){
216  c.mat[2 * (3*i + j)] += a.mat[2 * (3*i + k)] * b.mat[2 * (3*k + j)]
217  -a.mat[2 * (3*i + k) + 1] * b.mat[2 * (3*k + j) + 1];
218  c.mat[2 * (3*i + j) + 1] += a.mat[2 * (3*i + k)] * b.mat[2 * (3*k + j) + 1]
219  +a.mat[2 * (3*i + k) + 1] * b.mat[2 * (3*k + j)];
220 
221  }
222  }
223  }*/
224  return *this = result;
225 }
226 
227 
228 // sum with scalar
229 SU3& SU3::operator+=( const double scalar ) noexcept{
230  for(auto& i : mat)
231  i += scalar;
232  return *this;
233 }
234 
235 // subtraction with scalar
236 SU3& SU3::operator-=( const double scalar ) noexcept{
237  for(auto& i : mat)
238  i -= scalar;
239  return *this;
240 }
241 
242 // sum with scalar
243 SU3& SU3::operator*=( const double scalar ) noexcept{
244  for(auto& i : mat)
245  i *= scalar;
246  return *this;
247 }
248 
249 SU3 operator~(const SU3& a);
250 SU3 operator~(SU3&& a);
251 // Get the exponential of the matrix
252 SU3 exp(const SU3& Q);
253 SU3 exp(SU3&& Q);
254 
255 SU3 getRandomTransformation(double epsilon);
256 double getMultSumTrace(SU3& a, SU3& b);
257 double getMultSumTrace(SU3&& a, SU3& b);
Implementation of a class to perform arithmetics between links.
Definition: su3.h:53
double getMultSumTrace(SU3 &a, SU3 &b)
returns the real trace of the product of two SU3 matrices a and b
Definition: su3.cpp:118