FLA_Tevd_francis_v_opt_var1.c File Reference

(r)

Functions
FLA_Error	FLA_Tevd_francis_v_opt_var1 (FLA_Obj shift, FLA_Obj g, FLA_Obj d, FLA_Obj e)
FLA_Error	FLA_Tevd_francis_v_ops_var1 (int m_A, float *buff_shift, scomplex *buff_g, int inc_g, float *buff_d, int inc_d, float *buff_e, int inc_e)
FLA_Error	FLA_Tevd_francis_v_opd_var1 (int m_A, double *buff_shift, dcomplex *buff_g, int inc_g, double *buff_d, int inc_d, double *buff_e, int inc_e)

Function Documentation

FLA_Tevd_francis_v_opd_var1()

FLA_Error FLA_Tevd_francis_v_opd_var1	(	int	m_A,
		double *	buff_shift,
		dcomplex *	buff_g,
		int	inc_g,
		double *	buff_d,
		int	inc_d,
		double *	buff_e,
		int	inc_e
	)

{
    double    eps, safmin;
    double    temp0, temp1;
    double    bulge;
    int       ij_deflated;
    int       i;
 
    // Initialize the deflation index.
    ij_deflated = FLA_SUCCESS;
 
    // Initialize the bulge variable to zero.
    bulge = 0.0;
 
    // Query epsilon and safmin.
    eps    = FLA_Mach_params_opd( FLA_MACH_EPS );
    safmin = FLA_Mach_params_opd( FLA_MACH_SFMIN );
 
    // Apply the rotations in forward order.
    for ( i = 0; i < m_A - 1; ++i )
    {
        double*   alpha00  = buff_d + (i-1)*inc_d;
        double*   alpha10  = buff_e + (i-1)*inc_e;
        double*   alpha20  = &bulge;
 
        double*   alpha11  = buff_d + (i  )*inc_d;
        double*   alpha21  = buff_e + (i  )*inc_e;
        double*   alpha22  = buff_d + (i+1)*inc_d;
 
        double*   alpha31  = &bulge;
        double*   alpha32  = buff_e + (i+1)*inc_e;
 
        double*   gamma1   = &(buff_g[(i  )*inc_g].real);
        double*   sigma1   = &(buff_g[(i  )*inc_g].imag);
 
        double    alpha10_new;
 
        int       m_behind = i;
        int       m_ahead  = m_A - i - 2;
 
        /*------------------------------------------------------------*/
 
        if ( i == 0 )
        {
            // Induce an iteration that introduces the bulge by
            // changing the addresses of alpha10 and alpha20.
            temp0 = *buff_d - *buff_shift;
            temp1 = *buff_e;
            alpha10 = &temp0;
            alpha20 = &temp1;
 
            // Compute a new Givens rotation that introduces the bulge.
            MAC_Givens2_opd( alpha10,
            //FLA_Givens2_opd( alpha10,
                             alpha20,
                             gamma1,
                             sigma1,
                             &alpha10_new );
 
            // We don't apply the Givens rotation to the 2x1 vector at
            // alpha10 when introducing the bulge.
        }
        else
        {
            // Compute a new Givens rotation to push the bulge.
            MAC_Givens2_opd( alpha10,
            //FLA_Givens2_opd( alpha10,
                             alpha20,
                             gamma1,
                             sigma1,
                             &alpha10_new );
 
            // Apply the Givens rotation to the 2x1 vector from which it
            // was computed, which annihilates alpha20.
            *alpha10 = alpha10_new;
            *alpha20 = 0.0;
        }
 
        // Apply the Givens rotation to the 2x2 submatrix at alpha11.
        MAC_Apply_GTG_opd( gamma1,
                           sigma1,
                           alpha11,
                           alpha21,
                           alpha22 );
 
        if ( m_ahead > 0 )
        {
            // Apply the Givens rotation to the 1x2 vector below the 2x2
            // submatrix. This should move the bulge to alpha31.
            MAC_Apply_G_1x2_opd( gamma1,
                                 sigma1,
                                 alpha31,
                                 alpha32 );
 
            // Check for deflation after applying the rotations, except after
            // applying the initial bulge-introducing rotations.
            if ( m_behind > 0 )
            {
                // We check for deflation in the previous column now that we
                // are done modifying it. If deflation occurred, record the
                // index.
                if ( MAC_Tevd_eigval_converged_opd( eps, safmin, *alpha00, *alpha10, *alpha11 ) )
                {
                    ij_deflated = i - 1;
                }
            }
 
            // Sanity check. If the bulge ever disappears, something may be wrong.
            if ( *alpha31 == 0.0 )
            {
                printf( "FLA_Tevd_francis_v_opt_var1: bulge disappeared!\n" );
                if ( MAC_Tevd_eigval_converged_opd( eps, safmin, *alpha11, *alpha21, *alpha22 ) )
                {
                    printf( "FLA_Tevd_francis_v_opt_var1: deflation detected (col %d)\n", i );
                    printf( "FLA_Tevd_francis_v_opt_var1: alpha11         = %23.19e\n", *alpha11 );
                    printf( "FLA_Tevd_francis_v_opt_var1: alpha21 alpha22 = %23.19e %23.19e\n", *alpha21, *alpha22 );
                    return i;
                }
                else
                {
                    printf( "FLA_Tevd_francis_v_opt_var1: but NO deflation detected! (col %d)\n", i );
                    printf( "FLA_Tevd_francis_v_opt_var1: alpha11         = %23.19e\n", *alpha11 );
                    printf( "FLA_Tevd_francis_v_opt_var1: alpha21 alpha22 = %23.19e %23.19e\n", *alpha21, *alpha22 );
                    FLA_Abort();
                    return FLA_FAILURE;
                }
            }
        }
 
        /*------------------------------------------------------------*/
    }
 
    // Return the index of column where deflation most recently occurred,
    // or FLA_SUCCESS if no deflation was detected.
    return ij_deflated;
}

