#include <GPUClusterInfoAndMomentsCalculatorImpl.h>

Collaboration diagram for ClusterMomentsCalculator::RealSymmetricMatrixSolver:

Public Member Functions
CUDA_HOS_DEV	RealSymmetricMatrixSolver (const float a_orig, const float b_orig, const float c_orig, const float d_orig, const float e_orig, const float f_orig)
CUDA_HOS_DEV void	get_eigenvalues (float &e_1, float &e_2, float &e_3) const
	Calculate shifted and scaled eigenvalues of the matrix, in ascending value.
CUDA_HOS_DEV void	extract_one (const float eigenvalue, float(&res)[3], float(&representative)[3]) const
CUDA_HOS_DEV void	get_eigenvectors (float(&res)[3][3], const float e_1, const float e_2, const float e_3, const float epsilon=s_typical_epsilon) const
	Calculate the eigenvectors of the matrix, using the (possibly unscaled) eigenvalues `e_1`, `e_2`, `e_3` (in ascending order of magnitude) and `epsilon` to guard against undefined cases.
CUDA_HOS_DEV void	get_solution (float(&eigenvalues)[3], float(&eigenvectors)[3][3], bool rescale_and_reshift_values=true, const float epsilon=s_typical_epsilon)
	Get the full eigenvalues and eigenvectors for this matrix.

Public Attributes
float	a
float	b
float	c
float	d
float	e
float	f
float	shift
float	scale

Static Public Attributes
static constexpr float	s_typical_epsilon = std::numeric_limits<float>::epsilon()

Detailed Description

Definition at line 336 of file GPUClusterInfoAndMomentsCalculatorImpl.h.

Constructor & Destructor Documentation

◆ RealSymmetricMatrixSolver()

CUDA_HOS_DEV ClusterMomentsCalculator::RealSymmetricMatrixSolver::RealSymmetricMatrixSolver	(	const float	a_orig,
		const float	b_orig,
		const float	c_orig,
		const float	d_orig,
		const float	e_orig,
		const float	f_orig )

inline

Definition at line 346 of file GPUClusterInfoAndMomentsCalculatorImpl.h.

    {
      using namespace std;
 
      shift = CaloRecGPU::Helpers::sum_kahan_babushka_neumaier(a_orig, b_orig, c_orig) / 3.f;
      a = a_orig - shift;
      b = b_orig - shift;
      c = c_orig - shift;
      const float max_ab = max( fabsf(a),       fabsf(b)      );
      const float max_cd = max( fabsf(c),       fabsf(d_orig) );
      const float max_ef = max( fabsf(e_orig),  fabsf(f_orig) );
      scale = max(max_ab, max(max_cd, max_ef) );
      if (scale == 0.f)
        {
          scale = 1.f;
        }
      const float inv_scale = 1.0f / scale;
      a *= inv_scale;
      b *= inv_scale;
      c *= inv_scale;
      d = d_orig * inv_scale;
      e = e_orig * inv_scale;
      f = f_orig * inv_scale;
    }

Member Function Documentation

◆ extract_one()

CUDA_HOS_DEV void ClusterMomentsCalculator::RealSymmetricMatrixSolver::extract_one	(	const float	eigenvalue,
		float(&)	res[3],
		float(&)	representative[3] ) const

inline

Definition at line 438 of file GPUClusterInfoAndMomentsCalculatorImpl.h.

    {
      using namespace std;
 
      const float diag[3] = { a - eigenvalue,
                              b - eigenvalue,
                              c - eigenvalue
                            };
 
      float vec_1[3], vec_2[3];
 
      int   max_diag  = -1;
      float max_value = -1.f;
 
      for (int i = 0; i < 3; ++i)
        {
          const float this_diag = fabsf(diag[i]);
          if (this_diag >= max_value)
            {
              max_diag  = i;
              max_value = this_diag;
            }
        }
 
      if (max_diag == 0)
        {
          representative[0] = diag[0];
          representative[1] = d;
          representative[2] = f;
 
          vec_1[0] = d;
          vec_1[1] = diag[1];
          vec_1[2] = e;
 
          vec_2[0] = f;
          vec_2[1] = e;
          vec_2[2] = diag[2];
        }
      else if (max_diag == 1)
        {
          representative[0] = d;
          representative[1] = diag[1];
          representative[2] = e;
 
          vec_1[0] = f;
          vec_1[1] = e;
          vec_1[2] = diag[2];
 
          vec_2[0] = diag[0];
          vec_2[1] = d;
          vec_2[2] = f;
 
        }
      else /*if (max_diag == 2)*/
        {
          representative[0] = f;
          representative[1] = e;
          representative[2] = diag[2];
 
          vec_1[0] = diag[0];
          vec_1[1] = d;
          vec_1[2] = f;
 
          vec_2[0] = d;
          vec_2[1] = diag[1];
          vec_2[2] = e;
        }
 
      CaloRecGPU::Helpers::corrected_cross_product(res, representative, vec_1);
      CaloRecGPU::Helpers::corrected_cross_product(vec_1, representative, vec_2);
      //Can safely override previous value...
 
#ifdef __CUDA_ARCH__
      const float norm_1 = rnorm3df(res[0], res[1], res[2]);
      const float norm_2 = rnorm3df(vec_1[0], vec_1[1], vec_1[2]);
#else
      const float norm_1 = 1.f / hypot(res[0], res[1], res[2]);
      const float norm_2 = 1.f / hypot(vec_1[0], vec_1[1], vec_1[2]);
#endif
 
      if (norm_1 <= norm_2)
        //Greater magnitude -> multiply by a smaller value
        {
          res[0] *= norm_1;
          res[1] *= norm_1;
          res[2] *= norm_1;
        }
      else
        {
          res[0] = vec_1[0] * norm_2;
          res[1] = vec_1[1] * norm_2;
          res[2] = vec_1[2] * norm_2;
        }
    }

