use std::iter::{FromIterator, FusedIterator};
use std::{
ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign},
vec::Vec,
};
use na::{ComplexField, Matrix3, SVector};
use rayon::iter::FromParallelIterator;
use rayon::prelude::*;
#[cfg(feature = "serde-serialize")]
use serde::{Deserialize, Serialize};
use super::{
lattice::{Direction, LatticeCyclic, LatticeElementToIndex, LatticeLink, LatticePoint},
su3,
su3::GENERATORS,
thread::{run_pool_parallel_vec_with_initializations_mutable, ThreadError},
utils::levi_civita,
CMatrix3, Complex, Real, Vector8, I,
};
#[derive(Debug, Copy, Clone, PartialEq)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct Su3Adjoint {
data: Vector8<Real>,
}
#[allow(clippy::len_without_is_empty)]
impl Su3Adjoint {
pub const fn new(data: Vector8<Real>) -> Self {
Self { data }
}
pub fn new_from_array(data: [Real; 8]) -> Self {
Su3Adjoint::new(Vector8::from(data))
}
pub const fn data(&self) -> &Vector8<Real> {
&self.data
}
pub const fn as_vector(&self) -> &Vector8<Real> {
self.data()
}
pub fn as_vector_mut(&mut self) -> &mut Vector8<Real> {
self.data_mut()
}
pub fn to_matrix(self) -> Matrix3<na::Complex<Real>> {
self.data
.into_iter()
.enumerate()
.map(|(pos, el)| *GENERATORS[pos] * na::Complex::<Real>::from(el))
.sum()
}
pub fn to_su3(self) -> Matrix3<na::Complex<Real>> {
su3::su3_exp_i(self)
}
pub fn exp(self) -> Matrix3<na::Complex<Real>> {
su3::su3_exp_r(self)
}
pub fn random<Rng>(rng: &mut Rng, d: &impl rand_distr::Distribution<Real>) -> Self
where
Rng: rand::Rng + ?Sized,
{
Self {
data: Vector8::<Real>::from_fn(|_, _| d.sample(rng)),
}
}
#[inline]
pub fn trace_squared(&self) -> Real {
self.data.iter().map(|el| el * el).sum::<Real>() / 2_f64
}
#[inline]
pub fn t(&self) -> Real {
-0.5_f64 * self.trace_squared()
}
#[inline]
pub fn d(&self) -> na::Complex<Real> {
self.to_matrix().determinant() * I
}
#[allow(clippy::unused_self)]
pub fn len(&self) -> usize {
self.data.len()
}
pub const fn len_const() -> usize {
8
}
pub fn iter(&self) -> impl Iterator<Item = &Real> + ExactSizeIterator + FusedIterator {
self.data.iter()
}
pub fn iter_mut(
&mut self,
) -> impl Iterator<Item = &mut Real> + ExactSizeIterator + FusedIterator {
self.data.iter_mut()
}
pub fn data_mut(&mut self) -> &mut Vector8<Real> {
&mut self.data
}
}
impl AsRef<Vector8<f64>> for Su3Adjoint {
fn as_ref(&self) -> &Vector8<f64> {
self.as_vector()
}
}
impl AsMut<Vector8<f64>> for Su3Adjoint {
fn as_mut(&mut self) -> &mut Vector8<f64> {
self.as_vector_mut()
}
}
impl<'a> IntoIterator for &'a Su3Adjoint {
type IntoIter = <&'a Vector8<Real> as IntoIterator>::IntoIter;
type Item = &'a Real;
fn into_iter(self) -> Self::IntoIter {
self.data.iter()
}
}
impl<'a> IntoIterator for &'a mut Su3Adjoint {
type IntoIter = <&'a mut Vector8<Real> as IntoIterator>::IntoIter;
type Item = &'a mut Real;
fn into_iter(self) -> Self::IntoIter {
self.data.iter_mut()
}
}
impl AddAssign for Su3Adjoint {
fn add_assign(&mut self, other: Self) {
self.data += other.data();
}
}
impl Add<Su3Adjoint> for Su3Adjoint {
type Output = Self;
fn add(mut self, rhs: Self) -> Self::Output {
self += rhs;
self
}
}
impl Add<&Su3Adjoint> for Su3Adjoint {
type Output = Self;
fn add(self, rhs: &Self) -> Self::Output {
self + *rhs
}
}
impl Add<Su3Adjoint> for &Su3Adjoint {
type Output = Su3Adjoint;
fn add(self, rhs: Su3Adjoint) -> Self::Output {
rhs + self
}
}
impl Add<&Su3Adjoint> for &Su3Adjoint {
type Output = Su3Adjoint;
fn add(self, rhs: &Su3Adjoint) -> Self::Output {
self + *rhs
}
}
impl MulAssign<f64> for Su3Adjoint {
fn mul_assign(&mut self, rhs: f64) {
self.data *= rhs;
}
}
impl Mul<Real> for Su3Adjoint {
type Output = Self;
fn mul(mut self, rhs: Real) -> Self::Output {
self *= rhs;
self
}
}
impl Mul<&Real> for Su3Adjoint {
type Output = Self;
fn mul(self, rhs: &Real) -> Self::Output {
self * (*rhs)
}
}
impl Mul<Real> for &Su3Adjoint {
type Output = Su3Adjoint;
fn mul(self, rhs: Real) -> Self::Output {
*self * rhs
}
}
impl Mul<&Real> for &Su3Adjoint {
type Output = Su3Adjoint;
fn mul(self, rhs: &Real) -> Self::Output {
*self * rhs
}
}
impl Mul<Su3Adjoint> for Real {
type Output = Su3Adjoint;
fn mul(self, rhs: Su3Adjoint) -> Self::Output {
rhs * self
}
}
impl Mul<&Su3Adjoint> for Real {
type Output = Su3Adjoint;
fn mul(self, rhs: &Su3Adjoint) -> Self::Output {
rhs * self
}
}
impl Mul<Su3Adjoint> for &Real {
type Output = Su3Adjoint;
fn mul(self, rhs: Su3Adjoint) -> Self::Output {
rhs * self
}
}
impl Mul<&Su3Adjoint> for &Real {
type Output = Su3Adjoint;
fn mul(self, rhs: &Su3Adjoint) -> Self::Output {
rhs * self
}
}
impl DivAssign<f64> for Su3Adjoint {
fn div_assign(&mut self, rhs: f64) {
self.data /= rhs;
}
}
impl DivAssign<&f64> for Su3Adjoint {
fn div_assign(&mut self, rhs: &f64) {
self.data /= *rhs;
}
}
impl Div<Real> for Su3Adjoint {
type Output = Self;
fn div(mut self, rhs: Real) -> Self::Output {
self /= rhs;
self
}
}
impl Div<&Real> for Su3Adjoint {
type Output = Self;
fn div(self, rhs: &Real) -> Self::Output {
self / (*rhs)
}
}
impl Div<Real> for &Su3Adjoint {
type Output = Su3Adjoint;
fn div(self, rhs: Real) -> Self::Output {
*self / rhs
}
}
impl Div<&Real> for &Su3Adjoint {
type Output = Su3Adjoint;
fn div(self, rhs: &Real) -> Self::Output {
*self / rhs
}
}
impl SubAssign for Su3Adjoint {
fn sub_assign(&mut self, other: Self) {
self.data -= other.data();
}
}
impl Sub<Su3Adjoint> for Su3Adjoint {
type Output = Self;
fn sub(mut self, rhs: Self) -> Self::Output {
self -= rhs;
self
}
}
impl Sub<&Su3Adjoint> for Su3Adjoint {
type Output = Self;
fn sub(self, rhs: &Self) -> Self::Output {
self - *rhs
}
}
impl Sub<Su3Adjoint> for &Su3Adjoint {
type Output = Su3Adjoint;
fn sub(self, rhs: Su3Adjoint) -> Self::Output {
rhs - self
}
}
impl Sub<&Su3Adjoint> for &Su3Adjoint {
type Output = Su3Adjoint;
fn sub(self, rhs: &Su3Adjoint) -> Self::Output {
*self - rhs
}
}
impl Neg for Su3Adjoint {
type Output = Self;
fn neg(self) -> Self::Output {
Su3Adjoint::new(-self.data)
}
}
impl Neg for &Su3Adjoint {
type Output = Su3Adjoint;
fn neg(self) -> Self::Output {
Su3Adjoint::new(-self.data())
}
}
impl Default for Su3Adjoint {
fn default() -> Self {
Su3Adjoint::new(Vector8::from_element(0_f64))
}
}
impl std::fmt::Display for Su3Adjoint {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "{}", self.to_matrix())
}
}
impl Index<usize> for Su3Adjoint {
type Output = Real;
fn index(&self, pos: usize) -> &Self::Output {
&self.data[pos]
}
}
impl IndexMut<usize> for Su3Adjoint {
fn index_mut(&mut self, pos: usize) -> &mut Self::Output {
&mut self.data[pos]
}
}
impl From<Vector8<Real>> for Su3Adjoint {
fn from(v: Vector8<Real>) -> Self {
Su3Adjoint::new(v)
}
}
impl From<Su3Adjoint> for Vector8<Real> {
fn from(v: Su3Adjoint) -> Self {
v.data
}
}
impl From<&Su3Adjoint> for Vector8<Real> {
fn from(v: &Su3Adjoint) -> Self {
v.data
}
}
impl From<[Real; 8]> for Su3Adjoint {
fn from(v: [Real; 8]) -> Self {
Su3Adjoint::new_from_array(v)
}
}
#[derive(Debug, PartialEq, Clone)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct LinkMatrix {
data: Vec<Matrix3<na::Complex<Real>>>,
}
impl LinkMatrix {
pub const fn new(data: Vec<Matrix3<na::Complex<Real>>>) -> Self {
Self { data }
}
pub const fn data(&self) -> &Vec<Matrix3<na::Complex<Real>>> {
&self.data
}
pub fn data_mut(&mut self) -> &mut Vec<Matrix3<na::Complex<Real>>> {
&mut self.data
}
pub const fn as_vec(&self) -> &Vec<Matrix3<na::Complex<Real>>> {
self.data()
}
pub fn as_vec_mut(&mut self) -> &mut Vec<Matrix3<na::Complex<Real>>> {
self.data_mut()
}
pub fn as_slice(&self) -> &[Matrix3<na::Complex<Real>>] {
self.data()
}
pub fn as_slice_mut(&mut self) -> &mut [Matrix3<na::Complex<Real>>] {
&mut self.data
}
pub fn new_determinist<Rng: rand::Rng + ?Sized, const D: usize>(
l: &LatticeCyclic<D>,
rng: &mut Rng,
) -> Self {
let mut data = Vec::with_capacity(l.number_of_canonical_links_space());
for _ in l.get_links() {
let matrix = su3::random_su3(rng);
data.push(matrix);
}
Self { data }
}
pub fn new_random_threaded<const D: usize>(
l: &LatticeCyclic<D>,
number_of_thread: usize,
) -> Result<Self, ThreadError> {
if number_of_thread == 0 {
return Err(ThreadError::ThreadNumberIncorrect);
}
else if number_of_thread == 1 {
let mut rng = rand::thread_rng();
return Ok(LinkMatrix::new_determinist(l, &mut rng));
}
let data = run_pool_parallel_vec_with_initializations_mutable(
l.get_links(),
&(),
&|rng, _, _| su3::random_su3(rng),
rand::thread_rng,
number_of_thread,
l.number_of_canonical_links_space(),
l,
&CMatrix3::zeros(),
)?;
Ok(Self { data })
}
pub fn new_cold<const D: usize>(l: &LatticeCyclic<D>) -> Self {
Self {
data: vec![CMatrix3::identity(); l.number_of_canonical_links_space()],
}
}
pub fn matrix<const D: usize>(
&self,
link: &LatticeLink<D>,
l: &LatticeCyclic<D>,
) -> Option<Matrix3<na::Complex<Real>>> {
let link_c = l.into_canonical(*link);
let matrix = self.data.get(link_c.to_index(l))?;
if link.is_dir_negative() {
Some(matrix.adjoint())
}
else {
Some(*matrix)
}
}
pub fn sij<const D: usize>(
&self,
point: &LatticePoint<D>,
dir_i: &Direction<D>,
dir_j: &Direction<D>,
lattice: &LatticeCyclic<D>,
) -> Option<Matrix3<na::Complex<Real>>> {
let u_j = self.matrix(&LatticeLink::new(*point, *dir_j), lattice)?;
let point_pj = lattice.add_point_direction(*point, dir_j);
let u_i_p_j = self.matrix(&LatticeLink::new(point_pj, *dir_i), lattice)?;
let point_pi = lattice.add_point_direction(*point, dir_i);
let u_j_pi_d = self
.matrix(&LatticeLink::new(point_pi, *dir_j), lattice)?
.adjoint();
Some(u_j * u_i_p_j * u_j_pi_d)
}
pub fn pij<const D: usize>(
&self,
point: &LatticePoint<D>,
dir_i: &Direction<D>,
dir_j: &Direction<D>,
lattice: &LatticeCyclic<D>,
) -> Option<Matrix3<na::Complex<Real>>> {
let s_ij = self.sij(point, dir_i, dir_j, lattice)?;
let u_i = self.matrix(&LatticeLink::new(*point, *dir_i), lattice)?;
Some(u_i * s_ij.adjoint())
}
#[allow(clippy::as_conversions)] pub fn average_trace_plaquette<const D: usize>(
&self,
lattice: &LatticeCyclic<D>,
) -> Option<na::Complex<Real>> {
if lattice.number_of_canonical_links_space() != self.len() {
return None;
}
let sum = lattice
.get_points()
.par_bridge()
.map(|point| {
Direction::positive_directions()
.iter()
.map(|dir_i| {
Direction::positive_directions()
.iter()
.filter(|dir_j| dir_i.index() < dir_j.index())
.map(|dir_j| {
self.pij(&point, dir_i, dir_j, lattice).map(|el| el.trace())
})
.sum::<Option<na::Complex<Real>>>()
})
.sum::<Option<na::Complex<Real>>>()
})
.sum::<Option<na::Complex<Real>>>()?;
let number_of_directions = (D * (D - 1)) / 2;
let number_of_plaquette = (lattice.number_of_points() * number_of_directions) as f64;
Some(sum / number_of_plaquette)
}
pub fn clover<const D: usize>(
&self,
point: &LatticePoint<D>,
dir_i: &Direction<D>,
dir_j: &Direction<D>,
lattice: &LatticeCyclic<D>,
) -> Option<CMatrix3> {
Some(
self.pij(point, dir_i, dir_j, lattice)?
+ self.pij(point, dir_j, &-dir_i, lattice)?
+ self.pij(point, &-dir_i, &-dir_j, lattice)?
+ self.pij(point, &-dir_j, dir_i, lattice)?,
)
}
pub fn f_mu_nu<const D: usize>(
&self,
point: &LatticePoint<D>,
dir_i: &Direction<D>,
dir_j: &Direction<D>,
lattice: &LatticeCyclic<D>,
) -> Option<CMatrix3> {
let m = self.clover(point, dir_i, dir_j, lattice)?
- self.clover(point, dir_j, dir_i, lattice)?;
Some(m / Complex::from(8_f64 * lattice.size() * lattice.size()))
}
pub fn magnetic_field_vec<const D: usize>(
&self,
point: &LatticePoint<D>,
lattice: &LatticeCyclic<D>,
) -> Option<SVector<CMatrix3, D>> {
let mut vec = SVector::<CMatrix3, D>::zeros();
for dir in &Direction::<D>::positive_directions() {
vec[dir.index()] = self.magnetic_field(point, dir, lattice)?;
}
Some(vec)
}
pub fn magnetic_field<const D: usize>(
&self,
point: &LatticePoint<D>,
dir: &Direction<D>,
lattice: &LatticeCyclic<D>,
) -> Option<CMatrix3> {
let sum = Direction::<D>::positive_directions()
.iter()
.map(|dir_i| {
Direction::<D>::positive_directions()
.iter()
.map(|dir_j| {
let f_mn = self.f_mu_nu(point, dir_i, dir_j, lattice)?;
let lc = Complex::from(
levi_civita(&[dir.index(), dir_i.index(), dir_j.index()]).to_f64(),
);
Some(f_mn * lc)
})
.sum::<Option<CMatrix3>>()
})
.sum::<Option<CMatrix3>>()?;
Some(sum / Complex::new(0_f64, 2_f64))
}
pub fn magnetic_field_link<const D: usize>(
&self,
link: &LatticeLink<D>,
lattice: &LatticeCyclic<D>,
) -> Option<Matrix3<na::Complex<Real>>> {
self.magnetic_field(link.pos(), link.dir(), lattice)
}
pub fn len(&self) -> usize {
self.data.len()
}
pub fn is_empty(&self) -> bool {
self.data.is_empty()
}
pub fn normalize(&mut self) {
self.data.par_iter_mut().for_each(|el| {
su3::orthonormalize_matrix_mut(el);
});
}
pub fn iter(&self) -> impl Iterator<Item = &CMatrix3> + ExactSizeIterator + FusedIterator {
self.data.iter()
}
pub fn iter_mut(
&mut self,
) -> impl Iterator<Item = &mut CMatrix3> + ExactSizeIterator + FusedIterator {
self.data.iter_mut()
}
}
impl AsRef<Vec<CMatrix3>> for LinkMatrix {
fn as_ref(&self) -> &Vec<CMatrix3> {
self.as_vec()
}
}
impl AsMut<Vec<CMatrix3>> for LinkMatrix {
fn as_mut(&mut self) -> &mut Vec<CMatrix3> {
self.as_vec_mut()
}
}
impl AsRef<[CMatrix3]> for LinkMatrix {
fn as_ref(&self) -> &[CMatrix3] {
self.as_slice()
}
}
impl AsMut<[CMatrix3]> for LinkMatrix {
fn as_mut(&mut self) -> &mut [CMatrix3] {
self.as_slice_mut()
}
}
impl<'a> IntoIterator for &'a LinkMatrix {
type IntoIter = <&'a Vec<CMatrix3> as IntoIterator>::IntoIter;
type Item = &'a CMatrix3;
fn into_iter(self) -> Self::IntoIter {
self.data.iter()
}
}
impl<'a> IntoIterator for &'a mut LinkMatrix {
type IntoIter = <&'a mut Vec<CMatrix3> as IntoIterator>::IntoIter;
type Item = &'a mut CMatrix3;
fn into_iter(self) -> Self::IntoIter {
self.data.iter_mut()
}
}
impl Index<usize> for LinkMatrix {
type Output = CMatrix3;
fn index(&self, pos: usize) -> &Self::Output {
&self.data[pos]
}
}
impl IndexMut<usize> for LinkMatrix {
fn index_mut(&mut self, pos: usize) -> &mut Self::Output {
&mut self.data[pos]
}
}
impl<A> FromIterator<A> for LinkMatrix
where
Vec<CMatrix3>: FromIterator<A>,
{
fn from_iter<T>(iter: T) -> Self
where
T: IntoIterator<Item = A>,
{
Self::new(Vec::from_iter(iter))
}
}
impl<A> FromParallelIterator<A> for LinkMatrix
where
Vec<CMatrix3>: FromParallelIterator<A>,
A: Send,
{
fn from_par_iter<I>(par_iter: I) -> Self
where
I: IntoParallelIterator<Item = A>,
{
Self::new(Vec::from_par_iter(par_iter))
}
}
impl<T> ParallelExtend<T> for LinkMatrix
where
Vec<CMatrix3>: ParallelExtend<T>,
T: Send,
{
fn par_extend<I>(&mut self, par_iter: I)
where
I: IntoParallelIterator<Item = T>,
{
self.data.par_extend(par_iter);
}
}
impl<A> Extend<A> for LinkMatrix
where
Vec<CMatrix3>: Extend<A>,
{
fn extend<T>(&mut self, iter: T)
where
T: IntoIterator<Item = A>,
{
self.data.extend(iter);
}
}
#[derive(Debug, PartialEq, Clone)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct EField<const D: usize> {
data: Vec<SVector<Su3Adjoint, D>>, }
impl<const D: usize> EField<D> {
pub fn new(data: Vec<SVector<Su3Adjoint, D>>) -> Self {
Self { data }
}
pub const fn data(&self) -> &Vec<SVector<Su3Adjoint, D>> {
&self.data
}
pub fn data_mut(&mut self) -> &mut Vec<SVector<Su3Adjoint, D>> {
&mut self.data
}
pub const fn as_vec(&self) -> &Vec<SVector<Su3Adjoint, D>> {
self.data()
}
pub fn as_slice(&self) -> &[SVector<Su3Adjoint, D>] {
&self.data
}
pub fn as_slice_mut(&mut self) -> &mut [SVector<Su3Adjoint, D>] {
&mut self.data
}
pub fn len(&self) -> usize {
self.data.len()
}
pub fn new_determinist<Rng: rand::Rng + ?Sized>(
l: &LatticeCyclic<D>,
rng: &mut Rng,
d: &impl rand_distr::Distribution<Real>,
) -> Self {
let mut data = Vec::with_capacity(l.number_of_points());
for _ in l.get_points() {
data.push(SVector::<Su3Adjoint, D>::from_fn(|_, _| {
Su3Adjoint::random(rng, d)
}));
}
Self { data }
}
pub fn new_random(l: &LatticeCyclic<D>, d: &impl rand_distr::Distribution<Real>) -> Self {
let mut rng = rand::thread_rng();
EField::new_determinist(l, &mut rng, d)
}
pub fn new_cold(l: &LatticeCyclic<D>) -> Self {
let p1 = Su3Adjoint::new_from_array([0_f64; 8]);
Self {
data: vec![SVector::<Su3Adjoint, D>::from_element(p1); l.number_of_points()],
}
}
pub fn e_vec(
&self,
point: &LatticePoint<D>,
l: &LatticeCyclic<D>,
) -> Option<&SVector<Su3Adjoint, D>> {
self.data.get(point.to_index(l))
}
pub fn e_field(
&self,
point: &LatticePoint<D>,
dir: &Direction<D>,
l: &LatticeCyclic<D>,
) -> Option<&Su3Adjoint> {
let value = self.e_vec(point, l);
match value {
Some(vec) => vec.get(dir.index()),
None => None,
}
}
pub fn is_empty(&self) -> bool {
self.data.is_empty()
}
#[inline]
pub fn gauss(
&self,
link_matrix: &LinkMatrix,
point: &LatticePoint<D>,
lattice: &LatticeCyclic<D>,
) -> Option<CMatrix3> {
if lattice.number_of_points() != self.len()
|| lattice.number_of_canonical_links_space() != link_matrix.len()
{
return None;
}
Direction::positive_directions()
.iter()
.map(|dir| {
let e_i = self.e_field(point, dir, lattice)?;
let u_mi = link_matrix.matrix(&LatticeLink::new(*point, -*dir), lattice)?;
let p_mi = lattice.add_point_direction(*point, &-dir);
let e_m_i = self.e_field(&p_mi, dir, lattice)?;
Some(e_i.to_matrix() - u_mi * e_m_i.to_matrix() * u_mi.adjoint())
})
.sum::<Option<CMatrix3>>()
}
#[inline]
pub fn gauss_sum_div(
&self,
link_matrix: &LinkMatrix,
lattice: &LatticeCyclic<D>,
) -> Option<Real> {
if lattice.number_of_points() != self.len()
|| lattice.number_of_canonical_links_space() != link_matrix.len()
{
return None;
}
lattice
.get_points()
.par_bridge()
.map(|point| {
self.gauss(link_matrix, &point, lattice).map(|el| {
(su3::GENERATORS.iter().copied().sum::<CMatrix3>() * el)
.trace()
.abs()
})
})
.sum::<Option<Real>>()
}
#[allow(clippy::as_conversions)] #[inline]
pub fn project_to_gauss(
&self,
link_matrix: &LinkMatrix,
lattice: &LatticeCyclic<D>,
) -> Option<Self> {
const NUMBER_FOR_LOOP: usize = 4;
if lattice.number_of_points() != self.len()
|| lattice.number_of_canonical_links_space() != link_matrix.len()
{
return None;
}
let mut return_val = self.project_to_gauss_step(link_matrix, lattice);
loop {
let val_dif = return_val.gauss_sum_div(link_matrix, lattice)?;
if val_dif.is_nan() {
return None;
}
if val_dif <= f64::EPSILON * (lattice.number_of_points() * 4 * 8 * 10) as f64 {
break;
}
for _ in 0_usize..NUMBER_FOR_LOOP {
return_val = return_val.project_to_gauss_step(link_matrix, lattice);
}
}
Some(return_val)
}
#[inline]
fn project_to_gauss_step(&self, link_matrix: &LinkMatrix, lattice: &LatticeCyclic<D>) -> Self {
const K: na::Complex<f64> = na::Complex::new(0.12_f64, 0_f64);
let data = lattice
.get_points()
.collect::<Vec<LatticePoint<D>>>()
.par_iter()
.map(|point| {
let e = self.e_vec(point, lattice).unwrap();
SVector::<_, D>::from_fn(|index_dir, _| {
let dir = Direction::<D>::positive_directions()[index_dir];
let u = link_matrix
.matrix(&LatticeLink::new(*point, dir), lattice)
.unwrap();
let gauss = self.gauss(link_matrix, point, lattice).unwrap();
let gauss_p = self
.gauss(
link_matrix,
&lattice.add_point_direction(*point, &dir),
lattice,
)
.unwrap();
Su3Adjoint::new(Vector8::from_fn(|index, _| {
2_f64
* (su3::GENERATORS[index]
* ((u * gauss * u.adjoint() * gauss_p - gauss) * K
+ su3::GENERATORS[index]
* na::Complex::from(e[dir.index()][index])))
.trace()
.real()
}))
})
})
.collect();
Self::new(data)
}
}
impl<const D: usize> AsRef<Vec<SVector<Su3Adjoint, D>>> for EField<D> {
fn as_ref(&self) -> &Vec<SVector<Su3Adjoint, D>> {
self.as_vec()
}
}
impl<const D: usize> AsMut<Vec<SVector<Su3Adjoint, D>>> for EField<D> {
fn as_mut(&mut self) -> &mut Vec<SVector<Su3Adjoint, D>> {
self.data_mut()
}
}
impl<const D: usize> AsRef<[SVector<Su3Adjoint, D>]> for EField<D> {
fn as_ref(&self) -> &[SVector<Su3Adjoint, D>] {
self.as_slice()
}
}
impl<const D: usize> AsMut<[SVector<Su3Adjoint, D>]> for EField<D> {
fn as_mut(&mut self) -> &mut [SVector<Su3Adjoint, D>] {
self.as_slice_mut()
}
}
impl<'a, const D: usize> IntoIterator for &'a EField<D> {
type IntoIter = <&'a Vec<SVector<Su3Adjoint, D>> as IntoIterator>::IntoIter;
type Item = &'a SVector<Su3Adjoint, D>;
fn into_iter(self) -> Self::IntoIter {
self.data.iter()
}
}
impl<'a, const D: usize> IntoIterator for &'a mut EField<D> {
type IntoIter = <&'a mut Vec<SVector<Su3Adjoint, D>> as IntoIterator>::IntoIter;
type Item = &'a mut SVector<Su3Adjoint, D>;
fn into_iter(self) -> Self::IntoIter {
self.data.iter_mut()
}
}
impl<const D: usize> Index<usize> for EField<D> {
type Output = SVector<Su3Adjoint, D>;
fn index(&self, pos: usize) -> &Self::Output {
&self.data[pos]
}
}
impl<const D: usize> IndexMut<usize> for EField<D> {
fn index_mut(&mut self, pos: usize) -> &mut Self::Output {
&mut self.data[pos]
}
}
impl<A, const D: usize> FromIterator<A> for EField<D>
where
Vec<SVector<Su3Adjoint, D>>: FromIterator<A>,
{
fn from_iter<T>(iter: T) -> Self
where
T: IntoIterator<Item = A>,
{
Self::new(Vec::from_iter(iter))
}
}
impl<A, const D: usize> FromParallelIterator<A> for EField<D>
where
Vec<SVector<Su3Adjoint, D>>: FromParallelIterator<A>,
A: Send,
{
fn from_par_iter<I>(par_iter: I) -> Self
where
I: IntoParallelIterator<Item = A>,
{
Self::new(Vec::from_par_iter(par_iter))
}
}
impl<T, const D: usize> ParallelExtend<T> for EField<D>
where
Vec<SVector<Su3Adjoint, D>>: ParallelExtend<T>,
T: Send,
{
fn par_extend<I>(&mut self, par_iter: I)
where
I: IntoParallelIterator<Item = T>,
{
self.data.par_extend(par_iter);
}
}
impl<A, const D: usize> Extend<A> for EField<D>
where
Vec<SVector<Su3Adjoint, D>>: Extend<A>,
{
fn extend<T>(&mut self, iter: T)
where
T: IntoIterator<Item = A>,
{
self.data.extend(iter);
}
}
#[cfg(test)]
mod test {
use approx::*;
use rand::SeedableRng;
use super::super::{lattice::*, Complex};
use super::*;
const EPSILON: f64 = 0.000_000_001_f64;
const SEED_RNG: u64 = 0x45_78_93_f4_4a_b0_67_f0;
#[test]
fn test_get_e_field_pos_neg() {
use super::super::lattice;
let l = LatticeCyclic::new(1_f64, 4).unwrap();
let e = EField::new(vec![SVector::<_, 4>::from([
Su3Adjoint::from([1_f64; 8]),
Su3Adjoint::from([2_f64; 8]),
Su3Adjoint::from([3_f64; 8]),
Su3Adjoint::from([2_f64; 8]),
])]);
assert_eq!(
e.e_field(
&LatticePoint::new([0, 0, 0, 0].into()),
&lattice::DirectionEnum::XPos.into(),
&l
),
e.e_field(
&LatticePoint::new([0, 0, 0, 0].into()),
&lattice::DirectionEnum::XNeg.into(),
&l
)
);
}
#[test]
#[allow(clippy::eq_op)]
#[allow(clippy::op_ref)]
fn test_su3_adj() {
let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
let d = rand::distributions::Uniform::from(-1_f64..1_f64);
for _ in 0_u32..100_u32 {
let v = Su3Adjoint::random(&mut rng, &d);
let m = v.to_matrix();
assert_abs_diff_eq!(
v.trace_squared(),
(m * m).trace().modulus(),
epsilon = EPSILON
);
assert_eq_complex!(
v.d(),
nalgebra::Complex::new(0_f64, 1_f64) * m.determinant(),
EPSILON
);
assert_eq_complex!(v.t(), -(m * m).trace() / Complex::from(2_f64), EPSILON);
let adj_1 = Su3Adjoint::default();
let adj_2 = Su3Adjoint::new_from_array([1_f64; 8]);
assert_eq!(adj_2, adj_2 + adj_1);
assert_eq!(adj_2, &adj_2 + &adj_1);
assert_eq!(adj_2, &adj_2 - &adj_1);
assert_eq!(adj_1, &adj_2 - &adj_2);
assert_eq!(adj_1, &adj_2 - adj_2);
assert_eq!(adj_1, adj_2 - &adj_2);
assert_eq!(adj_1, -&adj_1);
let adj_3 = Su3Adjoint::new_from_array([2_f64; 8]);
assert_eq!(adj_3, &adj_2 + &adj_2);
assert_eq!(adj_3, &adj_2 * &2_f64);
assert_eq!(adj_3, &2_f64 * &adj_2);
assert_eq!(adj_3, 2_f64 * adj_2);
assert_eq!(adj_3, &2_f64 * adj_2);
assert_eq!(adj_3, 2_f64 * &adj_2);
assert_eq!(adj_2, &adj_3 / &2_f64);
assert_eq!(adj_2, &adj_3 / 2_f64);
let mut adj_5 = Su3Adjoint::new_from_array([2_f64; 8]);
adj_5 /= &2_f64;
assert_eq!(adj_2, adj_5);
let adj_4 = Su3Adjoint::new_from_array([-1_f64; 8]);
assert_eq!(adj_2, -adj_4);
}
use crate::su3::su3_exp_r;
let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
let d = rand::distributions::Uniform::from(-1_f64..1_f64);
for _ in 0_u32..10_u32 {
let v = Su3Adjoint::random(&mut rng, &d);
assert_eq!(su3_exp_r(v), v.exp());
}
}
#[test]
fn link_matrix() {
let lattice = LatticeCyclic::<3>::new(1_f64, 4).unwrap();
match LinkMatrix::new_random_threaded(&lattice, 0) {
Err(ThreadError::ThreadNumberIncorrect) => {}
_ => panic!("unexpected output"),
}
let link_s = LinkMatrix::new_random_threaded(&lattice, 2);
assert!(link_s.is_ok());
let mut link = link_s.unwrap();
assert!(!link.is_empty());
let l2 = LinkMatrix::new(vec![]);
assert!(l2.is_empty());
let _: &[_] = link.as_ref();
let _: &Vec<_> = link.as_ref();
let _: &mut [_] = link.as_mut();
let _: &mut Vec<_> = link.as_mut();
let _ = link.iter();
let _ = link.iter_mut();
let _ = (&link).into_iter();
let _ = (&mut link).into_iter();
}
#[test]
fn e_field() {
let lattice = LatticeCyclic::<3>::new(1_f64, 4).unwrap();
let e_field_s = LinkMatrix::new_random_threaded(&lattice, 2);
assert!(e_field_s.is_ok());
let mut e_field = e_field_s.unwrap();
let _: &[_] = e_field.as_ref();
let _: &Vec<_> = e_field.as_ref();
let _: &mut [_] = e_field.as_mut();
let _: &mut Vec<_> = e_field.as_mut();
let _ = e_field.iter();
let _ = e_field.iter_mut();
let _ = (&e_field).into_iter();
let _ = (&mut e_field).into_iter();
}
#[test]
fn magnetic_field() {
let lattice = LatticeCyclic::<3>::new(1_f64, 4).unwrap();
let mut link_matrix = LinkMatrix::new_cold(&lattice);
let point = LatticePoint::from([0, 0, 0]);
let dir_x = Direction::new(0, true).unwrap();
let dir_y = Direction::new(1, true).unwrap();
let dir_z = Direction::new(2, true).unwrap();
let clover = link_matrix
.clover(&point, &dir_x, &dir_y, &lattice)
.unwrap();
assert_eq_matrix!(CMatrix3::identity() * Complex::from(4_f64), clover, EPSILON);
let f = link_matrix
.f_mu_nu(&point, &dir_x, &dir_y, &lattice)
.unwrap();
assert_eq_matrix!(CMatrix3::zeros(), f, EPSILON);
let b = link_matrix
.magnetic_field(&point, &dir_x, &lattice)
.unwrap();
assert_eq_matrix!(CMatrix3::zeros(), b, EPSILON);
let b_vec = link_matrix.magnetic_field_vec(&point, &lattice).unwrap();
for i in &b_vec {
assert_eq_matrix!(CMatrix3::zeros(), i, EPSILON);
}
link_matrix[0] = CMatrix3::identity() * Complex::new(0_f64, 1_f64);
let clover = link_matrix
.clover(&point, &dir_x, &dir_y, &lattice)
.unwrap();
assert_eq_matrix!(
CMatrix3::identity() * Complex::new(2_f64, 0_f64),
clover,
EPSILON
);
let clover = link_matrix
.clover(&point, &dir_y, &dir_x, &lattice)
.unwrap();
assert_eq_matrix!(
CMatrix3::identity() * Complex::new(2_f64, 0_f64),
clover,
EPSILON
);
let f = link_matrix
.f_mu_nu(&point, &dir_x, &dir_y, &lattice)
.unwrap();
assert_eq_matrix!(
CMatrix3::identity() * Complex::new(0_f64, 0_f64),
f,
EPSILON
);
let b = link_matrix
.magnetic_field(&point, &dir_x, &lattice)
.unwrap();
assert_eq_matrix!(CMatrix3::zeros(), b, EPSILON);
let b_vec = link_matrix.magnetic_field_vec(&point, &lattice).unwrap();
for i in &b_vec {
assert_eq_matrix!(CMatrix3::zeros(), i, EPSILON);
}
assert_eq_matrix!(
link_matrix
.magnetic_field_link(&LatticeLink::new(point, dir_x), &lattice)
.unwrap(),
b,
EPSILON
);
let mut link_matrix = LinkMatrix::new_cold(&lattice);
let link = LatticeLinkCanonical::new([1, 0, 0].into(), dir_y).unwrap();
link_matrix[link.to_index(&lattice)] = CMatrix3::identity() * Complex::new(0_f64, 1_f64);
let clover = link_matrix
.clover(&point, &dir_x, &dir_y, &lattice)
.unwrap();
assert_eq_matrix!(
CMatrix3::identity() * Complex::new(3_f64, 1_f64),
clover,
EPSILON
);
let clover = link_matrix
.clover(&point, &dir_y, &dir_x, &lattice)
.unwrap();
assert_eq_matrix!(
CMatrix3::identity() * Complex::new(3_f64, -1_f64),
clover,
EPSILON
);
let f = link_matrix
.f_mu_nu(&point, &dir_x, &dir_y, &lattice)
.unwrap();
assert_eq_matrix!(
CMatrix3::identity() * Complex::new(0_f64, 0.25_f64),
f,
EPSILON
);
let b = link_matrix
.magnetic_field(&point, &dir_x, &lattice)
.unwrap();
assert_eq_matrix!(CMatrix3::zeros(), b, EPSILON);
assert_eq_matrix!(
link_matrix
.magnetic_field_link(&LatticeLink::new(point, dir_x), &lattice)
.unwrap(),
b,
EPSILON
);
let b = link_matrix
.magnetic_field(&point, &dir_z, &lattice)
.unwrap();
assert_eq_matrix!(
link_matrix
.magnetic_field_link(&LatticeLink::new(point, dir_z), &lattice)
.unwrap(),
b,
EPSILON
);
assert_eq_matrix!(
CMatrix3::identity() * Complex::new(0.25_f64, 0_f64),
b,
EPSILON
);
let b_2 = link_matrix
.magnetic_field(&[4, 0, 0].into(), &dir_z, &lattice)
.unwrap();
assert_eq_matrix!(b, b_2, EPSILON);
let b_vec = link_matrix.magnetic_field_vec(&point, &lattice).unwrap();
for (index, m) in b_vec.iter().enumerate() {
if index == 2 {
assert_eq_matrix!(m, b, EPSILON);
}
else {
assert_eq_matrix!(CMatrix3::zeros(), m, EPSILON);
}
}
}
#[test]
fn test_len() {
let su3 = Su3Adjoint::new(nalgebra::SVector::<f64, 8>::zeros());
assert_eq!(Su3Adjoint::len_const(), su3.len());
}
}