FLA_Wilkshift_tridiag.c File Reference

(r12600)

Functions
FLA_Error	FLA_Wilkshift_tridiag (FLA_Obj delta1, FLA_Obj epsilon, FLA_Obj delta2, FLA_Obj kappa)
FLA_Error	FLA_Wilkshift_tridiag_ops (float delta1, float epsilon, float delta2, float *kappa)
double	fla_dlapy2 (double x, double y)
FLA_Error	FLA_Wilkshift_tridiag_opd (double delta1, double epsilon, double delta2, double *kappa)

---

Function Documentation

double fla_dlapy2	(	double	x,
		double	y
	)

FLA_Error FLA_Wilkshift_tridiag	(	FLA_Obj	delta1,
		FLA_Obj	epsilon,
		FLA_Obj	delta2,
		FLA_Obj	kappa
	)

References FLA_Check_error_level(), FLA_Obj_datatype(), FLA_Wilkshift_tridiag_check(), FLA_Wilkshift_tridiag_opd(), and FLA_Wilkshift_tridiag_ops().

{
    FLA_Datatype datatype;

    datatype = FLA_Obj_datatype( delta1 );

    if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
        FLA_Wilkshift_tridiag_check( delta1, epsilon, delta2, kappa );

    switch ( datatype )
    {
        case FLA_FLOAT:
        {
            float* delta1_p  = ( float* ) FLA_FLOAT_PTR( delta1 );
            float* epsilon_p = ( float* ) FLA_FLOAT_PTR( epsilon );
            float* delta2_p  = ( float* ) FLA_FLOAT_PTR( delta2 );
            float* kappa_p   = ( float* ) FLA_FLOAT_PTR( kappa );

            FLA_Wilkshift_tridiag_ops( *delta1_p,
                                       *epsilon_p,
                                       *delta2_p,
                                       kappa_p );

            break;
        }

        case FLA_DOUBLE:
        {
            double* delta1_p  = ( double* ) FLA_DOUBLE_PTR( delta1 );
            double* epsilon_p = ( double* ) FLA_DOUBLE_PTR( epsilon );
            double* delta2_p  = ( double* ) FLA_DOUBLE_PTR( delta2 );
            double* kappa_p   = ( double* ) FLA_DOUBLE_PTR( kappa );

            FLA_Wilkshift_tridiag_opd( *delta1_p,
                                       *epsilon_p,
                                       *delta2_p,
                                       kappa_p );

            break;
        }
    }

    return FLA_SUCCESS;
}

FLA_Error FLA_Wilkshift_tridiag_opd	(	double	delta1,
		double	epsilon,
		double	delta2,
		double *	kappa
	)

Referenced by FLA_Tevd_eigval_n_opd_var1(), FLA_Tevd_eigval_v_opd_var1(), FLA_Tevd_eigval_v_opd_var3(), FLA_Tevd_find_perfshift_opd(), and FLA_Wilkshift_tridiag().

{
    double a = delta1;
    double c = epsilon;
    double d = delta2;
    double p, q, r, s, k;

    // Begin with kappa equal to d.
    k = d;

    // Compute a scaling factor to promote numerical stability.
    s = fabs( a ) + 2.0 * fabs( c ) + fabs( d );

    if ( s == 0.0 ) return FLA_SUCCESS;

    // Compute q with scaling applied.
    q = ( c / s ) * ( c / s );

    if ( q != 0.0 )
    {

        // Compute p = 0.5*( a - d ) with scaling applied.
        p = 0.5 * ( ( a / s ) - ( d / s ) );

        // Compute r = sqrt( p*p - q ).
        r = sqrt( p * p + q );

        // If p*r is negative, then we need to negate r to ensure we use the
        // larger of the two eigenvalues.
        if ( p * r < 0.0 ) r = -r;

        // Compute the Wilkinson shift with scaling removed:
        //   k = lambda_min + d
        //     = d + lambda_min
        //     = d + (-q / lambda_max)
        //     = d - q / lambda_max
        //     = d - q / (p + r)
        k = k - s * ( q / ( p + r ) );

/*
        // LAPACK method:

        p = 0.5 * ( ( a ) - ( d ) ) / c ;
        //r = sqrt( p * p + 1.0 );
        r = fla_dlapy2( p, 1.0 );
        if ( p < 0.0 ) r = -r;
        k = k - ( c / ( p + r ) );
*/
    }

    // Save the result.
    *kappa = k;

    return FLA_SUCCESS;
}

FLA_Error FLA_Wilkshift_tridiag_ops	(	float	delta1,
		float	epsilon,
		float	delta2,
		float *	kappa
	)

Referenced by FLA_Wilkshift_tridiag().

{
    float  a = delta1;
    float  c = epsilon;
    float  d = delta2;
    float  p, q, r, s, k;

    // Begin with kappa equal to d.
    k = d;

    // Compute a scaling factor to promote numerical stability.
    s = fabs( a ) + 2.0F * fabs( c ) + fabs( d );

    if ( s == 0.0F ) return FLA_SUCCESS;

    // Compute q with scaling applied.
    q = ( c / s ) * ( c / s );

    if ( q != 0.0F )
    {
        // Compute p = 0.5*( a - d ) with scaling applied.
        p = 0.5F * ( ( a / s ) - ( d / s ) );

        // Compute r = sqrt( p*p - q ).
        r = sqrt( p * p + q );

        // If p*r is negative, then we need to negate r to ensure we use the
        // larger of the two eigenvalues.
        if ( p * r < 0.0F ) r = -r;

        // Compute the Wilkinson shift with scaling removed:
        //   k = lambda_min + d
        //     = d + lambda_min
        //     = d + (-q / lambda_max)
        //     = d - q / lambda_max
        //     = d - q / (p + r)
        k = k - s * ( q / ( p + r ) );
    }

    // Save the result.
    *kappa = k;

    return FLA_SUCCESS;
}