From b1d19316100032cd889b949b02897feab685fbd5 Mon Sep 17 00:00:00 2001
From: Matthias Frey <matthias.frey@psi.ch>
Date: Wed, 15 Jul 2020 14:01:46 +0200
Subject: [PATCH] RealDiracMatrix: remove unused functions, change storage type
of ublas vector and fix bug in function 'transform'
---
src/Distribution/RealDiracMatrix.cpp | 124 ++++++---------------
src/Distribution/RealDiracMatrix.h | 50 ++++-----
src/Distribution/matrix_vector_operation.h | 6 +-
3 files changed, 56 insertions(+), 124 deletions(-)
diff --git a/src/Distribution/RealDiracMatrix.cpp b/src/Distribution/RealDiracMatrix.cpp
index 681527b9f..8b1cdd7bc 100644
--- a/src/Distribution/RealDiracMatrix.cpp
+++ b/src/Distribution/RealDiracMatrix.cpp
@@ -29,7 +29,7 @@
#include <string>
#include "matrix_vector_operation.h"
-RealDiracMatrix::RealDiracMatrix() : NumOfRDMs(16), DimOfRDMs(4) {};
+RealDiracMatrix::RealDiracMatrix() {};
typename RealDiracMatrix::sparse_matrix_t
RealDiracMatrix::getRDM(short i) {
@@ -59,56 +59,6 @@ RealDiracMatrix::getRDM(short i) {
}
-typename RealDiracMatrix::vector_t
-RealDiracMatrix::decompose(const matrix_t& M) {
- /*
- * Formula (11) from paper:
- * Geometrical method of decoupling
- */
- if (M.size1() != 4 && M.size1() != M.size2()) {
- throw OpalException("RealDiracMatrix::decompose()",
- "Wrong matrix dimensions.");
- }
-
- vector_t coeffs(16);
-
- matrix_t left(4,4), right(4,4), rdm2(4,4);
-
- double inv32 = 1.0 / 32.0;
-
- // do it with iterators
- for (short i = 0; i < 16; ++i) {
- sparse_matrix_t rdm = getRDM(i);
- left = boost::numeric::ublas::prod(M,rdm);
- right = boost::numeric::ublas::prod(rdm,M);
- rdm2 = boost::numeric::ublas::prod(rdm,rdm);
- left *= inv32;
- right *= inv32;
- left += right;
- coeffs(i) = matt_boost::trace(rdm2) * matt_boost::trace(left);
- }
- return coeffs;
-}
-
-
-typename RealDiracMatrix::matrix_t
-RealDiracMatrix::combine(const vector_t& coeffs) {
- if (coeffs.size() > 16) {
- throw OpalException("RealDiracMatrix::combine()",
- "Wrong size of coefficient vector.");
- }
-
- // initialize a 4x4 zero matrix
- matrix_t M = boost::numeric::ublas::zero_matrix<double>(4,4);
-
- // evaluate linear combination
- for (short i = 0; i < 16; ++i)
- M += coeffs(i)*getRDM(i);
-
- return M;
-}
-
-
void RealDiracMatrix::diagonalize(matrix_t& Ms, sparse_matrix_t& R, sparse_matrix_t& invR) {
// R and invR store the total transformation
@@ -129,39 +79,39 @@ void RealDiracMatrix::diagonalize(matrix_t& Ms, sparse_matrix_t& R, sparse_matri
};
// 1. Transformation with \gamma_{0}
- P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
- E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
- B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
- mr = boost::numeric::ublas::inner_prod(E,B); // formula (31), paper: Geometrical method of decoupling
- mg = boost::numeric::ublas::inner_prod(B,P); // formula (31), paper: Geometrical method of decoupling
+ P = vector_t({ mult(1), mult(2), mult(3)});
+ E = vector_t({ mult(4), mult(5), mult(6)});
+ B = vector_t({-mult(7), -mult(8), -mult(9)});
+ mr = boost::numeric::ublas::inner_prod(E, B); // formula (31), paper: Geometrical method of decoupling
+ mg = boost::numeric::ublas::inner_prod(B, P); // formula (31), paper: Geometrical method of decoupling
transform(Ms, 0, 0.5 * std::atan2(mg,mr), R, invR);
// 2. Transformation with \gamma_{7}
- eps = -mult(0);
- P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
- E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
- B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
- b = eps*B+matt_boost::cross_prod(E,P); // formula (32), paper: Geometrical method of decoupling
+ eps = -mult(0);
+ P = vector_t({ mult(1), mult(2), mult(3)});
+ E = vector_t({ mult(4), mult(5), mult(6)});
+ B = vector_t({-mult(7), -mult(8), -mult(9)});
+ b = eps * B + matt_boost::cross_prod(E, P); // formula (32), paper: Geometrical method of decoupling
transform(Ms, 7, 0.5 * std::atan2(b(2), b(1)), R, invR);
// 3. Transformation with \gamma_{9}
- eps = -mult(0);
- P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
- E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
- B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
- b = eps*B+matt_boost::cross_prod(E,P);
+ eps = -mult(0);
+ P = vector_t({ mult(1), mult(2), mult(3)});
+ E = vector_t({ mult(4), mult(5), mult(6)});
+ B = vector_t({-mult(7), -mult(8), -mult(9)});
+ b = eps * B + matt_boost::cross_prod(E, P);
transform(Ms, 9, -0.5 * std::atan2(b(0), b(1)), R, invR);
// 4. Transformation with \gamma_{2}
- eps = -mult(0);
- P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
- E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
- B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
- mr = boost::numeric::ublas::inner_prod(E,B);
- b = eps*B+matt_boost::cross_prod(E,P);
+ eps = -mult(0);
+ P = vector_t({ mult(1), mult(2), mult(3)});
+ E = vector_t({ mult(4), mult(5), mult(6)});
+ B = vector_t({-mult(7), -mult(8), -mult(9)});
+ mr = boost::numeric::ublas::inner_prod(E, B);
+ b = eps * B + matt_boost::cross_prod(E, P);
// Transformation distinction made according to function "rdm_Decouple_F"
// in rdm.c of Dr. Christian Baumgarten
@@ -172,18 +122,18 @@ void RealDiracMatrix::diagonalize(matrix_t& Ms, sparse_matrix_t& R, sparse_matri
transform(Ms, 2, 0.5 * std::atanh(b(1) / mr), R, invR);
}
- eps = -mult(0);
- P(0) = mult(1); P(1) = mult(2); P(2) = mult(3);
- E(0) = mult(4); E(1) = mult(5); E(2) = mult(6);
- B(0) = -mult(7); B(1) = -mult(8); B(2) = -mult(9);
+ eps = -mult(0);
+ P = vector_t({ mult(1), mult(2), mult(3)});
+ E = vector_t({ mult(4), mult(5), mult(6)});
+ B = vector_t({-mult(7), -mult(8), -mult(9)});
// formula (31), paper: Geometrical method of decoupling
- mr = boost::numeric::ublas::inner_prod(E,B);
- mg = boost::numeric::ublas::inner_prod(B,P);
- mb = boost::numeric::ublas::inner_prod(E,P);
+ mr = boost::numeric::ublas::inner_prod(E, B);
+ mg = boost::numeric::ublas::inner_prod(B, P);
+ mb = boost::numeric::ublas::inner_prod(E, P);
- double P2 = boost::numeric::ublas::inner_prod(P,P);
- double E2 = boost::numeric::ublas::inner_prod(E,E);
+ double P2 = boost::numeric::ublas::inner_prod(P, P);
+ double E2 = boost::numeric::ublas::inner_prod(E, E);
// 5. Transform with \gamma_{0}
transform(Ms, 0, 0.25 * std::atan2(mb,0.5 * (E2 - P2)), R, invR);
@@ -225,16 +175,6 @@ RealDiracMatrix::symplex(const matrix_t& M) {
}
-typename RealDiracMatrix::matrix_t
-RealDiracMatrix::cosymplex(const matrix_t& M) {
- /*
- * formula (16), p, 3
- */
- sparse_matrix_t rdm0 = getRDM(0);
- return 0.5 * (M - matt_boost::gemmm<matrix_t>(rdm0,boost::numeric::ublas::trans(M),rdm0));
-}
-
-
void RealDiracMatrix::transform(matrix_t& M, short i, double phi,
sparse_matrix_t& Rtot, sparse_matrix_t& invRtot)
{
@@ -258,7 +198,7 @@ void RealDiracMatrix::transform(short i, double phi,
if (phi) { // if phi == 0 --> nothing happens, since R and invR would be identity_matrix matrix
sparse_matrix_t I = boost::numeric::ublas::identity_matrix<double>(4);
- if (i < 7 && i != 0 && i > 10 && i != 14) {
+ if ((i < 7 && i != 0) || (i > 10 && i != 14)) {
R = I * std::cosh(phi) + getRDM(i) * std::sinh(phi);
invR = I * std::cosh(phi) - getRDM(i) * std::sinh(phi);
} else {
diff --git a/src/Distribution/RealDiracMatrix.h b/src/Distribution/RealDiracMatrix.h
index 1c21a74de..ae6fa40ca 100644
--- a/src/Distribution/RealDiracMatrix.h
+++ b/src/Distribution/RealDiracMatrix.h
@@ -40,7 +40,10 @@ public:
/// Dense matrix type definition
typedef boost::numeric::ublas::matrix<double> matrix_t;
/// Dense vector type definition
- typedef boost::numeric::ublas::vector<double> vector_t;
+ typedef boost::numeric::ublas::vector<
+ double,
+ std::vector<double>
+ > vector_t;
/// Default constructor (sets only NumOfRDMs and DimOfRDMs)
RealDiracMatrix();
@@ -51,21 +54,6 @@ public:
*/
sparse_matrix_t getRDM(short);
- /*!
- * Decomposes a real-valued 4x4 matrix into a linear combination
- * and returns a vector containing the coefficients
- * @param M an arbitrary real-valued 4x4 matrix
- */
- vector_t decompose(const matrix_t&);
-
- /*!
- * Takes a vector of coefficients, evaluates the linear combination of
- * RDMs with these coefficients and returns a 4x4 matrix
- *
- * @param coeffs is a vector of coefficients (at most length NumOfRDMs)
- */
- matrix_t combine(const vector_t&);
-
/*!
* Brings a 4x4 symplex matrix into Hamilton form and
* computes the transformation matrix and its inverse
@@ -91,22 +79,10 @@ public:
*/
matrix_t symplex(const matrix_t&);
- /*!
- * @param M 4x4 real-valued matrix
- * @returns the cosymplex part of a 4x4 real-valued matrix
- */
- matrix_t cosymplex(const matrix_t&);
-
- /// The number of real Dirac matrices
- short NumOfRDMs;
- /// The matrix dimension (4x4)
- short DimOfRDMs;
-
-
private:
- ///
/*!
* Applies a rotation to the matrix M by a given angle
+ *
* @param M is the matrix to be transformed
* @param i is the i-th RDM used for transformation
* @param phi is the angle of rotation
@@ -115,8 +91,24 @@ private:
*/
void transform(matrix_t&, short, double, sparse_matrix_t&, sparse_matrix_t&);
+ /*!
+ * Obtain transformation matrices.
+ *
+ * @param i is the i-th RDM used for transformation
+ * @param phi is the angle of rotation
+ * @param R is a reference to the transformation matrix
+ * @param invR is a reference to the inverse of the transformation matrix
+ */
void transform(short, double, sparse_matrix_t&, sparse_matrix_t&);
+ /*!
+ * Update quantites to decouple the sigma-matrix
+ *
+ * @param sigma current state of sigma-matrix
+ * @param i is the i-th step
+ * @param R is a reference to the transformation matrix
+ * @param invR is a reference to the inverse of the transformation matrix
+ */
void update(matrix_t& sigma, short i, sparse_matrix_t& R, sparse_matrix_t& invR);
};
diff --git a/src/Distribution/matrix_vector_operation.h b/src/Distribution/matrix_vector_operation.h
index 12be15405..45f7aec46 100644
--- a/src/Distribution/matrix_vector_operation.h
+++ b/src/Distribution/matrix_vector_operation.h
@@ -55,10 +55,10 @@ namespace matt_boost {
}
/// Computes the cross product \f$ v_{1}\times v_{2}\f$ of two vectors in \f$ \mathbb{R}^{3} \f$
- template<class V>
+ template<class V, class A>
BOOST_UBLAS_INLINE
- ublas::vector<V> cross_prod(ublas::vector<V>& v1, ublas::vector<V>& v2) {
- ublas::vector<V> v(v1.size());
+ ublas::vector<V, A> cross_prod(ublas::vector<V, A>& v1, ublas::vector<V, A>& v2) {
+ ublas::vector<V, A> v(v1.size());
if (v1.size() == v2.size() && v1.size() == 3) {
v(0) = v1(1) * v2(2) - v1(2) * v2(1);
v(1) = v1(2) * v2(0) - v1(0) * v2(2);
--
GitLab