|
|
|
@ -8,6 +8,7 @@
|
|
|
|
|
|
|
|
|
|
#include <float.h> // FLT_MAX
|
|
|
|
|
#include <vector>
|
|
|
|
|
#include <string.h>
|
|
|
|
|
|
|
|
|
|
using namespace nv;
|
|
|
|
|
|
|
|
|
@ -498,7 +499,7 @@ static void EigenSolver3_Tridiagonal(float mat[3][3], float * diag, float * subd
|
|
|
|
|
|
|
|
|
|
diag[0] = a;
|
|
|
|
|
subd[2] = 0.f;
|
|
|
|
|
if ( fabs(c) >= epsilon )
|
|
|
|
|
if (fabsf(c) >= epsilon)
|
|
|
|
|
{
|
|
|
|
|
const float ell = sqrtf(b*b+c*c);
|
|
|
|
|
b /= ell;
|
|
|
|
@ -538,8 +539,8 @@ static bool EigenSolver3_QLAlgorithm(float mat[3][3], float * diag, float * subd
|
|
|
|
|
int m;
|
|
|
|
|
for (m = ell; m <= 1; m++)
|
|
|
|
|
{
|
|
|
|
|
float dd = fabs(diag[m]) + fabs(diag[m+1]);
|
|
|
|
|
if ( fabs(subd[m]) + dd == dd )
|
|
|
|
|
float dd = fabsf(diag[m]) + fabsf(diag[m+1]);
|
|
|
|
|
if ( fabsf(subd[m]) + dd == dd )
|
|
|
|
|
break;
|
|
|
|
|
}
|
|
|
|
|
if ( m == ell )
|
|
|
|
@ -555,7 +556,7 @@ static bool EigenSolver3_QLAlgorithm(float mat[3][3], float * diag, float * subd
|
|
|
|
|
for (int i = m-1; i >= ell; i--)
|
|
|
|
|
{
|
|
|
|
|
float f = s*subd[i], b = c*subd[i];
|
|
|
|
|
if ( fabs(f) >= fabs(g) )
|
|
|
|
|
if ( fabsf(f) >= fabsf(g) )
|
|
|
|
|
{
|
|
|
|
|
c = g/f;
|
|
|
|
|
r = sqrtf(c*c+1);
|
|
|
|
@ -690,7 +691,7 @@ static void EigenSolver4_Tridiagonal(float mat[4][4], float * diag, float * subd
|
|
|
|
|
float maxElement = FLT_MAX;
|
|
|
|
|
for (int i = 0; i < n; ++i)
|
|
|
|
|
for (int j = 0; j < n; ++j)
|
|
|
|
|
maxElement = max(maxElement, fabs(mat[i][j]));
|
|
|
|
|
maxElement = max(maxElement, fabsf(mat[i][j]));
|
|
|
|
|
float epsilon = relEpsilon * maxElement;
|
|
|
|
|
|
|
|
|
|
// Iterative algorithm, works for any size of matrix but might be slower than
|
|
|
|
@ -711,7 +712,7 @@ static void EigenSolver4_Tridiagonal(float mat[4][4], float * diag, float * subd
|
|
|
|
|
float r = sqrtf(0.5f * (alpha*alpha - A(k+1,k)*alpha));
|
|
|
|
|
|
|
|
|
|
// If r is zero, skip this column - already in tridiagonal form
|
|
|
|
|
if (fabs(r) < epsilon)
|
|
|
|
|
if (fabsf(r) < epsilon)
|
|
|
|
|
continue;
|
|
|
|
|
|
|
|
|
|
float v[n] = {};
|
|
|
|
@ -728,12 +729,12 @@ static void EigenSolver4_Tridiagonal(float mat[4][4], float * diag, float * subd
|
|
|
|
|
Q = mul(Q, P);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
nvDebugCheck(fabs(A(2,0)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabs(A(0,2)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabs(A(3,0)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabs(A(0,3)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabs(A(3,1)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabs(A(1,3)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabsf(A(2,0)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabsf(A(0,2)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabsf(A(3,0)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabsf(A(0,3)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabsf(A(3,1)) < epsilon);
|
|
|
|
|
nvDebugCheck(fabsf(A(1,3)) < epsilon);
|
|
|
|
|
|
|
|
|
|
for (int i = 0; i < n; ++i)
|
|
|
|
|
diag[i] = A(i,i);
|
|
|
|
@ -758,8 +759,8 @@ static bool EigenSolver4_QLAlgorithm(float mat[4][4], float * diag, float * subd
|
|
|
|
|
int m;
|
|
|
|
|
for (m = ell; m < 3; m++)
|
|
|
|
|
{
|
|
|
|
|
float dd = fabs(diag[m]) + fabs(diag[m+1]);
|
|
|
|
|
if ( fabs(subd[m]) + dd == dd )
|
|
|
|
|
float dd = fabsf(diag[m]) + fabsf(diag[m+1]);
|
|
|
|
|
if ( fabsf(subd[m]) + dd == dd )
|
|
|
|
|
break;
|
|
|
|
|
}
|
|
|
|
|
if ( m == ell )
|
|
|
|
@ -775,7 +776,7 @@ static bool EigenSolver4_QLAlgorithm(float mat[4][4], float * diag, float * subd
|
|
|
|
|
for (int i = m-1; i >= ell; i--)
|
|
|
|
|
{
|
|
|
|
|
float f = s*subd[i], b = c*subd[i];
|
|
|
|
|
if ( fabs(f) >= fabs(g) )
|
|
|
|
|
if ( fabsf(f) >= fabsf(g) )
|
|
|
|
|
{
|
|
|
|
|
c = g/f;
|
|
|
|
|
r = sqrtf(c*c+1);
|
|
|
|
|