libflame  revision_anchor
Functions
FLA_Bsvd_compute_shift.c File Reference

(r)

Functions

FLA_Error FLA_Bsvd_compute_shift (FLA_Obj tol, FLA_Obj sminl, FLA_Obj smax, FLA_Obj d, FLA_Obj e, FLA_Obj shift)
FLA_Error FLA_Bsvd_compute_shift_ops (int m_A, float tol, float sminl, float smax, float *buff_d, int inc_d, float *buff_e, int inc_e, float *shift)
FLA_Error FLA_Bsvd_compute_shift_opd (int m_A, double tol, double sminl, double smax, double *buff_d, int inc_d, double *buff_e, int inc_e, double *shift)

Function Documentation

FLA_Error FLA_Bsvd_compute_shift ( FLA_Obj  tol,
FLA_Obj  sminl,
FLA_Obj  smax,
FLA_Obj  d,
FLA_Obj  e,
FLA_Obj  shift 
)

References FLA_Bsvd_compute_shift_opd(), FLA_Bsvd_compute_shift_ops(), FLA_Obj_datatype(), FLA_Obj_vector_dim(), and FLA_Obj_vector_inc().

{
    FLA_Datatype datatype;
    int          m_A;
    int          inc_d;
    int          inc_e;

    datatype = FLA_Obj_datatype( d );

    m_A      = FLA_Obj_vector_dim( d );

    inc_d    = FLA_Obj_vector_inc( d );
    inc_e    = FLA_Obj_vector_inc( e );

    switch ( datatype )
    {
    case FLA_FLOAT:
    {
        float*    buff_tol       = FLA_FLOAT_PTR( tol );
        float*    buff_sminl     = FLA_FLOAT_PTR( sminl );
        float*    buff_smax      = FLA_FLOAT_PTR( smax );
        float*    buff_d         = FLA_FLOAT_PTR( d );
        float*    buff_e         = FLA_FLOAT_PTR( e );
        float*    buff_shift     = FLA_FLOAT_PTR( shift );

        FLA_Bsvd_compute_shift_ops( m_A,
                                    *buff_tol,
                                    *buff_sminl,
                                    *buff_smax,
                                    buff_d, inc_d,
                                    buff_e, inc_e,
                                    buff_shift );

        break;
    }

    case FLA_DOUBLE:
    {
        double*   buff_tol       = FLA_DOUBLE_PTR( tol );
        double*   buff_sminl     = FLA_DOUBLE_PTR( sminl );
        double*   buff_smax      = FLA_DOUBLE_PTR( smax );
        double*   buff_d         = FLA_DOUBLE_PTR( d );
        double*   buff_e         = FLA_DOUBLE_PTR( e );
        double*   buff_shift     = FLA_DOUBLE_PTR( shift );

        FLA_Bsvd_compute_shift_opd( m_A,
                                    *buff_tol,
                                    *buff_sminl,
                                    *buff_smax,
                                    buff_d, inc_d,
                                    buff_e, inc_e,
                                    buff_shift );

        break;
    }
    }

    return FLA_SUCCESS;
}
FLA_Error FLA_Bsvd_compute_shift_opd ( int  m_A,
double  tol,
double  sminl,
double  smax,
double *  buff_d,
int  inc_d,
double *  buff_e,
int  inc_e,
double *  shift 
)

References FLA_Mach_params_opd(), and FLA_Sv_2x2_opd().

Referenced by FLA_Bsvd_compute_shift(), and FLA_Bsvd_sinval_v_opd_var1().

{
    double  hndrth = 0.01;
    double  eps;
    double* d_first;
    double* e_last;
    double* d_last_m1;
    double* d_last;
    double  sll, temp;

    eps = FLA_Mach_params_opd( FLA_MACH_EPS );

    d_first   = buff_d + (0    )*inc_d;
    e_last    = buff_e + (m_A-2)*inc_e;
    d_last_m1 = buff_d + (m_A-2)*inc_d;
    d_last    = buff_d + (m_A-1)*inc_d;

    // If the shift would ruin relative accuracy, set it to zero.
    if ( m_A * tol * ( sminl / smax ) <= max( eps, hndrth * tol ) )
    {
#ifdef PRINTF
        printf( "FLA_Bsvd_compute_shift_opd: shift would ruin accuracy; setting shift to 0.\n" );
        printf( "                   m_A = %d     \n", m_A );
        printf( "                   tol = %20.15e\n", tol );
        printf( "                 sminl = %20.15e\n", sminl );
        printf( "                  smax = %20.15e\n", smax );
        printf( "                   LHS = %20.15e\n", m_A * tol * ( sminl / smax ) );
        printf( "      max(eps,0.01*tol)= %20.15e\n", max( eps, hndrth * tol ) );
#endif
        *shift = 0.0;
    }
    else
    {
        // Compute the shift from the last 2x2 matrix.
        FLA_Sv_2x2_opd( d_last_m1,
                        e_last,
                        d_last,
                        shift,
                        &temp );

        sll = fabs( *d_first );

        // Check if the shift is negligible; if so, set it to zero.
        if ( sll > 0.0 )
        {
            temp = ( *shift / sll );
            if ( temp * temp < eps )
            {
#ifdef PRINTF
                printf( "FLA_Bsvd_compute_shift_opd: shift is negligible; setting shift to 0.\n" );
#endif
                *shift = 0.0;
            }
        }
    }

    return FLA_SUCCESS;
}
FLA_Error FLA_Bsvd_compute_shift_ops ( int  m_A,
float  tol,
float  sminl,
float  smax,
float *  buff_d,
int  inc_d,
float *  buff_e,
int  inc_e,
float *  shift 
)

References FLA_Mach_params_ops(), and FLA_Sv_2x2_ops().

Referenced by FLA_Bsvd_compute_shift(), and FLA_Bsvd_sinval_v_ops_var1().

{
    float  hndrth = 0.01;
    float  eps;
    float* d_first;
    float* e_last;
    float* d_last_m1;
    float* d_last;
    float  sll, temp;

    eps = FLA_Mach_params_ops( FLA_MACH_EPS );

    d_first   = buff_d + (0    )*inc_d;
    e_last    = buff_e + (m_A-2)*inc_e;
    d_last_m1 = buff_d + (m_A-2)*inc_d;
    d_last    = buff_d + (m_A-1)*inc_d;

    // If the shift would ruin relative accuracy, set it to zero.
    if ( m_A * tol * ( sminl / smax ) <= max( eps, hndrth * tol ) )
    {
        *shift = 0.0;
    }
    else
    {
        // Compute the shift from the last 2x2 matrix.
        FLA_Sv_2x2_ops( d_last_m1,
                        e_last,
                        d_last,
                        shift,
                        &temp );

        sll = fabsf( *d_first );

        // Check if the shift is negligible; if so, set it to zero.
        if ( sll > 0.0F )
        {
            temp = ( *shift / sll );
            if ( temp * temp < eps )
            {
                *shift = 0.0F;
            }
        }
    }

    return FLA_SUCCESS;
}