◆ get_eigenvalues()

CUDA_HOS_DEV void ClusterMomentsCalculator::RealSymmetricMatrixSolver::get_eigenvalues	(	float &	e_1,
		float &	e_2,
		float &	e_3 ) const

inline

Calculate shifted and scaled eigenvalues of the matrix, in ascending value.

To get the actual eigenvalues, you should multiply by scale and then add shift.

Definition at line 374 of file GPUClusterInfoAndMomentsCalculatorImpl.h.

    {
      using namespace std;
 
      const float ab      = a * b;
      const float corr_ab = fmaf(a, b, -ab);
      const float dd      = d * d;
      const float corr_dd = fmaf(d, d, -dd);
      const float ac      = a * c;
      const float corr_ac = fmaf(a, c, -ac);
      const float ff      = f * f;
      const float corr_ff = fmaf(f, f, -ff);
      const float bc      = b * c;
      const float corr_bc = fmaf(b, c, -bc);
      const float ee      = e * e;
      const float corr_ee = fmaf(e, e, -ee);
 
      const float c_0 = fmaf(2 * d, f * e,
                             CaloRecGPU::Helpers::product_sum_cornea_harrison_tang(ab,  c, -a, ee) +
                             CaloRecGPU::Helpers::product_sum_cornea_harrison_tang(-b, ff, -c, dd)   );
      //Consider using some other strategy to select the best pairs
      //to minimize error here?
 
      const float c_1 = CaloRecGPU::Helpers::sum_kahan_babushka_neumaier(ab, -dd, ac, -ff, bc, -ee, corr_ab, -corr_dd, corr_ac, -corr_ff, corr_bc, corr_ee);
 
      const float c_2 = CaloRecGPU::Helpers::sum_kahan_babushka_neumaier(a, b, c);
 
      constexpr float inv_3 = 1.f / 3.f;
 
      const float c_2_over_3 = c_2 * inv_3;
 
      const float a_over_3 = max(fma(c_2, c_2_over_3, -c_1) * inv_3, 0.f);
 
      const float half_b = 0.5f * (fma(c_2_over_3, fma(2.f * c_2_over_3, c_2_over_3, -c_1), c_0));
 
      const float q = max(CaloRecGPU::Helpers::product_sum_cornea_harrison_tang(a_over_3, a_over_3 * a_over_3, -half_b, half_b), 0.f);
 
      //None of what we are doing here is the best choice.
      //We would need to perform a proper, formal, numerical analysis
      //to figure out all possible errors and so on.
 
      const float rho = sqrtf(a_over_3);
 
#ifdef __CUDA_ARCH__
      const float theta = atan2f(1.0f, rsqrtf(q) * half_b) * inv_3;
#else
      const float theta = atan2f(sqrtf(q), half_b) * inv_3;
#endif
 
#ifdef __CUDA_ARCH__
      float sin_theta, cos_theta;
      sincosf(theta, &sin_theta, &cos_theta);
#else
      const float sin_theta = sinf(theta);
      const float cos_theta = cosf(theta);
#endif
 
      const float sqrt_3 = sqrtf(3.f);
 
      e_1 = fma(-rho, fma( sqrt_3, sin_theta, cos_theta), c_2_over_3);
      e_2 = fma(-rho, fma(-sqrt_3, sin_theta, cos_theta), c_2_over_3);
      e_3 = fma( rho, 2.0f * cos_theta,                   c_2_over_3);
    }

◆ get_eigenvectors()

