|
libflame
revision_anchor
|
Go to the source code of this file.
Functions | |
| FLA_Error | FLA_Bsvd_francis_v_opt_var1 (FLA_Obj shift, FLA_Obj g, FLA_Obj h, FLA_Obj d, FLA_Obj e) |
| FLA_Error | FLA_Bsvd_francis_v_ops_var1 (int m_A, float shift, scomplex *buff_g, int inc_g, scomplex *buff_h, int inc_h, float *buff_d, int inc_d, float *buff_e, int inc_e) |
| FLA_Error | FLA_Bsvd_francis_v_opd_var1 (int m_A, double shift, dcomplex *buff_g, int inc_g, dcomplex *buff_h, int inc_h, double *buff_d, int inc_d, double *buff_e, int inc_e) |
| FLA_Error FLA_Bsvd_francis_v_opd_var1 | ( | int | m_A, |
| double | shift, | ||
| dcomplex * | buff_g, | ||
| int | inc_g, | ||
| dcomplex * | buff_h, | ||
| int | inc_h, | ||
| double * | buff_d, | ||
| int | inc_d, | ||
| double * | buff_e, | ||
| int | inc_e | ||
| ) |
References bli_d1().
Referenced by FLA_Bsvd_francis_v_opt_var1(), and FLA_Bsvd_sinval_v_opd_var1().
{
double one = bli_d1();
double bulge = 0.0;
int i;
// If the shift is zero, perform a simplified Francis step.
if ( shift == 0.0 )
{
double cs, oldcs;
double sn, oldsn = 0;
double r, temp;
double a11temp, a22temp;
double* d_last = buff_d + (m_A-1)*inc_d;
double* e_last = buff_e + (m_A-2)*inc_e;
cs = one;
oldcs = one;
for ( i = 0; i < m_A - 1; ++i )
{
double* alpha01 = buff_e + (i-1)*inc_e;
double* alpha11 = buff_d + (i )*inc_d;
double* alpha12 = buff_e + (i )*inc_e;
double* alpha22 = buff_d + (i+1)*inc_d;
double* gammaL = &(buff_g[(i )*inc_g].real);
double* sigmaL = &(buff_g[(i )*inc_g].imag);
double* gammaR = &(buff_h[(i )*inc_h].real);
double* sigmaR = &(buff_h[(i )*inc_h].imag);
a11temp = *alpha11 * cs;
MAC_Givens2_opd( &a11temp,
alpha12,
&cs,
&sn,
&r );
if ( i > 0 ) *alpha01 = oldsn * r;
a11temp = oldcs * r;
a22temp = *alpha22 * sn;
MAC_Givens2_opd( &a11temp,
&a22temp,
&oldcs,
&oldsn,
alpha11 );
*gammaR = cs;
*sigmaR = sn;
*gammaL = oldcs;
*sigmaL = oldsn;
}
temp = *d_last * cs;
*d_last = temp * oldcs;
*e_last = temp * oldsn;
return FLA_SUCCESS;
}
// Apply Givens rotations in forward order.
for ( i = 0; i < m_A - 1; ++i )
{
double* alpha01 = buff_e + (i-1)*inc_e;
double* alpha11 = buff_d + (i )*inc_d;
double* alpha21 = &bulge;
double* alpha02 = &bulge;
double* alpha12 = buff_e + (i )*inc_e;
double* alpha22 = buff_d + (i+1)*inc_d;
double* alpha13 = &bulge;
double* alpha23 = buff_e + (i+1)*inc_e;
double* gammaL = &(buff_g[(i )*inc_g].real);
double* sigmaL = &(buff_g[(i )*inc_g].imag);
double* gammaR = &(buff_h[(i )*inc_h].real);
double* sigmaR = &(buff_h[(i )*inc_h].imag);
double alpha01_new;
double alpha11_new;
int mn_ahead = m_A - i - 2;
/*------------------------------------------------------------*/
if ( i == 0 )
{
double alpha11_temp;
double alpha12_temp;
// Induce an iteration that introduces the bulge (from the right).
//alpha11_temp = *buff_d - shift;
alpha11_temp = ( fabs( *buff_d ) - shift ) *
( signof( one, *buff_d ) + ( shift / *buff_d ) );
alpha12_temp = *buff_e;
// Compute a Givens rotation that introduces the bulge.
MAC_Givens2_opd( &alpha11_temp,
&alpha12_temp,
gammaR,
sigmaR,
&alpha11_new );
// Apply the bulge-introducting Givens rotation (from the right)
// to the top-left 2x2 matrix.
MAC_Apply_G_2x2_opd( gammaR,
sigmaR,
alpha11,
alpha21,
alpha12,
alpha22 );
}
else
{
// Compute a new Givens rotation to push the bulge (from the
// right).
MAC_Givens2_opd( alpha01,
alpha02,
gammaR,
sigmaR,
&alpha01_new );
// Apply the Givens rotation (from the right) to the 1x2 vector
// from which it was computed, which annihilates alpha02.
*alpha01 = alpha01_new;
*alpha02 = 0.0;
// Apply the Givens rotation to the 2x2 matrix below the 1x2
// vector that was just modified.
MAC_Apply_G_2x2_opd( gammaR,
sigmaR,
alpha11,
alpha21,
alpha12,
alpha22 );
}
// Compute a new Givens rotation to push the bulge (from the left).
MAC_Givens2_opd( alpha11,
alpha21,
gammaL,
sigmaL,
&alpha11_new );
// Apply the Givens rotation (from the left) to the 2x1 vector
// from which it was computed, which annihilates alpha11.
*alpha11 = alpha11_new;
*alpha21 = 0.0;
if ( mn_ahead > 0 )
{
// Apply the Givens rotation (from the left) to the 2x2 submatrix
// at alpha12.
MAC_Apply_GT_2x2_opd( gammaL,
sigmaL,
alpha12,
alpha22,
alpha13,
alpha23 );
}
else
{
// Apply the Givens rotation (from the left) to the last 2x1 vector
// at alpha12.
MAC_Apply_GT_2x1_opd( gammaL,
sigmaL,
alpha12,
alpha22 );
}
/*------------------------------------------------------------*/
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Bsvd_francis_v_ops_var1 | ( | int | m_A, |
| float | shift, | ||
| scomplex * | buff_g, | ||
| int | inc_g, | ||
| scomplex * | buff_h, | ||
| int | inc_h, | ||
| float * | buff_d, | ||
| int | inc_d, | ||
| float * | buff_e, | ||
| int | inc_e | ||
| ) |
Referenced by FLA_Bsvd_francis_v_opt_var1().
{
FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );
return FLA_SUCCESS;
}
References FLA_Bsvd_francis_v_opd_var1(), FLA_Bsvd_francis_v_ops_var1(), FLA_Obj_datatype(), FLA_Obj_vector_dim(), and FLA_Obj_vector_inc().
{
FLA_Datatype datatype;
int m_A;
int inc_g;
int inc_h;
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_h = FLA_Obj_vector_inc( h );
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 );
scomplex* buff_h = FLA_COMPLEX_PTR( h );
float* buff_d = FLA_FLOAT_PTR( d );
float* buff_e = FLA_FLOAT_PTR( e );
FLA_Bsvd_francis_v_ops_var1( m_A,
*buff_shift,
buff_g, inc_g,
buff_h, inc_h,
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 );
dcomplex* buff_h = FLA_DOUBLE_COMPLEX_PTR( h );
double* buff_d = FLA_DOUBLE_PTR( d );
double* buff_e = FLA_DOUBLE_PTR( e );
FLA_Bsvd_francis_v_opd_var1( m_A,
*buff_shift,
buff_g, inc_g,
buff_h, inc_h,
buff_d, inc_d,
buff_e, inc_e );
break;
}
}
return FLA_SUCCESS;
}
1.7.6.1