compactness.h

Go to the documentation of this file.

 
#ifndef VA_COMPACTNESS_H
#define VA_COMPACTNESS_H
 
#include "common_defs.h"
#include "linalg3.h"
 
/***************************************************************************/
VA_DEVICE_FUN void valib_compactness(VA_REAL *A, VA_REAL *compactness)
{
    // cast the eigenvalue problem to finding roots of the cubic characteristic
    // equation
    // ax^3 + bx^2 + cx + d = 0
    //  x^3 - trace(A) x^2 + II_A x - det(A) = 0
    VA_REAL a, b, c, d;
 
    a = 1.;
 
    valib_trace3(A, &b);
    b = -b;
 
    valib_second_invariant3(A, &c);
 
    valib_determinant3(A, &d);
    d = -d;
 
    // find the three complex roots
    int info;
    VA_COMPLEX x1, x2, x3;
    valib_solve_cubic_equation(a, b, c, d, &x1, &x2, &x3, &info);
 
    VA_REAL lambda_r, lambda_cr, lambda_ci;
 
    if (info == -1) {
        // only real eigenvalues exist, compactness cannot be determined
        *compactness = -VA_INFINITY;
        return;
    }
    else if (info == 0) {
        // search for eigenvalues was successful
        lambda_r  = creal(x1);
        lambda_cr = creal(x2);
        lambda_ci = fabs(cimag(x2));
    }
    else {
        // there was another error in the search for eigenvalues
        printf("The solve of the cubic equation has failed with info: %d \n", info);
        *compactness = -VA_INFINITY;
        return;
    }
 
    // now find the compactness criterion
    //*compactness = 1. - 1./3. * pow (lambda_cr / lambda_ci - lambda_r / lambda_ci, 2.);
    *compactness = 1. - 1./3. * valib_square(lambda_cr / lambda_ci - lambda_r / lambda_ci) 
                 - valib_sign(1.,2.*lambda_cr+lambda_r)*valib_square(2*lambda_cr+lambda_r) / 6. / valib_square(lambda_ci);
 
    return;
}
 
#endif // VA_COMPACTNESS_H