CUDA_HOS_DEV void ClusterMomentsCalculator::RealSymmetricMatrixSolver::get_eigenvectors	(	float(&)	res[3][3],
		const float	e_1,
		const float	e_2,
		const float	e_3,
		const float	epsilon = s_typical_epsilon ) const

inline

Calculate the eigenvectors of the matrix, using the (possibly unscaled) eigenvalues e_1, e_2, e_3 (in ascending order of magnitude) and epsilon to guard against undefined cases.

Definition at line 538 of file GPUClusterInfoAndMomentsCalculatorImpl.h.

    {
      using namespace std;
 
      if (e_3 - e_1 <= epsilon)
        {
          res[0][0] = 1.f;
          res[0][1] = 0.f;
          res[0][2] = 0.f;
 
          res[1][0] = 0.f;
          res[1][1] = 1.f;
          res[1][2] = 0.f;
 
          res[2][0] = 0.f;
          res[2][1] = 0.f;
          res[2][2] = 1.f;
        }
      else
        {
          const float d_0 = e_3 - e_2;
          const float d_1 = e_2 - e_1;
 
          const float d_min = min(d_0, d_1);
          const float d_max = max(d_0, d_1);
 
          int k, j;
          float first_e, second_e;
 
          if (d_0 > d_1)
            {
              k = 2;
              j = 0;
 
              first_e  = e_3;
              second_e = e_1;
            }
          else
            {
              k = 0;
              j = 2;
 
              first_e  = e_1;
              second_e = e_3;
            }
 
          extract_one(first_e, res[k], res[j]);
#if USE_ORIGINAL_EIGEN
          if (d_min <= 2 * epsilon * d1)
            {
#ifdef __CUDA_ARCH__
              const float base_norm = rnorm3df(res[j][0], res[j][1], res[j][2]);
#else
              const float base_norm = 1.f / hypot(res[j][0], res[j][1], res[j][2]);
#endif
              const float extra_factor = 1.f - CaloRecGPU::Helpers::corrected_dot_product(res[k], res[j]);
 
              const float norm = base_norm / extra_factor;
 
              res[j][0] *= norm;
              res[j][1] *= norm;
              res[j][2] *= norm;
            }
#else
          if (d_min <= 2 * epsilon * d_max)
            //Eigen has d0 <= 2 * eps <= d1, but d0 is overwritten while d1 isn't...
            {
              //Eigen speaks about ortho-normalization
              //but does eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l)
              //which... does not ortho-normalize anything.
 
              const float prod = CaloRecGPU::Helpers::corrected_dot_product(res[k], res[j]);
 
              res[j][0] -= res[k][0] * prod;
              res[j][1] -= res[k][1] * prod;
              res[j][2] -= res[k][2] * prod;
 
#ifdef __CUDA_ARCH__
              const float norm = rnorm3df(res[j][0], res[j][1], res[j][2]);
#else
              const float norm = 1.f / hypot(res[j][0], res[j][1], res[j][2]);
#endif
              res[j][0] *= norm;
              res[j][1] *= norm;
              res[j][2] *= norm;
            }
#endif
          else
            {
              float extra_vector[3];
 
              extract_one(second_e, res[j], extra_vector);
            }
 
          CaloRecGPU::Helpers::corrected_cross_product(res[1], res[2], res[0]);
 
#ifdef __CUDA_ARCH__
          const float norm = rnorm3df(res[1][0], res[1][1], res[1][2]);
#else
          const float norm = 1.f / hypot(res[1][0], res[1][1], res[1][2]);
#endif
 
          res[1][0] *= norm;
          res[1][1] *= norm;
          res[1][2] *= norm;
        }
    }

◆ get_solution()

CUDA_HOS_DEV void ClusterMomentsCalculator::RealSymmetricMatrixSolver::get_solution	(	float(&)	eigenvalues[3],
		float(&)	eigenvectors[3][3],
		bool	rescale_and_reshift_values = true,
		const float	epsilon = s_typical_epsilon )

inline

Get the full eigenvalues and eigenvectors for this matrix.

If rescale_and_reshift_values is true, the eigenvalues are scaled and shifted back to their proper value, given the original matrix.

Definition at line 650 of file GPUClusterInfoAndMomentsCalculatorImpl.h.

    {
      get_eigenvalues(eigenvalues[0], eigenvalues[1], eigenvalues[2]);
      get_eigenvectors(eigenvectors, eigenvalues[0], eigenvalues[1], eigenvalues[2], epsilon);
 
      if (rescale_and_reshift_values)
        {
          eigenvalues[0] = fmaf(eigenvalues[0], scale, shift);
          eigenvalues[1] = fmaf(eigenvalues[1], scale, shift);
          eigenvalues[2] = fmaf(eigenvalues[2], scale, shift);
        }
    }