References FLA_Abort(), FLA_Mach_params_opd(), i, imag, and temp1.

Referenced by FLA_Tevd_eigval_v_opd_var1(), FLA_Tevd_eigval_v_opd_var3(), and FLA_Tevd_francis_v_opt_var1().

FLA_Tevd_francis_v_ops_var1()

FLA_Error FLA_Tevd_francis_v_ops_var1	(	int	m_A,
		float *	buff_shift,
		scomplex *	buff_g,
		int	inc_g,
		float *	buff_d,
		int	inc_d,
		float *	buff_e,
		int	inc_e
	)

{
    return FLA_SUCCESS;
}

Referenced by FLA_Tevd_francis_v_opt_var1().

FLA_Tevd_francis_v_opt_var1()

FLA_Error FLA_Tevd_francis_v_opt_var1	(	FLA_Obj	shift,
		FLA_Obj	g,
		FLA_Obj	d,
		FLA_Obj	e
	)

{
    FLA_Datatype datatype;
    int          m_A;
    int          inc_g;
    int          inc_d;
    int          inc_e;
 
    datatype = FLA_Obj_datatype( d );
 
    m_A      = FLA_Obj_vector_dim( d );
 
    inc_g    = FLA_Obj_vector_inc( g );
    inc_d    = FLA_Obj_vector_inc( d );
    inc_e    = FLA_Obj_vector_inc( e );
    
 
    switch ( datatype )
    {
        case FLA_FLOAT:
        {
            float*    buff_shift = FLA_FLOAT_PTR( shift );
            scomplex* buff_g     = FLA_COMPLEX_PTR( g );
            float*    buff_d     = FLA_FLOAT_PTR( d );
            float*    buff_e     = FLA_FLOAT_PTR( e );
 
            FLA_Tevd_francis_v_ops_var1( m_A,
                                         buff_shift,
                                         buff_g, inc_g,
                                         buff_d, inc_d,
                                         buff_e, inc_e );
 
            break;
        }
 
        case FLA_DOUBLE:
        {
            double*   buff_shift = FLA_DOUBLE_PTR( shift );
            dcomplex* buff_g     = FLA_DOUBLE_COMPLEX_PTR( g );
            double*   buff_d     = FLA_DOUBLE_PTR( d );
            double*   buff_e     = FLA_DOUBLE_PTR( e );
 
            FLA_Tevd_francis_v_opd_var1( m_A,
                                         buff_shift,
                                         buff_g, inc_g,
                                         buff_d, inc_d,
                                         buff_e, inc_e );
 
            break;
        }
    }
 
    return FLA_SUCCESS;
}

References FLA_Obj_datatype(), FLA_Obj_vector_dim(), FLA_Obj_vector_inc(), FLA_Tevd_francis_v_opd_var1(), and FLA_Tevd_francis_v_ops_var1().