#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. More...

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. More...

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. More...

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 445 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 455 of file GPUClusterInfoAndMomentsCalculatorImpl.h.

     {
       using namespace std;
  
       shift = 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 547 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;
         }
  
       corrected_cross_product(res, representative, vec_1);
       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 483 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,
                              product_sum_cornea_harrison_tang(ab,  c, -a, ee) +
                              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 = 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 = 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(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 647 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 - 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 = 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);
             }
  
           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 759 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);
         }
     }