|
libflame
12600
|
Go to the source code of this file.
Functions | |
| FLA_Error | FLA_Househ2_UT (FLA_Side side, FLA_Obj chi_1, FLA_Obj x2, FLA_Obj tau) |
| FLA_Error | FLA_Househ2_UT_l_ops (int m_x2, float *chi_1, float *x2, int inc_x2, float *tau) |
| FLA_Error | FLA_Househ2_UT_l_opd (int m_x2, double *chi_1, double *x2, int inc_x2, double *tau) |
| FLA_Error | FLA_Househ2_UT_l_opc (int m_x2, scomplex *chi_1, scomplex *x2, int inc_x2, scomplex *tau) |
| FLA_Error | FLA_Househ2_UT_l_opz (int m_x2, dcomplex *chi_1, dcomplex *x2, int inc_x2, dcomplex *tau) |
| FLA_Error | FLA_Househ2_UT_r_ops (int m_x2, float *chi_1, float *x2, int inc_x2, float *tau) |
| FLA_Error | FLA_Househ2_UT_r_opd (int m_x2, double *chi_1, double *x2, int inc_x2, double *tau) |
| FLA_Error | FLA_Househ2_UT_r_opc (int m_x2, scomplex *chi_1, scomplex *x2, int inc_x2, scomplex *tau) |
| FLA_Error | FLA_Househ2_UT_r_opz (int m_x2, dcomplex *chi_1, dcomplex *x2, int inc_x2, dcomplex *tau) |
| FLA_Error | FLA_Househ3UD_UT (FLA_Obj chi_1, FLA_Obj x2, FLA_Obj y2, FLA_Obj tau) |
| FLA_Error | FLA_Househ3UD_UT_ops (int m_x2, int m_y2, float *chi_1, float *x2, int inc_x2, float *y2, int inc_y2, float *tau) |
| FLA_Error | FLA_Househ3UD_UT_opd (int m_x2, int m_y2, double *chi_1, double *x2, int inc_x2, double *y2, int inc_y2, double *tau) |
| FLA_Error | FLA_Househ3UD_UT_opc (int m_x2, int m_y2, scomplex *chi_1, scomplex *x2, int inc_x2, scomplex *y2, int inc_y2, scomplex *tau) |
| FLA_Error | FLA_Househ3UD_UT_opz (int m_x2, int m_y2, dcomplex *chi_1, dcomplex *x2, int inc_x2, dcomplex *y2, int inc_y2, dcomplex *tau) |
| FLA_Error | FLA_Househ2s_UT (FLA_Side side, FLA_Obj chi_1, FLA_Obj x2, FLA_Obj alpha, FLA_Obj chi_1_minus_alpha, FLA_Obj tau) |
| FLA_Error | FLA_Househ2s_UT_l_ops (int m_x2, float *chi_1, float *x2, int inc_x2, float *alpha, float *chi_1_minus_alpha, float *tau) |
| FLA_Error | FLA_Househ2s_UT_l_opd (int m_x2, double *chi_1, double *x2, int inc_x2, double *alpha, double *chi_1_minus_alpha, double *tau) |
| FLA_Error | FLA_Househ2s_UT_l_opc (int m_x2, scomplex *chi_1, scomplex *x2, int inc_x2, scomplex *alpha, scomplex *chi_1_minus_alpha, scomplex *tau) |
| FLA_Error | FLA_Househ2s_UT_l_opz (int m_x2, dcomplex *chi_1, dcomplex *x2, int inc_x2, dcomplex *alpha, dcomplex *chi_1_minus_alpha, dcomplex *tau) |
| FLA_Error | FLA_Househ2s_UT_r_ops (int m_x2, float *chi_1, float *x2, int inc_x2, float *alpha, float *chi_1_minus_alpha, float *tau) |
| FLA_Error | FLA_Househ2s_UT_r_opd (int m_x2, double *chi_1, double *x2, int inc_x2, double *alpha, double *chi_1_minus_alpha, double *tau) |
| FLA_Error | FLA_Househ2s_UT_r_opc (int m_x2, scomplex *chi_1, scomplex *x2, int inc_x2, scomplex *alpha, scomplex *chi_1_minus_alpha, scomplex *tau) |
| FLA_Error | FLA_Househ2s_UT_r_opz (int m_x2, dcomplex *chi_1, dcomplex *x2, int inc_x2, dcomplex *alpha, dcomplex *chi_1_minus_alpha, dcomplex *tau) |
| FLA_Error | FLA_Hev_2x2 (FLA_Obj alpha11, FLA_Obj alpha21, FLA_Obj alpha22, FLA_Obj lambda1, FLA_Obj lambda2) |
| FLA_Error | FLA_Hev_2x2_ops (float *buff_alpha11, float *buff_alpha21, float *buff_alpha22, float *buff_lambda1, float *buff_lambda2) |
| FLA_Error | FLA_Hev_2x2_opd (double *buff_alpha11, double *buff_alpha21, double *buff_alpha22, double *buff_lambda1, double *buff_lambda2) |
| FLA_Error | FLA_Hevv_2x2 (FLA_Obj alpha11, FLA_Obj alpha21, FLA_Obj alpha22, FLA_Obj lambda1, FLA_Obj lambda2, FLA_Obj gamma1, FLA_Obj sigma1) |
| FLA_Error | FLA_Hevv_2x2_ops (float *alpha11, float *alpha21, float *alpha22, float *lambda1, float *lambda2, float *gamma1, float *sigma1) |
| FLA_Error | FLA_Hevv_2x2_opd (double *alpha11, double *alpha21, double *alpha22, double *lambda1, double *lambda2, double *gamma1, double *sigma1) |
| FLA_Error | FLA_Hevv_2x2_opc (scomplex *alpha11, scomplex *alpha21, scomplex *alpha22, float *lambda1, float *lambda2, float *gamma1, scomplex *sigma1) |
| FLA_Error | FLA_Hevv_2x2_opz (dcomplex *alpha11, dcomplex *alpha21, dcomplex *alpha22, double *lambda1, double *lambda2, double *gamma1, dcomplex *sigma1) |
| 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) |
| FLA_Error | FLA_Wilkshift_tridiag_opd (double delta1, double epsilon, double delta2, double *kappa) |
| FLA_Error | FLA_Pythag2 (FLA_Obj chi, FLA_Obj psi, FLA_Obj rho) |
| FLA_Error | FLA_Pythag2_ops (float *chi, float *psi, float *rho) |
| FLA_Error | FLA_Pythag2_opd (double *chi, double *psi, double *rho) |
| FLA_Error | FLA_Pythag3 (FLA_Obj chi, FLA_Obj psi, FLA_Obj zeta, FLA_Obj rho) |
| FLA_Error | FLA_Pythag3_ops (float *chi, float *psi, float *zeta, float *rho) |
| FLA_Error | FLA_Pythag3_opd (double *chi, double *psi, double *zeta, double *rho) |
| FLA_Error | FLA_Sort_evd (FLA_Direct direct, FLA_Obj l, FLA_Obj V) |
| FLA_Error | FLA_Sort_evd_f_ops (int m_A, float *l, int inc_l, float *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_evd_b_ops (int m_A, float *l, int inc_l, float *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_evd_f_opd (int m_A, double *l, int inc_l, double *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_evd_b_opd (int m_A, double *l, int inc_l, double *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_evd_f_opc (int m_A, float *l, int inc_l, scomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_evd_b_opc (int m_A, float *l, int inc_l, scomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_evd_f_opz (int m_A, double *l, int inc_l, dcomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_evd_b_opz (int m_A, double *l, int inc_l, dcomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_bsvd_ext (FLA_Direct direct, FLA_Obj s, FLA_Bool apply_U, FLA_Obj U, FLA_Bool apply_V, FLA_Obj V, FLA_Bool apply_C, FLA_Obj C) |
| FLA_Error | FLA_Sort_bsvd_ext_f_ops (int m_s, float *s, int inc_s, int m_U, float *U, int rs_U, int cs_U, int m_V, float *V, int rs_V, int cs_V, int n_C, float *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_bsvd_ext_b_ops (int m_s, float *s, int inc_s, int m_U, float *U, int rs_U, int cs_U, int m_V, float *V, int rs_V, int cs_V, int n_C, float *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_bsvd_ext_f_opd (int m_s, double *s, int inc_s, int m_U, double *U, int rs_U, int cs_U, int m_V, double *V, int rs_V, int cs_V, int n_C, double *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_bsvd_ext_b_opd (int m_s, double *s, int inc_s, int m_U, double *U, int rs_U, int cs_U, int m_V, double *V, int rs_V, int cs_V, int n_C, double *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_bsvd_ext_f_opc (int m_s, float *s, int inc_s, int m_U, scomplex *U, int rs_U, int cs_U, int m_V, scomplex *V, int rs_V, int cs_V, int n_C, scomplex *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_bsvd_ext_b_opc (int m_s, float *s, int inc_s, int m_U, scomplex *U, int rs_U, int cs_U, int m_V, scomplex *V, int rs_V, int cs_V, int n_C, scomplex *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_bsvd_ext_f_opz (int m_s, double *s, int inc_s, int m_U, dcomplex *U, int rs_U, int cs_U, int m_V, dcomplex *V, int rs_V, int cs_V, int n_C, dcomplex *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_bsvd_ext_b_opz (int m_s, double *s, int inc_s, int m_U, dcomplex *U, int rs_U, int cs_U, int m_V, dcomplex *V, int rs_V, int cs_V, int n_C, dcomplex *C, int rs_C, int cs_C) |
| FLA_Error | FLA_Sort_svd (FLA_Direct direct, FLA_Obj s, FLA_Obj U, FLA_Obj V) |
| FLA_Error | FLA_Sort_svd_f_ops (int m_U, int n_V, float *s, int inc_s, float *U, int rs_U, int cs_U, float *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_svd_b_ops (int m_U, int n_V, float *s, int inc_s, float *U, int rs_U, int cs_U, float *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_svd_f_opd (int m_U, int n_V, double *s, int inc_s, double *U, int rs_U, int cs_U, double *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_svd_b_opd (int m_U, int n_V, double *s, int inc_s, double *U, int rs_U, int cs_U, double *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_svd_f_opc (int m_U, int n_V, float *s, int inc_s, scomplex *U, int rs_U, int cs_U, scomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_svd_b_opc (int m_U, int n_V, float *s, int inc_s, scomplex *U, int rs_U, int cs_U, scomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_svd_f_opz (int m_U, int n_V, double *s, int inc_s, dcomplex *U, int rs_U, int cs_U, dcomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sort_svd_b_opz (int m_U, int n_V, double *s, int inc_s, dcomplex *U, int rs_U, int cs_U, dcomplex *V, int rs_V, int cs_V) |
| FLA_Error | FLA_Sv_2x2 (FLA_Obj alpha11, FLA_Obj alpha12, FLA_Obj alpha22, FLA_Obj sigma1, FLA_Obj sigma2) |
| FLA_Error | FLA_Sv_2x2_ops (float *alpha11, float *alpha12, float *alpha22, float *sigma1, float *sigma2) |
| FLA_Error | FLA_Sv_2x2_opd (double *alpha11, double *alpha12, double *alpha22, double *sigma1, double *sigma2) |
| FLA_Error | FLA_Svv_2x2 (FLA_Obj alpha11, FLA_Obj alpha12, FLA_Obj alpha22, FLA_Obj sigma1, FLA_Obj sigma2, FLA_Obj gammaL, FLA_Obj sigmaL, FLA_Obj gammaR, FLA_Obj sigmaR) |
| FLA_Error | FLA_Svv_2x2_ops (float *alpha11, float *alpha12, float *alpha22, float *sigma1, float *sigma2, float *gammaL, float *sigmaL, float *gammaR, float *sigmaR) |
| FLA_Error | FLA_Svv_2x2_opd (double *alpha11, double *alpha12, double *alpha22, double *sigma1, double *sigma2, double *gammaL, double *sigmaL, double *gammaR, double *sigmaR) |
| FLA_Error | FLA_Mach_params (FLA_Machval machval, FLA_Obj val) |
| float | FLA_Mach_params_ops (FLA_Machval machval) |
| double | FLA_Mach_params_opd (FLA_Machval machval) |
| FLA_Error | FLA_Apply_diag_matrix (FLA_Side side, FLA_Conj conj, FLA_Obj x, FLA_Obj A) |
| FLA_Error | FLA_Shift_pivots_to (FLA_Pivot_type ptype, FLA_Obj p) |
| FLA_Error | FLA_Form_perm_matrix (FLA_Obj p, FLA_Obj A) |
| FLA_Error | FLA_LU_find_zero_on_diagonal (FLA_Obj A) |
| doublereal | fla_dlamch (char *cmach, ftnlen cmach_len) |
| real | fla_slamch (char *cmach, ftnlen cmach_len) |
| logical | fla_lsame (char *ca, char *cb, ftnlen ca_len, ftnlen cb_len) |
| double | fla_pow_di (doublereal *a, integer *n) |
| real | fla_pow_ri (real *a, integer *n) |
| FLA_Error | FLA_Househ2_UT_check (FLA_Side side, FLA_Obj chi_1, FLA_Obj x2, FLA_Obj tau) |
| FLA_Error | FLA_Househ3UD_UT_check (FLA_Obj chi_1, FLA_Obj x2, FLA_Obj y2, FLA_Obj tau) |
| FLA_Error | FLA_Househ2s_UT_check (FLA_Side side, FLA_Obj chi_1, FLA_Obj x2, FLA_Obj alpha, FLA_Obj chi_1_minus_alpha, FLA_Obj tau) |
| FLA_Error | FLA_Givens2_check (FLA_Obj chi_1, FLA_Obj chi_2, FLA_Obj gamma, FLA_Obj sigma, FLA_Obj chi_1_new) |
| FLA_Error | FLA_Apply_GTG_check (FLA_Obj gamma, FLA_Obj sigma, FLA_Obj delta1, FLA_Obj epsilon1, FLA_Obj delta2) |
| FLA_Error | FLA_Apply_G_1x2_check (FLA_Obj gamma, FLA_Obj sigma, FLA_Obj beta, FLA_Obj epsilon) |
| FLA_Error | FLA_Apply_G_mx2_check (FLA_Obj gamma, FLA_Obj sigma, FLA_Obj a1, FLA_Obj a2) |
| FLA_Error | FLA_Apply_G_check (FLA_Side side, FLA_Direct direct, FLA_Obj G, FLA_Obj A) |
| FLA_Error | FLA_Wilkshift_tridiag_check (FLA_Obj delta1, FLA_Obj epsilon, FLA_Obj delta2, FLA_Obj kappa) |
| FLA_Error | FLA_Wilkshift_bidiag_check (FLA_Obj epsilon1, FLA_Obj delta1, FLA_Obj epsilon2, FLA_Obj delta2, FLA_Obj kappa) |
| FLA_Error | FLA_Introduce_bulge_check (FLA_Obj shift, FLA_Obj gamma, FLA_Obj sigma, FLA_Obj delta1, FLA_Obj epsilon1, FLA_Obj delta2, FLA_Obj beta, FLA_Obj epsilon2) |
| FLA_Error | FLA_Mach_params_check (FLA_Machval machval, FLA_Obj val) |
| FLA_Error | FLA_Sort_evd_check (FLA_Direct direct, FLA_Obj l, FLA_Obj V) |
| FLA_Error | FLA_Sort_svd_check (FLA_Direct direct, FLA_Obj s, FLA_Obj U, FLA_Obj V) |
| FLA_Error | FLA_Apply_diag_matrix_check (FLA_Side side, FLA_Conj conj, FLA_Obj x, FLA_Obj A) |
| FLA_Error | FLA_Shift_pivots_to_check (FLA_Pivot_type ptype, FLA_Obj p) |
| FLA_Error | FLA_Form_perm_matrix_check (FLA_Obj p, FLA_Obj A) |
| FLA_Error | FLA_LU_find_zero_on_diagonal_check (FLA_Obj A) |
References bl1_capdiagmv(), bl1_csapdiagmv(), bl1_dapdiagmv(), bl1_sapdiagmv(), bl1_zapdiagmv(), bl1_zdapdiagmv(), FLA_Apply_diag_matrix_check(), FLA_Check_error_level(), FLA_Obj_col_stride(), FLA_Obj_datatype(), FLA_Obj_length(), FLA_Obj_row_stride(), FLA_Obj_vector_inc(), FLA_Obj_width(), FLA_Param_map_flame_to_blis_conj(), and FLA_Param_map_flame_to_blis_side().
Referenced by FLA_Hevd_lv_unb_var1(), FLA_Hevd_lv_unb_var2(), FLA_Svd_ext_u_unb_var1(), FLA_Svd_uv_unb_var1(), and FLA_Svd_uv_unb_var2().
{
FLA_Datatype dt_x, dt_A;
int m_A, n_A;
int rs_A, cs_A;
int inc_x;
side1_t blis_side;
conj1_t blis_conj;
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Apply_diag_matrix_check( side, conj, x, A );
dt_x = FLA_Obj_datatype( x );
dt_A = FLA_Obj_datatype( A );
m_A = FLA_Obj_length( A );
n_A = FLA_Obj_width( A );
rs_A = FLA_Obj_row_stride( A );
cs_A = FLA_Obj_col_stride( A );
inc_x = FLA_Obj_vector_inc( x );
FLA_Param_map_flame_to_blis_side( side, &blis_side );
FLA_Param_map_flame_to_blis_conj( conj, &blis_conj );
switch ( dt_A )
{
case FLA_FLOAT:
{
float* buff_x = ( float* ) FLA_FLOAT_PTR( x );
float* buff_A = ( float* ) FLA_FLOAT_PTR( A );
bl1_sapdiagmv( blis_side,
blis_conj,
m_A,
n_A,
buff_x, inc_x,
buff_A, rs_A, cs_A );
break;
}
case FLA_DOUBLE:
{
double* buff_x = ( double* ) FLA_DOUBLE_PTR( x );
double* buff_A = ( double* ) FLA_DOUBLE_PTR( A );
bl1_dapdiagmv( blis_side,
blis_conj,
m_A,
n_A,
buff_x, inc_x,
buff_A, rs_A, cs_A );
break;
}
case FLA_COMPLEX:
{
if ( dt_x == FLA_FLOAT )
{
float* buff_x = ( float* ) FLA_FLOAT_PTR( x );
scomplex* buff_A = ( scomplex* ) FLA_COMPLEX_PTR( A );
bl1_csapdiagmv( blis_side,
blis_conj,
m_A,
n_A,
buff_x, inc_x,
buff_A, rs_A, cs_A );
}
else if ( dt_x == FLA_COMPLEX )
{
scomplex* buff_x = ( scomplex* ) FLA_COMPLEX_PTR( x );
scomplex* buff_A = ( scomplex* ) FLA_COMPLEX_PTR( A );
bl1_capdiagmv( blis_side,
blis_conj,
m_A,
n_A,
buff_x, inc_x,
buff_A, rs_A, cs_A );
}
break;
}
case FLA_DOUBLE_COMPLEX:
{
if ( dt_x == FLA_DOUBLE )
{
double* buff_x = ( double* ) FLA_DOUBLE_PTR( x );
dcomplex* buff_A = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( A );
bl1_zdapdiagmv( blis_side,
blis_conj,
m_A,
n_A,
buff_x, inc_x,
buff_A, rs_A, cs_A );
}
else if ( dt_x == FLA_DOUBLE_COMPLEX )
{
dcomplex* buff_x = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( x );
dcomplex* buff_A = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( A );
bl1_zapdiagmv( blis_side,
blis_conj,
m_A,
n_A,
buff_x, inc_x,
buff_A, rs_A, cs_A );
}
break;
}
}
return FLA_SUCCESS;
}
References FLA_Check_floating_object(), FLA_Check_identical_object_precision(), FLA_Check_nonconstant_object(), FLA_Check_object_length_equals(), FLA_Check_object_width_equals(), FLA_Check_valid_conj(), FLA_Check_valid_leftright_side(), and FLA_Obj_vector_dim().
Referenced by FLA_Apply_diag_matrix().
{
FLA_Error e_val;
e_val = FLA_Check_valid_leftright_side( side );
FLA_Check_error_code( e_val );
e_val = FLA_Check_valid_conj( conj );
FLA_Check_error_code( e_val );
e_val = FLA_Check_floating_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_precision( A, x );
FLA_Check_error_code( e_val );
if ( side == FLA_LEFT )
{
e_val = FLA_Check_object_length_equals( A, FLA_Obj_vector_dim( x ) );
FLA_Check_error_code( e_val );
}
else // if ( side == FLA_RIGHT )
{
e_val = FLA_Check_object_width_equals( A, FLA_Obj_vector_dim( x ) );
FLA_Check_error_code( e_val );
}
return FLA_SUCCESS;
}
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), FLA_Check_nonconstant_object(), and FLA_Check_real_object().
{
FLA_Error e_val;
e_val = FLA_Check_nonconstant_object( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, beta );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, epsilon );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( beta );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( epsilon );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Apply_G_check | ( | FLA_Side | side, |
| FLA_Direct | direct, | ||
| FLA_Obj | G, | ||
| FLA_Obj | A | ||
| ) |
References FLA_Check_complex_object(), FLA_Check_floating_object(), FLA_Check_identical_object_precision(), FLA_Check_nonconstant_object(), FLA_Check_object_length_equals(), FLA_Check_valid_direct(), FLA_Check_valid_leftright_side(), FLA_Obj_length(), and FLA_Obj_width().
Referenced by FLA_Apply_G().
{
FLA_Error e_val;
e_val = FLA_Check_valid_leftright_side( side );
FLA_Check_error_code( e_val );
e_val = FLA_Check_valid_direct( direct );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( G );
FLA_Check_error_code( e_val );
e_val = FLA_Check_complex_object( G );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_floating_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_precision( G, A );
FLA_Check_error_code( e_val );
if ( side == FLA_LEFT )
{
e_val = FLA_Check_object_length_equals( G, FLA_Obj_length( A ) - 1 );
FLA_Check_error_code( e_val );
}
else // if ( side == FLA_RIGHT )
{
e_val = FLA_Check_object_length_equals( G, FLA_Obj_width( A ) - 1 );
FLA_Check_error_code( e_val );
}
return FLA_SUCCESS;
}
References FLA_Check_equal_vector_dims(), FLA_Check_identical_object_datatype(), FLA_Check_identical_object_precision(), FLA_Check_if_scalar(), FLA_Check_if_vector(), FLA_Check_nonconstant_object(), and FLA_Check_real_object().
{
FLA_Error e_val;
e_val = FLA_Check_nonconstant_object( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( a1, a2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_precision( gamma, a1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_vector( a1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_vector( a2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_equal_vector_dims( a1, a2 );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Apply_GTG_check | ( | FLA_Obj | gamma, |
| FLA_Obj | sigma, | ||
| FLA_Obj | delta1, | ||
| FLA_Obj | epsilon1, | ||
| FLA_Obj | delta2 | ||
| ) |
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), FLA_Check_nonconstant_object(), and FLA_Check_real_object().
{
FLA_Error e_val;
e_val = FLA_Check_nonconstant_object( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, epsilon );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( gamma, delta2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( epsilon );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( delta2 );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| doublereal fla_dlamch | ( | char * | cmach, |
| ftnlen | cmach_len | ||
| ) |
References fla_dlamc2(), fla_lsame(), and fla_pow_di().
Referenced by FLA_Mach_params_opd().
{
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Builtin functions */
double fla_pow_di(doublereal *, integer *);
/* Local variables */
static doublereal base;
static integer beta;
static doublereal emin, prec, emax;
static integer imin, imax;
static logical lrnd;
static doublereal rmin, rmax, t, rmach;
extern logical fla_lsame(char *, char *, ftnlen, ftnlen);
static doublereal small, sfmin;
extern /* Subroutine */ int fla_dlamc2(integer *, integer *, logical *,
doublereal *, integer *, doublereal *, integer *, doublereal *);
static integer it;
static doublereal rnd, eps;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLAMCH determines double precision machine parameters. */
/* Arguments */
/* ========= */
/* CMACH (input) CHARACTER*1 */
/* Specifies the value to be returned by DLAMCH: */
/* = 'E' or 'e', DLAMCH := eps */
/* = 'S' or 's , DLAMCH := sfmin */
/* = 'B' or 'b', DLAMCH := base */
/* = 'P' or 'p', DLAMCH := eps*base */
/* = 'N' or 'n', DLAMCH := t */
/* = 'R' or 'r', DLAMCH := rnd */
/* = 'M' or 'm', DLAMCH := emin */
/* = 'U' or 'u', DLAMCH := rmin */
/* = 'L' or 'l', DLAMCH := emax */
/* = 'O' or 'o', DLAMCH := rmax */
/* where */
/* eps = relative machine precision */
/* sfmin = safe minimum, such that 1/sfmin does not overflow */
/* base = base of the machine */
/* prec = eps*base */
/* t = number of (base) digits in the mantissa */
/* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
/* emin = minimum exponent before (gradual) underflow */
/* rmin = underflow threshold - base**(emin-1) */
/* emax = largest exponent before overflow */
/* rmax = overflow threshold - (base**emax)*(1-eps) */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Save statement .. */
/* .. */
/* .. Data statements .. */
/* .. */
/* .. Executable Statements .. */
if (first) {
fla_dlamc2(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
base = (doublereal) beta;
t = (doublereal) it;
if (lrnd) {
rnd = 1.;
i__1 = 1 - it;
eps = fla_pow_di(&base, &i__1) / 2;
} else {
rnd = 0.;
i__1 = 1 - it;
eps = fla_pow_di(&base, &i__1);
}
prec = eps * base;
emin = (doublereal) imin;
emax = (doublereal) imax;
sfmin = rmin;
small = 1. / rmax;
if (small >= sfmin) {
/* Use SMALL plus a bit, to avoid the possibility of rounding */
/* causing overflow when computing 1/sfmin. */
sfmin = small * (eps + 1.);
}
}
if (fla_lsame(cmach, "E", (ftnlen)1, (ftnlen)1)) {
rmach = eps;
} else if (fla_lsame(cmach, "S", (ftnlen)1, (ftnlen)1)) {
rmach = sfmin;
} else if (fla_lsame(cmach, "B", (ftnlen)1, (ftnlen)1)) {
rmach = base;
} else if (fla_lsame(cmach, "P", (ftnlen)1, (ftnlen)1)) {
rmach = prec;
} else if (fla_lsame(cmach, "N", (ftnlen)1, (ftnlen)1)) {
rmach = t;
} else if (fla_lsame(cmach, "R", (ftnlen)1, (ftnlen)1)) {
rmach = rnd;
} else if (fla_lsame(cmach, "M", (ftnlen)1, (ftnlen)1)) {
rmach = emin;
} else if (fla_lsame(cmach, "U", (ftnlen)1, (ftnlen)1)) {
rmach = rmin;
} else if (fla_lsame(cmach, "L", (ftnlen)1, (ftnlen)1)) {
rmach = emax;
} else if (fla_lsame(cmach, "O", (ftnlen)1, (ftnlen)1)) {
rmach = rmax;
}
ret_val = rmach;
first = FALSE_;
return ret_val;
/* End of DLAMCH */
} /* fla_dlamch_ */
| FLA_Error FLA_Form_perm_matrix | ( | FLA_Obj | p, |
| FLA_Obj | A | ||
| ) |
References FLA_Apply_pivots(), FLA_Check_error_level(), FLA_Form_perm_matrix_check(), and FLA_Set_to_identity().
{
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Form_perm_matrix_check( p, A );
// We assume that A is correctly sized, m x m, where m is the row
// dimension of the matrix given to FLA_LU_piv() or similar function.
FLA_Set_to_identity( A );
// We assume that p contains pivots in native FLAME format. That is,
// we assume the pivot type is FLA_NATIVE_PIVOTS. This is not a huge
// assumption since the user has to go out of his way to shift the
// pivots into LAPACK-indexed pivots.
FLA_Apply_pivots( FLA_LEFT, FLA_NO_TRANSPOSE, p, A );
return FLA_SUCCESS;
}
References FLA_Check_floating_object(), FLA_Check_if_vector(), FLA_Check_int_object(), FLA_Check_matrix_vector_dims(), FLA_Check_nonconstant_object(), and FLA_Check_square().
Referenced by FLA_Form_perm_matrix().
{
FLA_Error e_val;
e_val = FLA_Check_int_object( p );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( p );
FLA_Check_error_code( e_val );
e_val = FLA_Check_floating_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_vector( p );
FLA_Check_error_code( e_val );
e_val = FLA_Check_square( A );
FLA_Check_error_code( e_val );
FLA_Check_matrix_vector_dims( FLA_NO_TRANSPOSE, A, p, p );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Givens2_check | ( | FLA_Obj | chi_1, |
| FLA_Obj | chi_2, | ||
| FLA_Obj | gamma, | ||
| FLA_Obj | sigma, | ||
| FLA_Obj | chi_1_new | ||
| ) |
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), FLA_Check_nonconstant_object(), and FLA_Check_real_object().
Referenced by FLA_Givens2().
{
FLA_Error e_val;
e_val = FLA_Check_nonconstant_object( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, chi_2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, chi_1_new );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( chi_2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( chi_1_new );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Hev_2x2 | ( | FLA_Obj | alpha11, |
| FLA_Obj | alpha21, | ||
| FLA_Obj | alpha22, | ||
| FLA_Obj | lambda1, | ||
| FLA_Obj | lambda2 | ||
| ) |
References FLA_Hev_2x2_opd(), FLA_Hev_2x2_ops(), and FLA_Obj_datatype().
{
FLA_Datatype datatype;
datatype = FLA_Obj_datatype( alpha11 );
switch ( datatype )
{
case FLA_FLOAT:
{
float* buff_alpha11 = FLA_FLOAT_PTR( alpha11 );
float* buff_alpha21 = FLA_FLOAT_PTR( alpha21 );
float* buff_alpha22 = FLA_FLOAT_PTR( alpha22 );
float* buff_lambda1 = FLA_FLOAT_PTR( lambda1 );
float* buff_lambda2 = FLA_FLOAT_PTR( lambda2 );
FLA_Hev_2x2_ops( buff_alpha11,
buff_alpha21,
buff_alpha22,
buff_lambda1,
buff_lambda2 );
break;
}
case FLA_DOUBLE:
{
double* buff_alpha11 = FLA_DOUBLE_PTR( alpha11 );
double* buff_alpha21 = FLA_DOUBLE_PTR( alpha21 );
double* buff_alpha22 = FLA_DOUBLE_PTR( alpha22 );
double* buff_lambda1 = FLA_DOUBLE_PTR( lambda1 );
double* buff_lambda2 = FLA_DOUBLE_PTR( lambda2 );
FLA_Hev_2x2_opd( buff_alpha11,
buff_alpha21,
buff_alpha22,
buff_lambda1,
buff_lambda2 );
break;
}
case FLA_COMPLEX:
{
FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );
break;
}
case FLA_DOUBLE_COMPLEX:
{
FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Hev_2x2_opd | ( | double * | buff_alpha11, |
| double * | buff_alpha21, | ||
| double * | buff_alpha22, | ||
| double * | buff_lambda1, | ||
| double * | buff_lambda2 | ||
| ) |
Referenced by FLA_Hev_2x2(), and FLA_Tevd_iteracc_n_opd_var1().
{
double a11, a21, a22;
double l1, l2;
double ab, acmn, acmx, adf, df, rt, sm, tb;
a11 = *alpha11;
a21 = *alpha21;
a22 = *alpha22;
sm = a11 + a22;
df = a11 - a22;
adf = fabs( df );
tb = a21 + a21;
ab = fabs( tb );
if ( fabs( a11 ) > fabs( a22 ) )
{
acmx = a11;
acmn = a22;
}
else
{
acmx = a22;
acmn = a11;
}
if ( adf > ab ) rt = adf * sqrt( 1.0 + ( ab / adf ) * ( ab / adf ) );
else if ( adf < ab ) rt = ab * sqrt( 1.0 + ( adf / ab ) * ( adf / ab ) );
else rt = ab * sqrt( 2.0 );
if ( sm < 0.0 )
{
l1 = 0.5 * ( sm - rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
}
else if ( sm > 0.0 )
{
l1 = 0.5 * ( sm + rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
}
else
{
l1 = 0.5 * rt;
l2 = -0.5 * rt;
}
*lambda1 = l1;
*lambda2 = l2;
return FLA_SUCCESS;
}
| FLA_Error FLA_Hev_2x2_ops | ( | float * | buff_alpha11, |
| float * | buff_alpha21, | ||
| float * | buff_alpha22, | ||
| float * | buff_lambda1, | ||
| float * | buff_lambda2 | ||
| ) |
Referenced by FLA_Hev_2x2().
{
float a11, a21, a22;
float l1, l2;
float ab, acmn, acmx, adf, df, rt, sm, tb;
a11 = *alpha11;
a21 = *alpha21;
a22 = *alpha22;
sm = a11 + a22;
df = a11 - a22;
adf = fabs( df );
tb = a21 + a21;
ab = fabs( tb );
if ( fabs( a11 ) > fabs( a22 ) )
{
acmx = a11;
acmn = a22;
}
else
{
acmx = a22;
acmn = a11;
}
if ( adf > ab ) rt = adf * sqrt( 1.0F + ( ab / adf ) * ( ab / adf ) );
else if ( adf < ab ) rt = ab * sqrt( 1.0F + ( adf / ab ) * ( adf / ab ) );
else rt = ab * sqrt( 2.0F );
if ( sm < 0.0F )
{
l1 = 0.5F * ( sm - rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
}
else if ( sm > 0.0F )
{
l1 = 0.5F * ( sm + rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
}
else
{
l1 = 0.5F * rt;
l2 = -0.5F * rt;
}
*lambda1 = l1;
*lambda2 = l2;
return FLA_SUCCESS;
}
| FLA_Error FLA_Hevv_2x2 | ( | FLA_Obj | alpha11, |
| FLA_Obj | alpha21, | ||
| FLA_Obj | alpha22, | ||
| FLA_Obj | lambda1, | ||
| FLA_Obj | lambda2, | ||
| FLA_Obj | gamma1, | ||
| FLA_Obj | sigma1 | ||
| ) |
References FLA_Hevv_2x2_opc(), FLA_Hevv_2x2_opd(), FLA_Hevv_2x2_ops(), FLA_Hevv_2x2_opz(), and FLA_Obj_datatype().
{
FLA_Datatype datatype;
datatype = FLA_Obj_datatype( alpha11 );
switch ( datatype )
{
case FLA_FLOAT:
{
float* buff_alpha11 = FLA_FLOAT_PTR( alpha11 );
float* buff_alpha21 = FLA_FLOAT_PTR( alpha21 );
float* buff_alpha22 = FLA_FLOAT_PTR( alpha22 );
float* buff_lambda1 = FLA_FLOAT_PTR( lambda1 );
float* buff_lambda2 = FLA_FLOAT_PTR( lambda2 );
float* buff_gamma1 = FLA_FLOAT_PTR( gamma1 );
float* buff_sigma1 = FLA_FLOAT_PTR( sigma1 );
FLA_Hevv_2x2_ops( buff_alpha11,
buff_alpha21,
buff_alpha22,
buff_lambda1,
buff_lambda2,
buff_gamma1,
buff_sigma1 );
break;
}
case FLA_DOUBLE:
{
double* buff_alpha11 = FLA_DOUBLE_PTR( alpha11 );
double* buff_alpha21 = FLA_DOUBLE_PTR( alpha21 );
double* buff_alpha22 = FLA_DOUBLE_PTR( alpha22 );
double* buff_lambda1 = FLA_DOUBLE_PTR( lambda1 );
double* buff_lambda2 = FLA_DOUBLE_PTR( lambda2 );
double* buff_gamma1 = FLA_DOUBLE_PTR( gamma1 );
double* buff_sigma1 = FLA_DOUBLE_PTR( sigma1 );
FLA_Hevv_2x2_opd( buff_alpha11,
buff_alpha21,
buff_alpha22,
buff_lambda1,
buff_lambda2,
buff_gamma1,
buff_sigma1 );
break;
}
case FLA_COMPLEX:
{
scomplex* buff_alpha11 = FLA_COMPLEX_PTR( alpha11 );
scomplex* buff_alpha21 = FLA_COMPLEX_PTR( alpha21 );
scomplex* buff_alpha22 = FLA_COMPLEX_PTR( alpha22 );
float* buff_lambda1 = FLA_FLOAT_PTR( lambda1 );
float* buff_lambda2 = FLA_FLOAT_PTR( lambda2 );
float* buff_gamma1 = FLA_FLOAT_PTR( gamma1 );
scomplex* buff_sigma1 = FLA_COMPLEX_PTR( sigma1 );
FLA_Hevv_2x2_opc( buff_alpha11,
buff_alpha21,
buff_alpha22,
buff_lambda1,
buff_lambda2,
buff_gamma1,
buff_sigma1 );
break;
}
case FLA_DOUBLE_COMPLEX:
{
dcomplex* buff_alpha11 = FLA_DOUBLE_COMPLEX_PTR( alpha11 );
dcomplex* buff_alpha21 = FLA_DOUBLE_COMPLEX_PTR( alpha21 );
dcomplex* buff_alpha22 = FLA_DOUBLE_COMPLEX_PTR( alpha22 );
double* buff_lambda1 = FLA_DOUBLE_PTR( lambda1 );
double* buff_lambda2 = FLA_DOUBLE_PTR( lambda2 );
double* buff_gamma1 = FLA_DOUBLE_PTR( gamma1 );
dcomplex* buff_sigma1 = FLA_DOUBLE_COMPLEX_PTR( sigma1 );
FLA_Hevv_2x2_opz( buff_alpha11,
buff_alpha21,
buff_alpha22,
buff_lambda1,
buff_lambda2,
buff_gamma1,
buff_sigma1 );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Hevv_2x2_opc | ( | scomplex * | alpha11, |
| scomplex * | alpha21, | ||
| scomplex * | alpha22, | ||
| float * | lambda1, | ||
| float * | lambda2, | ||
| float * | gamma1, | ||
| scomplex * | sigma1 | ||
| ) |
Referenced by FLA_Hevv_2x2().
{
FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );
return FLA_SUCCESS;
}
| FLA_Error FLA_Hevv_2x2_opd | ( | double * | alpha11, |
| double * | alpha21, | ||
| double * | alpha22, | ||
| double * | lambda1, | ||
| double * | lambda2, | ||
| double * | gamma1, | ||
| double * | sigma1 | ||
| ) |
Referenced by FLA_Hevv_2x2(), FLA_Tevd_iteracc_v_opd_var1(), and FLA_Tevd_iteracc_v_opd_var3().
{
double a11, a21, a22;
double l1, l2;
double g1, s1;
double ab, acmn, acmx, acs, adf, cs, ct, df, rt, sm, tb, tn;
int sgn1, sgn2;
a11 = *alpha11;
a21 = *alpha21;
a22 = *alpha22;
// Compute the eigenvalues.
sm = a11 + a22;
df = a11 - a22;
adf = fabs( df );
tb = a21 + a21;
ab = fabs( tb );
if ( fabs( a11 ) > fabs( a22 ) )
{
acmx = a11;
acmn = a22;
}
else
{
acmx = a22;
acmn = a11;
}
if ( adf > ab ) rt = adf * sqrt( 1.0 + pow( ( ab / adf ), 2.0 ) );
else if ( adf < ab ) rt = ab * sqrt( 1.0 + pow( ( adf / ab ), 2.0 ) );
else rt = ab * sqrt( 2.0 );
if ( sm < 0.0 )
{
l1 = 0.5 * ( sm - rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
sgn1 = -1;
}
else if ( sm > 0.0 )
{
l1 = 0.5 * ( sm + rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
sgn1 = 1;
}
else
{
l1 = 0.5 * rt;
l2 = -0.5 * rt;
sgn1 = 1;
}
*lambda1 = l1;
*lambda2 = l2;
// Compute the eigenvector.
if ( df >= 0.0 )
{
cs = df + rt;
sgn2 = 1;
}
else
{
cs = df - rt;
sgn2 = -1;
}
acs = fabs( cs );
if ( acs > ab )
{
ct = -tb / cs;
s1 = 1.0 / sqrt( 1.0 + ct*ct );
g1 = ct * s1;
}
else
{
if ( ab == 0.0 )
{
g1 = 1.0;
s1 = 0.0;
}
else
{
tn = -cs / tb;
g1 = 1.0 / sqrt( 1.0 + tn*tn );
s1 = tn * g1;
}
}
if ( sgn1 == sgn2 )
{
tn = g1;
g1 = -s1;
s1 = tn;
}
*gamma1 = g1;
*sigma1 = s1;
return FLA_SUCCESS;
}
| FLA_Error FLA_Hevv_2x2_ops | ( | float * | alpha11, |
| float * | alpha21, | ||
| float * | alpha22, | ||
| float * | lambda1, | ||
| float * | lambda2, | ||
| float * | gamma1, | ||
| float * | sigma1 | ||
| ) |
Referenced by FLA_Hevv_2x2().
{
float a11, a21, a22;
float l1, l2;
float g1, s1;
float ab, acmn, acmx, acs, adf, cs, ct, df, rt, sm, tb, tn;
int sgn1, sgn2;
a11 = *alpha11;
a21 = *alpha21;
a22 = *alpha22;
// Compute the eigenvalues.
sm = a11 + a22;
df = a11 - a22;
adf = fabs( df );
tb = a21 + a21;
ab = fabs( tb );
if ( fabs( a11 ) > fabs( a22 ) )
{
acmx = a11;
acmn = a22;
}
else
{
acmx = a22;
acmn = a11;
}
if ( adf > ab ) rt = adf * sqrt( 1.0F + ( ab / adf ) * ( ab / adf ) );
else if ( adf < ab ) rt = ab * sqrt( 1.0F + ( adf / ab ) * ( adf / ab ) );
else rt = ab * sqrt( 2.0F );
if ( sm < 0.0F )
{
l1 = 0.5F * ( sm - rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
sgn1 = -1;
}
else if ( sm > 0.0F )
{
l1 = 0.5F * ( sm + rt );
l2 = ( acmx / l1 ) * acmn - ( a21 / l1 ) * a21;
sgn1 = 1;
}
else
{
l1 = 0.5F * rt;
l2 = -0.5F * rt;
sgn1 = 1;
}
*lambda1 = l1;
*lambda2 = l2;
// Compute the eigenvector.
if ( df >= 0.0F )
{
cs = df + rt;
sgn2 = 1;
}
else
{
cs = df - rt;
sgn2 = -1;
}
acs = fabs( cs );
if ( acs > ab )
{
ct = -tb / cs;
s1 = 1.0F / sqrt( 1.0F + ct*ct );
g1 = ct * s1;
}
else
{
if ( ab == 0.0F )
{
g1 = 1.0F;
s1 = 0.0F;
}
else
{
tn = -cs / tb;
g1 = 1.0F / sqrt( 1.0F + tn*tn );
s1 = tn * g1;
}
}
if ( sgn1 == sgn2 )
{
tn = g1;
g1 = -s1;
s1 = tn;
}
*gamma1 = g1;
*sigma1 = s1;
return FLA_SUCCESS;
}
| FLA_Error FLA_Hevv_2x2_opz | ( | dcomplex * | alpha11, |
| dcomplex * | alpha21, | ||
| dcomplex * | alpha22, | ||
| double * | lambda1, | ||
| double * | lambda2, | ||
| double * | gamma1, | ||
| dcomplex * | sigma1 | ||
| ) |
Referenced by FLA_Hevv_2x2().
{
FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );
return FLA_SUCCESS;
}
References FLA_Check_error_level(), FLA_Househ2_UT_check(), FLA_Househ2_UT_l_opc(), FLA_Househ2_UT_l_opd(), FLA_Househ2_UT_l_ops(), FLA_Househ2_UT_l_opz(), FLA_Househ2_UT_r_opc(), FLA_Househ2_UT_r_opd(), FLA_Househ2_UT_r_ops(), FLA_Househ2_UT_r_opz(), FLA_Obj_datatype(), FLA_Obj_vector_dim(), and FLA_Obj_vector_inc().
Referenced by FLA_Bidiag_UT_u_step_unb_var1(), FLA_Bidiag_UT_u_step_unb_var2(), FLA_Bidiag_UT_u_step_unb_var3(), FLA_Bidiag_UT_u_step_unb_var4(), FLA_Bidiag_UT_u_step_unb_var5(), FLA_CAQR2_UT_unb_var1(), FLA_Hess_UT_step_unb_var1(), FLA_Hess_UT_step_unb_var2(), FLA_Hess_UT_step_unb_var3(), FLA_Hess_UT_step_unb_var4(), FLA_Hess_UT_step_unb_var5(), FLA_LQ_UT_unb_var1(), FLA_LQ_UT_unb_var2(), FLA_QR2_UT_unb_var1(), FLA_QR_UT_piv_unb_var1(), FLA_QR_UT_piv_unb_var2(), FLA_QR_UT_unb_var1(), FLA_QR_UT_unb_var2(), FLA_Tridiag_UT_l_step_unb_var1(), FLA_Tridiag_UT_l_step_unb_var2(), and FLA_Tridiag_UT_l_step_unb_var3().
{
FLA_Datatype datatype;
int m_x2;
int inc_x2;
datatype = FLA_Obj_datatype( x2 );
m_x2 = FLA_Obj_vector_dim( x2 );
inc_x2 = FLA_Obj_vector_inc( x2 );
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Househ2_UT_check( side, chi_1, x2, tau );
switch ( datatype )
{
case FLA_FLOAT:
{
float* chi_1_p = ( float* ) FLA_FLOAT_PTR( chi_1 );
float* x2_p = ( float* ) FLA_FLOAT_PTR( x2 );
float* tau_p = ( float* ) FLA_FLOAT_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2_UT_l_ops( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2_UT_r_ops( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
break;
}
case FLA_DOUBLE:
{
double* chi_1_p = ( double* ) FLA_DOUBLE_PTR( chi_1 );
double* x2_p = ( double* ) FLA_DOUBLE_PTR( x2 );
double* tau_p = ( double* ) FLA_DOUBLE_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2_UT_l_opd( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2_UT_r_opd( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
break;
}
case FLA_COMPLEX:
{
scomplex* chi_1_p = ( scomplex* ) FLA_COMPLEX_PTR( chi_1 );
scomplex* x2_p = ( scomplex* ) FLA_COMPLEX_PTR( x2 );
scomplex* tau_p = ( scomplex* ) FLA_COMPLEX_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2_UT_l_opc( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2_UT_r_opc( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
break;
}
case FLA_DOUBLE_COMPLEX:
{
dcomplex* chi_1_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( chi_1 );
dcomplex* x2_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( x2 );
dcomplex* tau_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2_UT_l_opz( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2_UT_r_opz( m_x2,
chi_1_p,
x2_p, inc_x2,
tau_p );
break;
}
}
return FLA_SUCCESS;
}
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), FLA_Check_nonconstant_object(), and FLA_Check_valid_leftright_side().
Referenced by FLA_Househ2_UT().
{
FLA_Error e_val;
e_val = FLA_Check_valid_leftright_side( side );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, x2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, tau );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( tau );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_l_opc | ( | int | m_x2, |
| scomplex * | chi_1, | ||
| scomplex * | x2, | ||
| int | inc_x2, | ||
| scomplex * | tau | ||
| ) |
References bl1_cinvscalv(), bl1_cnrm2(), BLIS1_NO_CONJUGATE, FLA_ONE_HALF, scomplex::imag, and scomplex::real.
Referenced by FLA_Bidiag_UT_u_step_ofc_var2(), FLA_Bidiag_UT_u_step_ofc_var3(), FLA_Bidiag_UT_u_step_ofc_var4(), FLA_Bidiag_UT_u_step_opc_var1(), FLA_Bidiag_UT_u_step_opc_var2(), FLA_Bidiag_UT_u_step_opc_var3(), FLA_Bidiag_UT_u_step_opc_var4(), FLA_Bidiag_UT_u_step_opc_var5(), FLA_CAQR2_UT_opc_var1(), FLA_Hess_UT_step_ofc_var2(), FLA_Hess_UT_step_ofc_var3(), FLA_Hess_UT_step_ofc_var4(), FLA_Hess_UT_step_opc_var1(), FLA_Hess_UT_step_opc_var2(), FLA_Hess_UT_step_opc_var3(), FLA_Hess_UT_step_opc_var4(), FLA_Hess_UT_step_opc_var5(), FLA_Househ2_UT(), FLA_Househ2_UT_r_opc(), FLA_QR2_UT_opc_var1(), FLA_QR_UT_opc_var1(), FLA_QR_UT_opc_var2(), FLA_Tridiag_UT_l_step_ofc_var2(), FLA_Tridiag_UT_l_step_ofc_var3(), FLA_Tridiag_UT_l_step_opc_var1(), FLA_Tridiag_UT_l_step_opc_var2(), and FLA_Tridiag_UT_l_step_opc_var3().
{
scomplex one_half = *FLA_COMPLEX_PTR( FLA_ONE_HALF );
scomplex y[2];
scomplex alpha;
scomplex chi_1_minus_alpha;
float abs_chi_1;
float norm_x_2;
float norm_x;
float abs_chi_1_minus_alpha;
float norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_cnrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0F )
{
chi_1->real = -(chi_1->real);
chi_1->imag = -(chi_1->imag);
tau->real = one_half.real;
tau->imag = one_half.imag;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_cnrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0].real = abs_chi_1;
y[0].imag = 0.0F;
y[1].real = norm_x_2;
y[1].imag = 0.0F;
bl1_cnrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
//
if ( abs_chi_1 == 0.0F )
{
alpha.real = norm_x * ( -1.0F );
alpha.imag = norm_x * ( -1.0F );
}
else
{
alpha.real = norm_x * ( -chi_1->real / abs_chi_1 );
alpha.imag = norm_x * ( -chi_1->imag / abs_chi_1 );
}
//
// Overwrite x_2 with u_2:
//
// x_2 := x_2 / ( chi_1 - alpha )
//
chi_1_minus_alpha.real = chi_1->real - alpha.real;
chi_1_minus_alpha.imag = chi_1->imag - alpha.imag;
bl1_cinvscalv( BLIS1_NO_CONJUGATE,
m_x2,
&chi_1_minus_alpha,
x2, inc_x2 );
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_csabsval2( &chi_1_minus_alpha, &abs_chi_1_minus_alpha );
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
tau->real = one_half.real + one_half.real*(norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha);
tau->imag = 0.0F;
//
// Overwrite chi_1 with alpha:
//
// chi_1 := alpha
//
chi_1->real = alpha.real;
chi_1->imag = alpha.imag;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_l_opd | ( | int | m_x2, |
| double * | chi_1, | ||
| double * | x2, | ||
| int | inc_x2, | ||
| double * | tau | ||
| ) |
References bl1_dinvscalv(), bl1_dnrm2(), BLIS1_NO_CONJUGATE, and FLA_ONE_HALF.
Referenced by FLA_Bidiag_UT_u_step_ofd_var2(), FLA_Bidiag_UT_u_step_ofd_var3(), FLA_Bidiag_UT_u_step_ofd_var4(), FLA_Bidiag_UT_u_step_opd_var1(), FLA_Bidiag_UT_u_step_opd_var2(), FLA_Bidiag_UT_u_step_opd_var3(), FLA_Bidiag_UT_u_step_opd_var4(), FLA_Bidiag_UT_u_step_opd_var5(), FLA_CAQR2_UT_opd_var1(), FLA_Hess_UT_step_ofd_var2(), FLA_Hess_UT_step_ofd_var3(), FLA_Hess_UT_step_ofd_var4(), FLA_Hess_UT_step_opd_var1(), FLA_Hess_UT_step_opd_var2(), FLA_Hess_UT_step_opd_var3(), FLA_Hess_UT_step_opd_var4(), FLA_Hess_UT_step_opd_var5(), FLA_Househ2_UT(), FLA_Househ2_UT_r_opd(), FLA_QR2_UT_opd_var1(), FLA_QR_UT_opd_var1(), FLA_QR_UT_opd_var2(), FLA_Tridiag_UT_l_step_ofd_var2(), FLA_Tridiag_UT_l_step_ofd_var3(), FLA_Tridiag_UT_l_step_opd_var1(), FLA_Tridiag_UT_l_step_opd_var2(), and FLA_Tridiag_UT_l_step_opd_var3().
{
double one_half = *FLA_DOUBLE_PTR( FLA_ONE_HALF );
double y[2];
double alpha;
double chi_1_minus_alpha;
double abs_chi_1;
double norm_x_2;
double norm_x;
double abs_chi_1_minus_alpha;
double norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_dnrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0 )
{
*chi_1 = -(*chi_1);
*tau = one_half;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_dnrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0] = abs_chi_1;
y[1] = norm_x_2;
bl1_dnrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
// = -sign( chi_1 ) * || x ||_2
//
alpha = -dsign( *chi_1 ) * norm_x;
//
// Overwrite x_2 with u_2:
//
// x_2 := x_2 / ( chi_1 - alpha )
//
chi_1_minus_alpha = *chi_1 - alpha;
bl1_dinvscalv( BLIS1_NO_CONJUGATE,
m_x2,
&chi_1_minus_alpha,
x2, inc_x2 );
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_dabsval2( &chi_1_minus_alpha, &abs_chi_1_minus_alpha );
abs_chi_1_minus_alpha = (double)fabs(chi_1_minus_alpha);
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
*tau = one_half + one_half*(norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha);
//
// Overwrite chi_1 with alpha:
//
// chi_1 := alpha
//
*chi_1 = alpha;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_l_ops | ( | int | m_x2, |
| float * | chi_1, | ||
| float * | x2, | ||
| int | inc_x2, | ||
| float * | tau | ||
| ) |
References bl1_sinvscalv(), bl1_snrm2(), BLIS1_NO_CONJUGATE, and FLA_ONE_HALF.
Referenced by FLA_Bidiag_UT_u_step_ofs_var2(), FLA_Bidiag_UT_u_step_ofs_var3(), FLA_Bidiag_UT_u_step_ofs_var4(), FLA_Bidiag_UT_u_step_ops_var1(), FLA_Bidiag_UT_u_step_ops_var2(), FLA_Bidiag_UT_u_step_ops_var3(), FLA_Bidiag_UT_u_step_ops_var4(), FLA_Bidiag_UT_u_step_ops_var5(), FLA_CAQR2_UT_ops_var1(), FLA_Hess_UT_step_ofs_var2(), FLA_Hess_UT_step_ofs_var3(), FLA_Hess_UT_step_ofs_var4(), FLA_Hess_UT_step_ops_var1(), FLA_Hess_UT_step_ops_var2(), FLA_Hess_UT_step_ops_var3(), FLA_Hess_UT_step_ops_var4(), FLA_Hess_UT_step_ops_var5(), FLA_Househ2_UT(), FLA_Househ2_UT_r_ops(), FLA_QR2_UT_ops_var1(), FLA_QR_UT_ops_var1(), FLA_QR_UT_ops_var2(), FLA_Tridiag_UT_l_step_ofs_var2(), FLA_Tridiag_UT_l_step_ofs_var3(), FLA_Tridiag_UT_l_step_ops_var1(), FLA_Tridiag_UT_l_step_ops_var2(), and FLA_Tridiag_UT_l_step_ops_var3().
{
float one_half = *FLA_FLOAT_PTR( FLA_ONE_HALF );
float y[2];
float alpha;
float chi_1_minus_alpha;
float abs_chi_1;
float norm_x_2;
float norm_x;
float abs_chi_1_minus_alpha;
float norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_snrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0F )
{
*chi_1 = -(*chi_1);
*tau = one_half;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_snrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0] = abs_chi_1;
y[1] = norm_x_2;
bl1_snrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
// = -sign( chi_1 ) * || x ||_2
//
alpha = -ssign( *chi_1 ) * norm_x;
//
// Overwrite x_2 with u_2:
//
// x_2 := x_2 / ( chi_1 - alpha )
//
chi_1_minus_alpha = *chi_1 - alpha;
bl1_sinvscalv( BLIS1_NO_CONJUGATE,
m_x2,
&chi_1_minus_alpha,
x2, inc_x2 );
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_sabsval2( &chi_1_minus_alpha, &abs_chi_1_minus_alpha );
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
*tau = one_half + one_half*(norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha);
//
// Overwrite chi_1 with alpha:
//
// chi_1 := alpha
//
*chi_1 = alpha;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_l_opz | ( | int | m_x2, |
| dcomplex * | chi_1, | ||
| dcomplex * | x2, | ||
| int | inc_x2, | ||
| dcomplex * | tau | ||
| ) |
References bl1_zinvscalv(), bl1_znrm2(), BLIS1_NO_CONJUGATE, FLA_ONE_HALF, dcomplex::imag, and dcomplex::real.
Referenced by FLA_Bidiag_UT_u_step_ofz_var2(), FLA_Bidiag_UT_u_step_ofz_var3(), FLA_Bidiag_UT_u_step_ofz_var4(), FLA_Bidiag_UT_u_step_opz_var1(), FLA_Bidiag_UT_u_step_opz_var2(), FLA_Bidiag_UT_u_step_opz_var3(), FLA_Bidiag_UT_u_step_opz_var4(), FLA_Bidiag_UT_u_step_opz_var5(), FLA_CAQR2_UT_opz_var1(), FLA_Hess_UT_step_ofz_var2(), FLA_Hess_UT_step_ofz_var3(), FLA_Hess_UT_step_ofz_var4(), FLA_Hess_UT_step_opz_var1(), FLA_Hess_UT_step_opz_var2(), FLA_Hess_UT_step_opz_var3(), FLA_Hess_UT_step_opz_var4(), FLA_Hess_UT_step_opz_var5(), FLA_Househ2_UT(), FLA_Househ2_UT_r_opz(), FLA_QR2_UT_opz_var1(), FLA_QR_UT_opz_var1(), FLA_QR_UT_opz_var2(), FLA_Tridiag_UT_l_step_ofz_var2(), FLA_Tridiag_UT_l_step_ofz_var3(), FLA_Tridiag_UT_l_step_opz_var1(), FLA_Tridiag_UT_l_step_opz_var2(), and FLA_Tridiag_UT_l_step_opz_var3().
{
dcomplex one_half = *FLA_DOUBLE_COMPLEX_PTR( FLA_ONE_HALF );
dcomplex y[2];
dcomplex alpha;
dcomplex chi_1_minus_alpha;
double abs_chi_1;
double norm_x_2;
double norm_x;
double abs_chi_1_minus_alpha;
double norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_znrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0 )
{
chi_1->real = -(chi_1->real);
chi_1->imag = -(chi_1->imag);
tau->real = one_half.real;
tau->imag = one_half.imag;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_znrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0].real = abs_chi_1;
y[0].imag = 0.0;
y[1].real = norm_x_2;
y[1].imag = 0.0;
bl1_znrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
//
if ( abs_chi_1 == 0.0 )
{
alpha.real = norm_x * ( -1.0 );
alpha.imag = norm_x * ( -1.0 );
}
else
{
alpha.real = norm_x * ( -chi_1->real / abs_chi_1 );
alpha.imag = norm_x * ( -chi_1->imag / abs_chi_1 );
}
//
// Overwrite x_2 with u_2:
//
// x_2 := x_2 / ( chi_1 - alpha )
//
chi_1_minus_alpha.real = chi_1->real - alpha.real;
chi_1_minus_alpha.imag = chi_1->imag - alpha.imag;
bl1_zinvscalv( BLIS1_NO_CONJUGATE,
m_x2,
&chi_1_minus_alpha,
x2, inc_x2 );
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_zdabsval2( &chi_1_minus_alpha, &abs_chi_1_minus_alpha );
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
tau->real = one_half.real + one_half.real*(norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha);
tau->imag = 0.0;
//
// Overwrite chi_1 with alpha:
//
// chi_1 := alpha
//
chi_1->real = alpha.real;
chi_1->imag = alpha.imag;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_r_opc | ( | int | m_x2, |
| scomplex * | chi_1, | ||
| scomplex * | x2, | ||
| int | inc_x2, | ||
| scomplex * | tau | ||
| ) |
References bl1_cconjv(), and FLA_Househ2_UT_l_opc().
Referenced by FLA_Bidiag_UT_u_step_ofc_var2(), FLA_Bidiag_UT_u_step_opc_var1(), FLA_Bidiag_UT_u_step_opc_var2(), FLA_Bidiag_UT_u_step_opc_var5(), FLA_Househ2_UT(), FLA_LQ_UT_opc_var1(), and FLA_LQ_UT_opc_var2().
{
FLA_Househ2_UT_l_opc( m_x2,
chi_1,
x2, inc_x2,
tau );
bl1_cconjv( m_x2,
x2, inc_x2 );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_r_opd | ( | int | m_x2, |
| double * | chi_1, | ||
| double * | x2, | ||
| int | inc_x2, | ||
| double * | tau | ||
| ) |
References FLA_Househ2_UT_l_opd().
Referenced by FLA_Bidiag_UT_u_step_ofd_var2(), FLA_Bidiag_UT_u_step_opd_var1(), FLA_Bidiag_UT_u_step_opd_var2(), FLA_Bidiag_UT_u_step_opd_var5(), FLA_Househ2_UT(), FLA_LQ_UT_opd_var1(), and FLA_LQ_UT_opd_var2().
{
FLA_Househ2_UT_l_opd( m_x2,
chi_1,
x2, inc_x2,
tau );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_r_ops | ( | int | m_x2, |
| float * | chi_1, | ||
| float * | x2, | ||
| int | inc_x2, | ||
| float * | tau | ||
| ) |
References FLA_Househ2_UT_l_ops().
Referenced by FLA_Bidiag_UT_u_step_ofs_var2(), FLA_Bidiag_UT_u_step_ops_var1(), FLA_Bidiag_UT_u_step_ops_var2(), FLA_Bidiag_UT_u_step_ops_var5(), FLA_Househ2_UT(), FLA_LQ_UT_ops_var1(), and FLA_LQ_UT_ops_var2().
{
FLA_Househ2_UT_l_ops( m_x2,
chi_1,
x2, inc_x2,
tau );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2_UT_r_opz | ( | int | m_x2, |
| dcomplex * | chi_1, | ||
| dcomplex * | x2, | ||
| int | inc_x2, | ||
| dcomplex * | tau | ||
| ) |
References bl1_zconjv(), and FLA_Househ2_UT_l_opz().
Referenced by FLA_Bidiag_UT_u_step_ofz_var2(), FLA_Bidiag_UT_u_step_opz_var1(), FLA_Bidiag_UT_u_step_opz_var2(), FLA_Bidiag_UT_u_step_opz_var5(), FLA_Househ2_UT(), FLA_LQ_UT_opz_var1(), and FLA_LQ_UT_opz_var2().
{
FLA_Househ2_UT_l_opz( m_x2,
chi_1,
x2, inc_x2,
tau );
bl1_zconjv( m_x2,
x2, inc_x2 );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT | ( | FLA_Side | side, |
| FLA_Obj | chi_1, | ||
| FLA_Obj | x2, | ||
| FLA_Obj | alpha, | ||
| FLA_Obj | chi_1_minus_alpha, | ||
| FLA_Obj | tau | ||
| ) |
References FLA_Check_error_level(), FLA_Househ2s_UT_check(), FLA_Househ2s_UT_l_opc(), FLA_Househ2s_UT_l_opd(), FLA_Househ2s_UT_l_ops(), FLA_Househ2s_UT_l_opz(), FLA_Househ2s_UT_r_opc(), FLA_Househ2s_UT_r_opd(), FLA_Househ2s_UT_r_ops(), FLA_Househ2s_UT_r_opz(), FLA_Obj_datatype(), FLA_Obj_vector_dim(), and FLA_Obj_vector_inc().
Referenced by FLA_Bidiag_UT_u_step_unb_var3(), and FLA_Bidiag_UT_u_step_unb_var4().
{
FLA_Datatype datatype;
int m_x2;
int inc_x2;
datatype = FLA_Obj_datatype( x2 );
m_x2 = FLA_Obj_vector_dim( x2 );
inc_x2 = FLA_Obj_vector_inc( x2 );
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Househ2s_UT_check( side, chi_1, x2, alpha, chi_1_minus_alpha, tau );
switch ( datatype )
{
case FLA_FLOAT:
{
float* chi_1_p = ( float* ) FLA_FLOAT_PTR( chi_1 );
float* x2_p = ( float* ) FLA_FLOAT_PTR( x2 );
float* alpha_p = ( float* ) FLA_FLOAT_PTR( alpha );
float* chi_1_minus_alpha_p = ( float* ) FLA_FLOAT_PTR( chi_1_minus_alpha );
float* tau_p = ( float* ) FLA_FLOAT_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2s_UT_l_ops( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2s_UT_r_ops( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
break;
}
case FLA_DOUBLE:
{
double* chi_1_p = ( double* ) FLA_DOUBLE_PTR( chi_1 );
double* x2_p = ( double* ) FLA_DOUBLE_PTR( x2 );
double* alpha_p = ( double* ) FLA_DOUBLE_PTR( alpha );
double* chi_1_minus_alpha_p = ( double* ) FLA_DOUBLE_PTR( chi_1_minus_alpha );
double* tau_p = ( double* ) FLA_DOUBLE_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2s_UT_l_opd( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2s_UT_r_opd( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
break;
}
case FLA_COMPLEX:
{
scomplex* chi_1_p = ( scomplex* ) FLA_COMPLEX_PTR( chi_1 );
scomplex* x2_p = ( scomplex* ) FLA_COMPLEX_PTR( x2 );
scomplex* alpha_p = ( scomplex* ) FLA_COMPLEX_PTR( alpha );
scomplex* chi_1_minus_alpha_p = ( scomplex* ) FLA_COMPLEX_PTR( chi_1_minus_alpha );
scomplex* tau_p = ( scomplex* ) FLA_COMPLEX_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2s_UT_l_opc( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2s_UT_r_opc( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
break;
}
case FLA_DOUBLE_COMPLEX:
{
dcomplex* chi_1_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( chi_1 );
dcomplex* x2_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( x2 );
dcomplex* alpha_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( alpha );
dcomplex* chi_1_minus_alpha_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( chi_1_minus_alpha );
dcomplex* tau_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( tau );
if ( side == FLA_LEFT )
FLA_Househ2s_UT_l_opz( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
else // if ( side == FLA_RIGHT )
FLA_Househ2s_UT_r_opz( m_x2,
chi_1_p,
x2_p, inc_x2,
alpha_p,
chi_1_minus_alpha_p,
tau_p );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_check | ( | FLA_Side | side, |
| FLA_Obj | chi_1, | ||
| FLA_Obj | x2, | ||
| FLA_Obj | alpha, | ||
| FLA_Obj | chi_1_minus_alpha, | ||
| FLA_Obj | tau | ||
| ) |
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), FLA_Check_if_vector(), FLA_Check_nonconstant_object(), and FLA_Check_valid_leftright_side().
Referenced by FLA_Househ2s_UT().
{
FLA_Error e_val;
e_val = FLA_Check_valid_leftright_side( side );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, x2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, alpha );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, chi_1_minus_alpha );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, tau );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_vector( x2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( alpha );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( chi_1_minus_alpha );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( tau );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_l_opc | ( | int | m_x2, |
| scomplex * | chi_1, | ||
| scomplex * | x2, | ||
| int | inc_x2, | ||
| scomplex * | alpha, | ||
| scomplex * | chi_1_minus_alpha, | ||
| scomplex * | tau | ||
| ) |
References bl1_cnrm2(), FLA_ONE_HALF, scomplex::imag, and scomplex::real.
Referenced by FLA_Househ2s_UT(), and FLA_Househ2s_UT_r_opc().
{
scomplex one_half = *FLA_COMPLEX_PTR( FLA_ONE_HALF );
scomplex y[2];
float abs_chi_1;
float norm_x_2;
float norm_x;
float abs_chi_1_minus_alpha;
float norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_cnrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0F )
{
alpha->real = -(chi_1->real);
alpha->imag = -(chi_1->imag);
chi_1_minus_alpha->real = 2.0F * chi_1->real;
chi_1_minus_alpha->imag = 2.0F * chi_1->imag;
tau->real = one_half.real;
tau->imag = one_half.imag;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_cnrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0].real = abs_chi_1;
y[0].imag = 0.0;
y[1].real = norm_x_2;
y[1].imag = 0.0F;
bl1_cnrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
//
if ( abs_chi_1 == 0.0F )
{
alpha->real = norm_x * ( -1.0F );
alpha->imag = norm_x * ( -1.0F );
}
else
{
alpha->real = norm_x * ( -chi_1->real / abs_chi_1 );
alpha->imag = norm_x * ( -chi_1->imag / abs_chi_1 );
}
chi_1_minus_alpha->real = chi_1->real - alpha->real;
chi_1_minus_alpha->imag = chi_1->imag - alpha->imag;
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_csabsval2( chi_1_minus_alpha, &abs_chi_1_minus_alpha );
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
tau->real = one_half.real + one_half.real*( norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha );
tau->imag = 0.0F;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_l_opd | ( | int | m_x2, |
| double * | chi_1, | ||
| double * | x2, | ||
| int | inc_x2, | ||
| double * | alpha, | ||
| double * | chi_1_minus_alpha, | ||
| double * | tau | ||
| ) |
References bl1_dnrm2(), and FLA_ONE_HALF.
Referenced by FLA_Househ2s_UT(), and FLA_Househ2s_UT_r_opd().
{
double one_half = *FLA_DOUBLE_PTR( FLA_ONE_HALF );
double y[2];
double abs_chi_1;
double norm_x_2;
double norm_x;
double abs_chi_1_minus_alpha;
double norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_dnrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0 )
{
*alpha = -(*chi_1);
*chi_1_minus_alpha = 2.0 * (*chi_1);
*tau = one_half;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_dnrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0] = abs_chi_1;
y[1] = norm_x_2;
bl1_dnrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
// = -sign( chi_1 ) * || x ||_2
//
*alpha = -dsign( *chi_1 ) * norm_x;
*chi_1_minus_alpha = (*chi_1) - (*alpha);
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_dabsval2( chi_1_minus_alpha, &abs_chi_1_minus_alpha );
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
*tau = one_half + one_half*( norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_l_ops | ( | int | m_x2, |
| float * | chi_1, | ||
| float * | x2, | ||
| int | inc_x2, | ||
| float * | alpha, | ||
| float * | chi_1_minus_alpha, | ||
| float * | tau | ||
| ) |
References bl1_snrm2(), and FLA_ONE_HALF.
Referenced by FLA_Househ2s_UT(), and FLA_Househ2s_UT_r_ops().
{
float one_half = *FLA_FLOAT_PTR( FLA_ONE_HALF );
float y[2];
float abs_chi_1;
float norm_x_2;
float norm_x;
float abs_chi_1_minus_alpha;
float norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_snrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0F )
{
*alpha = -(*chi_1);
*chi_1_minus_alpha = 2.0F * (*chi_1);
*tau = one_half;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_snrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0] = abs_chi_1;
y[1] = norm_x_2;
bl1_snrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
// = -sign( chi_1 ) * || x ||_2
//
*alpha = -ssign( *chi_1 ) * norm_x;
*chi_1_minus_alpha = (*chi_1) - (*alpha);
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_sabsval2( chi_1_minus_alpha, &abs_chi_1_minus_alpha );
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
*tau = one_half + one_half*( norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_l_opz | ( | int | m_x2, |
| dcomplex * | chi_1, | ||
| dcomplex * | x2, | ||
| int | inc_x2, | ||
| dcomplex * | alpha, | ||
| dcomplex * | chi_1_minus_alpha, | ||
| dcomplex * | tau | ||
| ) |
References bl1_znrm2(), FLA_ONE_HALF, dcomplex::imag, and dcomplex::real.
Referenced by FLA_Househ2s_UT(), and FLA_Househ2s_UT_r_opz().
{
dcomplex one_half = *FLA_DOUBLE_COMPLEX_PTR( FLA_ONE_HALF );
dcomplex y[2];
double abs_chi_1;
double norm_x_2;
double norm_x;
double abs_chi_1_minus_alpha;
double norm_x_2_div_abs_chi_1_minus_alpha;
int i_one = 1;
int i_two = 2;
//
// Compute the 2-norm of x_2:
//
// norm_x_2 := || x_2 ||_2
//
bl1_znrm2( m_x2,
x2, inc_x2,
&norm_x_2 );
//
// If 2-norm of x_2 is zero, then return with trivial values.
//
if ( norm_x_2 == 0.0 )
{
alpha->real = -(chi_1->real);
alpha->imag = -(chi_1->imag);
chi_1_minus_alpha->real = 2.0 * chi_1->real;
chi_1_minus_alpha->imag = 2.0 * chi_1->imag;
tau->real = one_half.real;
tau->imag = one_half.imag;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_1, which equals the 2-norm
// of chi_1:
//
// abs_chi_1 := | chi_1 | = || chi_1 ||_2
//
bl1_znrm2( i_one,
chi_1, i_one,
&abs_chi_1 );
//
// Compute the 2-norm of x via the two norms previously computed above:
//
// norm_x := || x ||_2 = || / chi_1 \ || = || / || chi_1 ||_2 \ ||
// || \ x_2 / ||_2 || \ || x_2 ||_2 / ||_2
//
y[0].real = abs_chi_1;
y[0].imag = 0.0;
y[1].real = norm_x_2;
y[1].imag = 0.0;
bl1_znrm2( i_two,
y, i_one,
&norm_x );
//
// Compute alpha:
//
// alpha := - || x ||_2 * chi_1 / | chi_1 |
//
if ( abs_chi_1 == 0.0 )
{
alpha->real = norm_x * ( -1.0 );
alpha->imag = norm_x * ( -1.0 );
}
else
{
alpha->real = norm_x * ( -chi_1->real / abs_chi_1 );
alpha->imag = norm_x * ( -chi_1->imag / abs_chi_1 );
}
chi_1_minus_alpha->real = chi_1->real - alpha->real;
chi_1_minus_alpha->imag = chi_1->imag - alpha->imag;
//
// Compute tau:
//
// tau := ( 1 + u_2' * u_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_2' * x_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = 1/2 + ( || x ||_2 / | chi_1 - alpha | )^2
//
bl1_zdabsval2( chi_1_minus_alpha, &abs_chi_1_minus_alpha );
norm_x_2_div_abs_chi_1_minus_alpha = norm_x_2 / abs_chi_1_minus_alpha;
tau->real = one_half.real + one_half.real*( norm_x_2_div_abs_chi_1_minus_alpha *
norm_x_2_div_abs_chi_1_minus_alpha );
tau->imag = 0.0;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_r_opc | ( | int | m_x2, |
| scomplex * | chi_1, | ||
| scomplex * | x2, | ||
| int | inc_x2, | ||
| scomplex * | alpha, | ||
| scomplex * | chi_1_minus_alpha, | ||
| scomplex * | tau | ||
| ) |
References FLA_Househ2s_UT_l_opc().
Referenced by FLA_Bidiag_UT_u_step_ofc_var3(), FLA_Bidiag_UT_u_step_ofc_var4(), FLA_Bidiag_UT_u_step_opc_var3(), FLA_Bidiag_UT_u_step_opc_var4(), and FLA_Househ2s_UT().
{
FLA_Househ2s_UT_l_opc( m_x2,
chi_1,
x2, inc_x2,
alpha,
chi_1_minus_alpha,
tau );
//chi_1_minus_alpha->real = chi_1->real - alpha->real;
//chi_1_minus_alpha->imag = chi_1->imag - -alpha->imag;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_r_opd | ( | int | m_x2, |
| double * | chi_1, | ||
| double * | x2, | ||
| int | inc_x2, | ||
| double * | alpha, | ||
| double * | chi_1_minus_alpha, | ||
| double * | tau | ||
| ) |
References FLA_Househ2s_UT_l_opd().
Referenced by FLA_Bidiag_UT_u_step_ofd_var3(), FLA_Bidiag_UT_u_step_ofd_var4(), FLA_Bidiag_UT_u_step_opd_var3(), FLA_Bidiag_UT_u_step_opd_var4(), and FLA_Househ2s_UT().
{
FLA_Househ2s_UT_l_opd( m_x2,
chi_1,
x2, inc_x2,
alpha,
chi_1_minus_alpha,
tau );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_r_ops | ( | int | m_x2, |
| float * | chi_1, | ||
| float * | x2, | ||
| int | inc_x2, | ||
| float * | alpha, | ||
| float * | chi_1_minus_alpha, | ||
| float * | tau | ||
| ) |
References FLA_Househ2s_UT_l_ops().
Referenced by FLA_Bidiag_UT_u_step_ofs_var3(), FLA_Bidiag_UT_u_step_ofs_var4(), FLA_Bidiag_UT_u_step_ops_var3(), FLA_Bidiag_UT_u_step_ops_var4(), and FLA_Househ2s_UT().
{
FLA_Househ2s_UT_l_ops( m_x2,
chi_1,
x2, inc_x2,
alpha,
chi_1_minus_alpha,
tau );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ2s_UT_r_opz | ( | int | m_x2, |
| dcomplex * | chi_1, | ||
| dcomplex * | x2, | ||
| int | inc_x2, | ||
| dcomplex * | alpha, | ||
| dcomplex * | chi_1_minus_alpha, | ||
| dcomplex * | tau | ||
| ) |
References FLA_Househ2s_UT_l_opz().
Referenced by FLA_Bidiag_UT_u_step_ofz_var3(), FLA_Bidiag_UT_u_step_ofz_var4(), FLA_Bidiag_UT_u_step_opz_var3(), FLA_Bidiag_UT_u_step_opz_var4(), and FLA_Househ2s_UT().
{
FLA_Househ2s_UT_l_opz( m_x2,
chi_1,
x2, inc_x2,
alpha,
chi_1_minus_alpha,
tau );
//chi_1_minus_alpha->real = chi_1->real - alpha->real;
//chi_1_minus_alpha->imag = chi_1->imag - -alpha->imag;
return FLA_SUCCESS;
}
References FLA_Check_error_level(), FLA_Househ3UD_UT_check(), FLA_Househ3UD_UT_opc(), FLA_Househ3UD_UT_opd(), FLA_Househ3UD_UT_ops(), FLA_Househ3UD_UT_opz(), FLA_Obj_datatype(), FLA_Obj_vector_dim(), and FLA_Obj_vector_inc().
Referenced by FLA_UDdate_UT_unb_var1().
{
FLA_Datatype datatype;
int m_x1;
int m_y2;
int inc_x1;
int inc_y2;
datatype = FLA_Obj_datatype( x1 );
m_x1 = FLA_Obj_vector_dim( x1 );
m_y2 = FLA_Obj_vector_dim( y2 );
inc_x1 = FLA_Obj_vector_inc( x1 );
inc_y2 = FLA_Obj_vector_inc( y2 );
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Househ3UD_UT_check( chi_0, x1, y2, tau );
switch ( datatype )
{
case FLA_FLOAT:
{
float* chi_0_p = ( float* ) FLA_FLOAT_PTR( chi_0 );
float* x1_p = ( float* ) FLA_FLOAT_PTR( x1 );
float* y2_p = ( float* ) FLA_FLOAT_PTR( y2 );
float* tau_p = ( float* ) FLA_FLOAT_PTR( tau );
FLA_Househ3UD_UT_ops( m_x1,
m_y2,
chi_0_p,
x1_p, inc_x1,
y2_p, inc_y2,
tau_p );
break;
}
case FLA_DOUBLE:
{
double* chi_0_p = ( double* ) FLA_DOUBLE_PTR( chi_0 );
double* x1_p = ( double* ) FLA_DOUBLE_PTR( x1 );
double* y2_p = ( double* ) FLA_DOUBLE_PTR( y2 );
double* tau_p = ( double* ) FLA_DOUBLE_PTR( tau );
FLA_Househ3UD_UT_opd( m_x1,
m_y2,
chi_0_p,
x1_p, inc_x1,
y2_p, inc_y2,
tau_p );
break;
}
case FLA_COMPLEX:
{
scomplex* chi_0_p = ( scomplex* ) FLA_COMPLEX_PTR( chi_0 );
scomplex* x1_p = ( scomplex* ) FLA_COMPLEX_PTR( x1 );
scomplex* y2_p = ( scomplex* ) FLA_COMPLEX_PTR( y2 );
scomplex* tau_p = ( scomplex* ) FLA_COMPLEX_PTR( tau );
FLA_Househ3UD_UT_opc( m_x1,
m_y2,
chi_0_p,
x1_p, inc_x1,
y2_p, inc_y2,
tau_p );
break;
}
case FLA_DOUBLE_COMPLEX:
{
dcomplex* chi_0_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( chi_0 );
dcomplex* x1_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( x1 );
dcomplex* y2_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( y2 );
dcomplex* tau_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( tau );
FLA_Househ3UD_UT_opz( m_x1,
m_y2,
chi_0_p,
x1_p, inc_x1,
y2_p, inc_y2,
tau_p );
break;
}
}
return FLA_SUCCESS;
}
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), and FLA_Check_nonconstant_object().
Referenced by FLA_Househ3UD_UT().
{
FLA_Error e_val;
e_val = FLA_Check_nonconstant_object( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, x2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, y2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( chi_1, tau );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( chi_1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( tau );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ3UD_UT_opc | ( | int | m_x2, |
| int | m_y2, | ||
| scomplex * | chi_1, | ||
| scomplex * | x2, | ||
| int | inc_x2, | ||
| scomplex * | y2, | ||
| int | inc_y2, | ||
| scomplex * | tau | ||
| ) |
References bl1_cinvscalv(), bl1_cnrm2(), BLIS1_NO_CONJUGATE, FLA_ONE_HALF, scomplex::imag, and scomplex::real.
Referenced by FLA_Househ3UD_UT(), and FLA_UDdate_UT_opc_var1().
{
scomplex one_half = *FLA_COMPLEX_PTR( FLA_ONE_HALF );
scomplex alpha;
scomplex chi_0_minus_alpha;
scomplex neg_chi_0_minus_alpha;
float abs_chi_0;
float norm_x_1;
float norm_y_2;
float lambda;
float abs_sq_chi_0_minus_alpha;
int i_one = 1;
//
// Compute the 2-norms of x_1 and y_2:
//
// norm_x_1 := || x_1 ||_2
// norm_y_2 := || y_2 ||_2
//
bl1_cnrm2( m_x1,
x1, inc_x1,
&norm_x_1 );
bl1_cnrm2( m_y2,
y2, inc_y2,
&norm_y_2 );
//
// If 2-norms of x_1, y_2 are zero, then return with trivial tau, chi_0 values.
//
if ( norm_x_1 == 0.0F &&
norm_y_2 == 0.0F )
{
chi_0->real = -(chi_0->real);
chi_0->imag = -(chi_0->imag);
tau->real = one_half.real;
tau->imag = one_half.imag;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_0, which equals the 2-norm
// of chi_0:
//
// abs_chi_0 := | chi_0 | = || chi_0 ||_2
//
bl1_cnrm2( i_one,
chi_0, i_one,
&abs_chi_0 );
//
// Compute lambda:
//
// lambda := sqrt( conj(chi0) chi0 + x1' x1 - y2' y2 )
//
lambda = ( float ) sqrt( abs_chi_0 * abs_chi_0 +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 );
//
// Compute alpha:
//
// alpha := - lambda * chi_0 / | chi_0 |
//
alpha.real = -chi_0->real * lambda / abs_chi_0;
alpha.imag = -chi_0->imag * lambda / abs_chi_0;
//
// Overwrite x_1 and y_2 with u_1 and v_2, respectively:
//
// x_1 := x_1 / ( chi_0 - alpha )
// y_2 := y_2 / -( chi_0 - alpha )
//
chi_0_minus_alpha.real = chi_0->real - alpha.real;
chi_0_minus_alpha.imag = chi_0->imag - alpha.imag;
bl1_cinvscalv( BLIS1_NO_CONJUGATE,
m_x1,
&chi_0_minus_alpha,
x1, inc_x1 );
neg_chi_0_minus_alpha.real = -chi_0_minus_alpha.real;
neg_chi_0_minus_alpha.imag = -chi_0_minus_alpha.imag;
bl1_cinvscalv( BLIS1_NO_CONJUGATE,
m_y2,
&neg_chi_0_minus_alpha,
y2, inc_y2 );
//
// Compute tau:
//
// tau := ( 1 + u_1' * u_1 - v_2' * v_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_1' * x_1 - y_2' * y_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 - || y_2 ||_2^2 ) /
// ( 2 * | chi_1 - alpha |^2 )
//
abs_sq_chi_0_minus_alpha = chi_0_minus_alpha.real * chi_0_minus_alpha.real +
chi_0_minus_alpha.imag * chi_0_minus_alpha.imag;
tau->real = ( abs_sq_chi_0_minus_alpha +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 ) /
( 2.0F * abs_sq_chi_0_minus_alpha );
tau->imag = 0.0F;
//
// Overwrite chi_0 with alpha:
//
// chi_0 := alpha
//
chi_0->real = alpha.real;
chi_0->imag = alpha.imag;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ3UD_UT_opd | ( | int | m_x2, |
| int | m_y2, | ||
| double * | chi_1, | ||
| double * | x2, | ||
| int | inc_x2, | ||
| double * | y2, | ||
| int | inc_y2, | ||
| double * | tau | ||
| ) |
References bl1_dinvscalv(), bl1_dnrm2(), BLIS1_NO_CONJUGATE, and FLA_ONE_HALF.
Referenced by FLA_Househ3UD_UT(), and FLA_UDdate_UT_opd_var1().
{
double one_half = *FLA_DOUBLE_PTR( FLA_ONE_HALF );
double alpha;
double chi_0_minus_alpha;
double neg_chi_0_minus_alpha;
double abs_chi_0;
double norm_x_1;
double norm_y_2;
double lambda;
double abs_sq_chi_0_minus_alpha;
int i_one = 1;
//
// Compute the 2-norms of x_1 and y_2:
//
// norm_x_1 := || x_1 ||_2
// norm_y_2 := || y_2 ||_2
//
bl1_dnrm2( m_x1,
x1, inc_x1,
&norm_x_1 );
bl1_dnrm2( m_y2,
y2, inc_y2,
&norm_y_2 );
//
// If 2-norms of x_1, y_2 are zero, then return with trivial tau, chi_0 values.
//
if ( norm_x_1 == 0.0 &&
norm_y_2 == 0.0 )
{
*chi_0 = -(*chi_0);
*tau = one_half;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_0, which equals the 2-norm
// of chi_0:
//
// abs_chi_0 := | chi_0 | = || chi_0 ||_2
//
bl1_dnrm2( i_one,
chi_0, i_one,
&abs_chi_0 );
//
// Compute lambda:
//
// lambda := sqrt( conj(chi0) chi0 + x1' x1 - y2' y2 )
//
lambda = sqrt( abs_chi_0 * abs_chi_0 +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 );
// Compute alpha:
//
// alpha := - lambda * chi_0 / | chi_0 |
// = -sign( chi_0 ) * lambda
//
alpha = -dsign( *chi_0 ) * lambda;
//
// Overwrite x_1 and y_2 with u_1 and v_2, respectively:
//
// x_1 := x_1 / ( chi_0 - alpha )
// y_2 := y_2 / -( chi_0 - alpha )
//
chi_0_minus_alpha = (*chi_0) - alpha;
bl1_dinvscalv( BLIS1_NO_CONJUGATE,
m_x1,
&chi_0_minus_alpha,
x1, inc_x1 );
neg_chi_0_minus_alpha = -chi_0_minus_alpha;
bl1_dinvscalv( BLIS1_NO_CONJUGATE,
m_y2,
&neg_chi_0_minus_alpha,
y2, inc_y2 );
//
// Compute tau:
//
// tau := ( 1 + u_1' * u_1 - v_2' * v_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_1' * x_1 - y_2' * y_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 - || y_2 ||_2^2 ) /
// ( 2 * | chi_1 - alpha |^2 )
//
abs_sq_chi_0_minus_alpha = chi_0_minus_alpha * chi_0_minus_alpha;
*tau = ( abs_sq_chi_0_minus_alpha +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 ) /
( 2.0 * abs_sq_chi_0_minus_alpha );
//
// Overwrite chi_0 with alpha:
//
// chi_0 := alpha
//
*chi_0 = alpha;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ3UD_UT_ops | ( | int | m_x2, |
| int | m_y2, | ||
| float * | chi_1, | ||
| float * | x2, | ||
| int | inc_x2, | ||
| float * | y2, | ||
| int | inc_y2, | ||
| float * | tau | ||
| ) |
References bl1_sinvscalv(), bl1_snrm2(), BLIS1_NO_CONJUGATE, and FLA_ONE_HALF.
Referenced by FLA_Househ3UD_UT(), and FLA_UDdate_UT_ops_var1().
{
float one_half = *FLA_FLOAT_PTR( FLA_ONE_HALF );
float alpha;
float chi_0_minus_alpha;
float neg_chi_0_minus_alpha;
float abs_chi_0;
float norm_x_1;
float norm_y_2;
float lambda;
float abs_sq_chi_0_minus_alpha;
int i_one = 1;
//
// Compute the 2-norms of x_1 and y_2:
//
// norm_x_1 := || x_1 ||_2
// norm_y_2 := || y_2 ||_2
//
bl1_snrm2( m_x1,
x1, inc_x1,
&norm_x_1 );
bl1_snrm2( m_y2,
y2, inc_y2,
&norm_y_2 );
//
// If 2-norms of x_1, y_2 are zero, then return with trivial tau, chi_0 values.
//
if ( norm_x_1 == 0.0F &&
norm_y_2 == 0.0F )
{
*chi_0 = -(*chi_0);
*tau = one_half;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_0, which equals the 2-norm
// of chi_0:
//
// abs_chi_0 := | chi_0 | = || chi_0 ||_2
//
bl1_snrm2( i_one,
chi_0, i_one,
&abs_chi_0 );
//
// Compute lambda:
//
// lambda := sqrt( conj(chi0) chi0 + x1' x1 - y2' y2 )
//
lambda = ( float ) sqrt( abs_chi_0 * abs_chi_0 +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 );
// Compute alpha:
//
// alpha := - lambda * chi_0 / | chi_0 |
// = -sign( chi_0 ) * lambda
//
alpha = -ssign( *chi_0 ) * lambda;
//
// Overwrite x_1 and y_2 with u_1 and v_2, respectively:
//
// x_1 := x_1 / ( chi_0 - alpha )
// y_2 := y_2 / -( chi_0 - alpha )
//
chi_0_minus_alpha = (*chi_0) - alpha;
bl1_sinvscalv( BLIS1_NO_CONJUGATE,
m_x1,
&chi_0_minus_alpha,
x1, inc_x1 );
neg_chi_0_minus_alpha = -chi_0_minus_alpha;
bl1_sinvscalv( BLIS1_NO_CONJUGATE,
m_y2,
&neg_chi_0_minus_alpha,
y2, inc_y2 );
//
// Compute tau:
//
// tau := ( 1 + u_1' * u_1 - v_2' * v_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_1' * x_1 - y_2' * y_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 - || y_2 ||_2^2 ) /
// ( 2 * | chi_1 - alpha |^2 )
//
abs_sq_chi_0_minus_alpha = chi_0_minus_alpha * chi_0_minus_alpha;
*tau = ( abs_sq_chi_0_minus_alpha +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 ) /
( 2.0F * abs_sq_chi_0_minus_alpha );
//
// Overwrite chi_0 with alpha:
//
// chi_0 := alpha
//
*chi_0 = alpha;
return FLA_SUCCESS;
}
| FLA_Error FLA_Househ3UD_UT_opz | ( | int | m_x2, |
| int | m_y2, | ||
| dcomplex * | chi_1, | ||
| dcomplex * | x2, | ||
| int | inc_x2, | ||
| dcomplex * | y2, | ||
| int | inc_y2, | ||
| dcomplex * | tau | ||
| ) |
References bl1_zinvscalv(), bl1_znrm2(), BLIS1_NO_CONJUGATE, FLA_ONE_HALF, dcomplex::imag, and dcomplex::real.
Referenced by FLA_Househ3UD_UT(), and FLA_UDdate_UT_opz_var1().
{
dcomplex one_half = *FLA_DOUBLE_COMPLEX_PTR( FLA_ONE_HALF );
dcomplex alpha;
dcomplex chi_0_minus_alpha;
dcomplex neg_chi_0_minus_alpha;
double abs_chi_0;
double norm_x_1;
double norm_y_2;
double lambda;
double abs_sq_chi_0_minus_alpha;
int i_one = 1;
//
// Compute the 2-norms of x_1 and y_2:
//
// norm_x_1 := || x_1 ||_2
// norm_y_2 := || y_2 ||_2
//
bl1_znrm2( m_x1,
x1, inc_x1,
&norm_x_1 );
bl1_znrm2( m_y2,
y2, inc_y2,
&norm_y_2 );
//
// If 2-norms of x_1, y_2 are zero, then return with trivial tau, chi_0 values.
//
if ( norm_x_1 == 0.0 &&
norm_y_2 == 0.0 )
{
chi_0->real = -(chi_0->real);
chi_0->imag = -(chi_0->imag);
tau->real = one_half.real;
tau->imag = one_half.imag;
return FLA_SUCCESS;
}
//
// Compute the absolute value (magnitude) of chi_0, which equals the 2-norm
// of chi_0:
//
// abs_chi_0 := | chi_0 | = || chi_0 ||_2
//
bl1_znrm2( i_one,
chi_0, i_one,
&abs_chi_0 );
//
// Compute lambda:
//
// lambda := sqrt( conj(chi0) chi0 + x1' x1 - y2' y2 )
//
lambda = sqrt( abs_chi_0 * abs_chi_0 +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 );
//
// Compute alpha:
//
// alpha := - lambda * chi_0 / | chi_0 |
//
alpha.real = -chi_0->real * lambda / abs_chi_0;
alpha.imag = -chi_0->imag * lambda / abs_chi_0;
//
// Overwrite x_1 and y_2 with u_1 and v_2, respectively:
//
// x_1 := x_1 / ( chi_0 - alpha )
// y_2 := y_2 / -( chi_0 - alpha )
//
chi_0_minus_alpha.real = chi_0->real - alpha.real;
chi_0_minus_alpha.imag = chi_0->imag - alpha.imag;
bl1_zinvscalv( BLIS1_NO_CONJUGATE,
m_x1,
&chi_0_minus_alpha,
x1, inc_x1 );
neg_chi_0_minus_alpha.real = -chi_0_minus_alpha.real;
neg_chi_0_minus_alpha.imag = -chi_0_minus_alpha.imag;
bl1_zinvscalv( BLIS1_NO_CONJUGATE,
m_y2,
&neg_chi_0_minus_alpha,
y2, inc_y2 );
//
// Compute tau:
//
// tau := ( 1 + u_1' * u_1 - v_2' * v_2 ) / 2
// = ( ( chi_1 - alpha ) * conj( chi_1 - alpha ) + x_1' * x_1 - y_2' * y_2 ) /
// ( 2 * ( chi_1 - alpha ) * conj( chi_1 - alpha ) )
// = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 - || y_2 ||_2^2 ) /
// ( 2 * | chi_1 - alpha |^2 )
//
abs_sq_chi_0_minus_alpha = chi_0_minus_alpha.real * chi_0_minus_alpha.real +
chi_0_minus_alpha.imag * chi_0_minus_alpha.imag;
tau->real = ( abs_sq_chi_0_minus_alpha +
norm_x_1 * norm_x_1 -
norm_y_2 * norm_y_2 ) /
( 2.0 * abs_sq_chi_0_minus_alpha );
tau->imag = 0.0;
//
// Overwrite chi_0 with alpha:
//
// chi_0 := alpha
//
chi_0->real = alpha.real;
chi_0->imag = alpha.imag;
return FLA_SUCCESS;
}
| FLA_Error FLA_Introduce_bulge_check | ( | FLA_Obj | shift, |
| FLA_Obj | gamma, | ||
| FLA_Obj | sigma, | ||
| FLA_Obj | delta1, | ||
| FLA_Obj | epsilon1, | ||
| FLA_Obj | delta2, | ||
| FLA_Obj | beta, | ||
| FLA_Obj | epsilon2 | ||
| ) |
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), FLA_Check_nonconstant_object(), and FLA_Check_real_object().
{
FLA_Error e_val;
e_val = FLA_Check_nonconstant_object( delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, shift );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, epsilon1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, delta2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, beta );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, epsilon2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( shift );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( gamma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( sigma );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( epsilon1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( delta2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( beta );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( epsilon2 );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
Referenced by fla_dlamch(), and fla_slamch().
{
/* System generated locals */
logical ret_val;
/* Local variables */
static integer inta, intb, zcode;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* LSAME returns .TRUE. if CA is the same letter as CB regardless of */
/* case. */
/* Arguments */
/* ========= */
/* CA (input) CHARACTER*1 */
/* CB (input) CHARACTER*1 */
/* CA and CB specify the single characters to be compared. */
/* ===================================================================== */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. Executable Statements .. */
/* Test if the characters are equal */
ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
if (ret_val) {
return ret_val;
}
/* Now test for equivalence if both characters are alphabetic. */
zcode = 'Z';
/* Use 'Z' rather than 'A' so that ASCII can be detected on Prime */
/* machines, on which ICHAR returns a value with bit 8 set. */
/* ICHAR('A') on Prime machines returns 193 which is the same as */
/* ICHAR('A') on an EBCDIC machine. */
inta = *(unsigned char *)ca;
intb = *(unsigned char *)cb;
if (zcode == 90 || zcode == 122) {
/* ASCII is assumed - ZCODE is the ASCII code of either lower or */
/* upper case 'Z'. */
if (inta >= 97 && inta <= 122) {
inta += -32;
}
if (intb >= 97 && intb <= 122) {
intb += -32;
}
} else if (zcode == 233 || zcode == 169) {
/* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or */
/* upper case 'Z'. */
if ((inta >= 129 && inta <= 137) || (inta >= 145 && inta <= 153) || (inta
>= 162 && inta <= 169)) {
inta += 64;
}
if ((intb >= 129 && intb <= 137) || (intb >= 145 && intb <= 153) || (intb
>= 162 && intb <= 169)) {
intb += 64;
}
} else if (zcode == 218 || zcode == 250) {
/* ASCII is assumed, on Prime machines - ZCODE is the ASCII code */
/* plus 128 of either lower or upper case 'Z'. */
if (inta >= 225 && inta <= 250) {
inta += -32;
}
if (intb >= 225 && intb <= 250) {
intb += -32;
}
}
ret_val = inta == intb;
/* RETURN */
/* End of LSAME */
return ret_val;
} /* fla_lsame */
References FLA_Check_error_level(), FLA_Cont_with_3x3_to_2x2(), FLA_LU_find_zero_on_diagonal_check(), FLA_Obj_equals(), FLA_Obj_length(), FLA_Obj_min_dim(), FLA_Part_2x2(), FLA_Repart_2x2_to_3x3(), and FLA_ZERO.
Referenced by FLA_LU_nopiv(), and FLASH_LU_find_zero_on_diagonal().
{
FLA_Obj ATL, ATR, A00, a01, A02,
ABL, ABR, a10t, alpha11, a12t,
A20, a21, A22;
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_LU_find_zero_on_diagonal_check( A );
FLA_Part_2x2( A, &ATL, &ATR,
&ABL, &ABR, 0, 0, FLA_TL );
while ( FLA_Obj_length( ATL ) < FLA_Obj_min_dim( A ) ){
FLA_Repart_2x2_to_3x3( ATL, /**/ ATR, &A00, /**/ &a01, &A02,
/* ************* */ /* ************************** */
&a10t, /**/ &alpha11, &a12t,
ABL, /**/ ABR, &A20, /**/ &a21, &A22,
1, 1, FLA_BR );
/*------------------------------------------------------------*/
if ( FLA_Obj_equals( alpha11, FLA_ZERO ) ) return FLA_Obj_length( A00 );
/*------------------------------------------------------------*/
FLA_Cont_with_3x3_to_2x2( &ATL, /**/ &ATR, A00, a01, /**/ A02,
a10t, alpha11, /**/ a12t,
/* ************** */ /* ************************ */
&ABL, /**/ &ABR, A20, a21, /**/ A22,
FLA_TL );
}
return FLA_SUCCESS;
}
References FLA_Check_floating_object(), FLA_Check_nonconstant_object(), and FLA_Check_object_scalar_elemtype().
Referenced by FLA_LU_find_zero_on_diagonal().
{
FLA_Error e_val;
e_val = FLA_Check_floating_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( A );
FLA_Check_error_code( e_val );
e_val = FLA_Check_object_scalar_elemtype( A );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Mach_params | ( | FLA_Machval | machval, |
| FLA_Obj | val | ||
| ) |
References FLA_Check_error_level(), FLA_Mach_params_check(), FLA_Mach_params_opd(), FLA_Mach_params_ops(), and FLA_Obj_datatype().
Referenced by FLA_Hevd_compute_scaling(), FLA_Hevdr_external(), and FLA_Svd_compute_scaling().
{
FLA_Datatype datatype;
datatype = FLA_Obj_datatype( val );
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Mach_params_check( machval, val );
switch ( datatype )
{
case FLA_FLOAT:
{
float* val_p = ( float* ) FLA_FLOAT_PTR( val );
*val_p = FLA_Mach_params_ops( machval );
break;
}
case FLA_DOUBLE:
{
double* val_p = ( double* ) FLA_DOUBLE_PTR( val );
*val_p = FLA_Mach_params_opd( machval );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Mach_params_check | ( | FLA_Machval | machval, |
| FLA_Obj | val | ||
| ) |
References FLA_Check_real_object(), and FLA_Check_valid_machval().
Referenced by FLA_Mach_params().
{
FLA_Error e_val;
e_val = FLA_Check_valid_machval( machval );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( val );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| double FLA_Mach_params_opd | ( | FLA_Machval | machval | ) |
References fla_dlamch(), and FLA_Param_map_flame_to_netlib_machval().
Referenced by FLA_Bsvd_compute_shift_opd(), FLA_Bsvd_compute_tol_thresh_opd(), FLA_Bsvd_ext_opd_var1(), FLA_Bsvd_ext_opz_var1(), FLA_Bsvd_v_opd_var1(), FLA_Bsvd_v_opd_var2(), FLA_Bsvd_v_opz_var1(), FLA_Bsvd_v_opz_var2(), FLA_Givens2_opd(), FLA_Mach_params(), FLA_Svv_2x2_opd(), FLA_Tevd_compute_scaling_opd(), FLA_Tevd_eigval_n_opd_var1(), FLA_Tevd_eigval_v_opd_var1(), FLA_Tevd_eigval_v_opd_var3(), FLA_Tevd_find_submatrix_opd(), FLA_Tevd_francis_n_opd_var1(), FLA_Tevd_francis_v_opd_var1(), and FLA_Tevd_n_opz_var1().
{
static int first_time = TRUE;
static double vals[FLA_MACH_N_VALS];
if ( first_time )
{
char lapack_machval;
int i;
for( i = 0; i < FLA_MACH_N_VALS - 1; ++i )
{
FLA_Param_map_flame_to_netlib_machval( FLA_MACH_START + i, &lapack_machval );
//printf( "querying %d %c\n", FLA_MACH_START + i, lapack_machval );
vals[i] = fla_dlamch( &lapack_machval, 1 );
//printf( "got back %34.29e\n", vals[i] );
}
// Store epsilon^2 in the last element.
vals[i] = vals[0] * vals[0];
first_time = FALSE;
}
return vals[ machval - FLA_MACH_START ];
}
| float FLA_Mach_params_ops | ( | FLA_Machval | machval | ) |
References FLA_Param_map_flame_to_netlib_machval(), and fla_slamch().
Referenced by FLA_Bsvd_compute_shift_ops(), FLA_Bsvd_compute_tol_thresh_ops(), FLA_Bsvd_ext_opc_var1(), FLA_Bsvd_ext_ops_var1(), FLA_Bsvd_v_opc_var1(), FLA_Bsvd_v_ops_var1(), FLA_Mach_params(), FLA_Svv_2x2_ops(), and FLA_Tevd_compute_scaling_ops().
{
static int first_time = TRUE;
static float vals[FLA_MACH_N_VALS];
if ( first_time )
{
char lapack_machval;
int i;
for( i = 0; i < FLA_MACH_N_VALS - 1; ++i )
{
FLA_Param_map_flame_to_netlib_machval( FLA_MACH_START + i, &lapack_machval );
//printf( "querying %d %c\n", FLA_MACH_START + i, lapack_machval );
vals[i] = fla_slamch( &lapack_machval, 1 );
//printf( "got back %34.29e\n", vals[i] );
}
// Store epsilon^2 in the last element.
vals[i] = vals[0] * vals[0];
first_time = FALSE;
}
return vals[ machval - FLA_MACH_START ];
}
| double fla_pow_di | ( | doublereal * | a, |
| integer * | n | ||
| ) |
Referenced by fla_dlamc2(), and fla_dlamch().
{
double pow, x;
integer n;
unsigned long u;
pow = 1;
x = *ap;
n = *bp;
if( n != 0 )
{
if( n < 0 )
{
n = -n;
x = 1/x;
}
for( u = n; ; )
{
if( u & 01 )
pow *= x;
if( u >>= 1 )
x *= x;
else
break;
}
}
return pow;
}
| real fla_pow_ri | ( | real * | a, |
| integer * | n | ||
| ) |
Referenced by fla_slamc2(), and fla_slamch().
{
double pow, x;
integer n;
unsigned long u;
pow = 1;
x = *ap;
n = *bp;
if( n != 0 )
{
if( n < 0 )
{
n = -n;
x = 1/x;
}
for( u = n; ; )
{
if( u & 01 )
pow *= x;
if( u >>= 1 )
x *= x;
else
break;
}
}
return pow;
}
| FLA_Error FLA_Pythag2 | ( | FLA_Obj | chi, |
| FLA_Obj | psi, | ||
| FLA_Obj | rho | ||
| ) |
References FLA_Obj_datatype(), FLA_Pythag2_opd(), and FLA_Pythag2_ops().
{
FLA_Datatype datatype;
datatype = FLA_Obj_datatype( chi );
switch ( datatype )
{
case FLA_FLOAT:
{
float* buff_chi = FLA_FLOAT_PTR( chi );
float* buff_psi = FLA_FLOAT_PTR( psi );
float* buff_rho = FLA_FLOAT_PTR( rho );
FLA_Pythag2_ops( buff_chi,
buff_psi,
buff_rho );
break;
}
case FLA_DOUBLE:
{
double* buff_chi = FLA_DOUBLE_PTR( chi );
double* buff_psi = FLA_DOUBLE_PTR( psi );
double* buff_rho = FLA_DOUBLE_PTR( rho );
FLA_Pythag2_opd( buff_chi,
buff_psi,
buff_rho );
break;
}
case FLA_COMPLEX:
{
FLA_Check_error_code( FLA_OBJECT_NOT_REAL );
break;
}
case FLA_DOUBLE_COMPLEX:
{
FLA_Check_error_code( FLA_OBJECT_NOT_REAL );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Pythag2_opd | ( | double * | chi, |
| double * | psi, | ||
| double * | rho | ||
| ) |
References bl1_d0(), and bl1_d1().
Referenced by FLA_Pythag2().
| FLA_Error FLA_Pythag2_ops | ( | float * | chi, |
| float * | psi, | ||
| float * | rho | ||
| ) |
References bl1_s0(), and bl1_s1().
Referenced by FLA_Pythag2().
References FLA_Obj_datatype(), FLA_Pythag3_opd(), and FLA_Pythag3_ops().
{
FLA_Datatype datatype;
datatype = FLA_Obj_datatype( chi );
switch ( datatype )
{
case FLA_FLOAT:
{
float* buff_chi = FLA_FLOAT_PTR( chi );
float* buff_psi = FLA_FLOAT_PTR( psi );
float* buff_zeta = FLA_FLOAT_PTR( zeta );
float* buff_rho = FLA_FLOAT_PTR( rho );
FLA_Pythag3_ops( buff_chi,
buff_psi,
buff_zeta,
buff_rho );
break;
}
case FLA_DOUBLE:
{
double* buff_chi = FLA_DOUBLE_PTR( chi );
double* buff_psi = FLA_DOUBLE_PTR( psi );
double* buff_zeta = FLA_DOUBLE_PTR( zeta );
double* buff_rho = FLA_DOUBLE_PTR( rho );
FLA_Pythag3_opd( buff_chi,
buff_psi,
buff_zeta,
buff_rho );
break;
}
case FLA_COMPLEX:
{
FLA_Check_error_code( FLA_OBJECT_NOT_REAL );
break;
}
case FLA_DOUBLE_COMPLEX:
{
FLA_Check_error_code( FLA_OBJECT_NOT_REAL );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Pythag3_opd | ( | double * | chi, |
| double * | psi, | ||
| double * | zeta, | ||
| double * | rho | ||
| ) |
References bl1_d0().
Referenced by FLA_Pythag3().
{
double zero = bl1_d0();
double xabs, yabs, zabs;
double w;
double xabsdivw;
double yabsdivw;
double zabsdivw;
xabs = fabs( *chi );
yabs = fabs( *psi );
zabs = fabs( *zeta );
w = max( xabs, max( yabs, zabs ) );
if ( w == zero )
{
// From netlib dlapy3:
// W can be zero for max(0,nan,0). Adding all three entries
// together will make sure NaN will not disappear.
*rho = xabs + yabs + zabs;
}
else
{
xabsdivw = xabs / w;
yabsdivw = yabs / w;
zabsdivw = zabs / w;
*rho = w * sqrt( xabsdivw * xabsdivw +
yabsdivw * yabsdivw +
zabsdivw * zabsdivw );
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Pythag3_ops | ( | float * | chi, |
| float * | psi, | ||
| float * | zeta, | ||
| float * | rho | ||
| ) |
References bl1_s0().
Referenced by FLA_Pythag3().
{
float zero = bl1_s0();
float xabs, yabs, zabs;
float w;
float xabsdivw;
float yabsdivw;
float zabsdivw;
xabs = fabsf( *chi );
yabs = fabsf( *psi );
zabs = fabsf( *zeta );
w = max( xabs, max( yabs, zabs ) );
if ( w == zero )
{
// From netlib dlapy3:
// W can be zero for max(0,nan,0). Adding all three entries
// together will make sure NaN will not disappear.
*rho = xabs + yabs + zabs;
}
else
{
xabsdivw = xabs / w;
yabsdivw = yabs / w;
zabsdivw = zabs / w;
*rho = w * sqrt( xabsdivw * xabsdivw +
yabsdivw * yabsdivw +
zabsdivw * zabsdivw );
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Shift_pivots_to | ( | FLA_Pivot_type | ptype, |
| FLA_Obj | p | ||
| ) |
References FLA_Check_error_level(), FLA_Obj_length(), FLA_Obj_width(), and FLA_Shift_pivots_to_check().
Referenced by FLA_LU_piv_blk_external(), and FLA_LU_piv_unb_external().
{
int m_p, n_p;
int* buff_p;
int i;
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Shift_pivots_to_check( ptype, p );
m_p = FLA_Obj_length( p );
n_p = FLA_Obj_width( p );
buff_p = FLA_INT_PTR( p );
if ( m_p < 1 || n_p < 1 ) return FLA_SUCCESS;
if ( ptype == FLA_LAPACK_PIVOTS )
{
// Shift FLAME pivots to LAPACK pivots.
for ( i = 0; i < m_p; i++ )
buff_p[ i ] += i + 1;
}
else
{
// Otherwise, shift LAPACK pivots back to FLAME.
for ( i = 0; i < m_p; i++ )
buff_p[ i ] -= i + 1;
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Shift_pivots_to_check | ( | FLA_Pivot_type | ptype, |
| FLA_Obj | p | ||
| ) |
References FLA_Check_col_vector(), FLA_Check_int_object(), FLA_Check_nonconstant_object(), and FLA_Check_valid_pivot_type().
Referenced by FLA_Shift_pivots_to().
{
FLA_Error e_val;
e_val = FLA_Check_valid_pivot_type( ptype );
FLA_Check_error_code( e_val );
e_val = FLA_Check_int_object( p );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( p );
FLA_Check_error_code( e_val );
e_val = FLA_Check_col_vector( p );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| real fla_slamch | ( | char * | cmach, |
| ftnlen | cmach_len | ||
| ) |
References fla_lsame(), fla_pow_ri(), and fla_slamc2().
Referenced by FLA_Mach_params_ops().
{
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
integer i__1;
real ret_val;
/* Builtin functions */
double fla_pow_ri(real *, integer *);
/* Local variables */
static real base;
static integer beta;
static real emin, prec, emax;
static integer imin, imax;
static logical lrnd;
static real rmin, rmax, t, rmach;
extern logical fla_lsame(char *, char *, ftnlen, ftnlen);
static real small, sfmin;
extern /* Subroutine */ int fla_slamc2(integer *, integer *, logical *, real
*, integer *, real *, integer *, real *);
static integer it;
static real rnd, eps;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* SLAMCH determines single precision machine parameters. */
/* Arguments */
/* ========= */
/* CMACH (input) CHARACTER*1 */
/* Specifies the value to be returned by SLAMCH: */
/* = 'E' or 'e', SLAMCH := eps */
/* = 'S' or 's , SLAMCH := sfmin */
/* = 'B' or 'b', SLAMCH := base */
/* = 'P' or 'p', SLAMCH := eps*base */
/* = 'N' or 'n', SLAMCH := t */
/* = 'R' or 'r', SLAMCH := rnd */
/* = 'M' or 'm', SLAMCH := emin */
/* = 'U' or 'u', SLAMCH := rmin */
/* = 'L' or 'l', SLAMCH := emax */
/* = 'O' or 'o', SLAMCH := rmax */
/* where */
/* eps = relative machine precision */
/* sfmin = safe minimum, such that 1/sfmin does not overflow */
/* base = base of the machine */
/* prec = eps*base */
/* t = number of (base) digits in the mantissa */
/* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
/* emin = minimum exponent before (gradual) underflow */
/* rmin = underflow threshold - base**(emin-1) */
/* emax = largest exponent before overflow */
/* rmax = overflow threshold - (base**emax)*(1-eps) */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Save statement .. */
/* .. */
/* .. Data statements .. */
/* .. */
/* .. Executable Statements .. */
if (first) {
fla_slamc2(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
base = (real) beta;
t = (real) it;
if (lrnd) {
rnd = (float)1.;
i__1 = 1 - it;
eps = fla_pow_ri(&base, &i__1) / 2;
} else {
rnd = (float)0.;
i__1 = 1 - it;
eps = fla_pow_ri(&base, &i__1);
}
prec = eps * base;
emin = (real) imin;
emax = (real) imax;
sfmin = rmin;
small = (float)1. / rmax;
if (small >= sfmin) {
/* Use SMALL plus a bit, to avoid the possibility of rounding */
/* causing overflow when computing 1/sfmin. */
sfmin = small * (eps + (float)1.);
}
}
if (fla_lsame(cmach, "E", (ftnlen)1, (ftnlen)1)) {
rmach = eps;
} else if (fla_lsame(cmach, "S", (ftnlen)1, (ftnlen)1)) {
rmach = sfmin;
} else if (fla_lsame(cmach, "B", (ftnlen)1, (ftnlen)1)) {
rmach = base;
} else if (fla_lsame(cmach, "P", (ftnlen)1, (ftnlen)1)) {
rmach = prec;
} else if (fla_lsame(cmach, "N", (ftnlen)1, (ftnlen)1)) {
rmach = t;
} else if (fla_lsame(cmach, "R", (ftnlen)1, (ftnlen)1)) {
rmach = rnd;
} else if (fla_lsame(cmach, "M", (ftnlen)1, (ftnlen)1)) {
rmach = emin;
} else if (fla_lsame(cmach, "U", (ftnlen)1, (ftnlen)1)) {
rmach = rmin;
} else if (fla_lsame(cmach, "L", (ftnlen)1, (ftnlen)1)) {
rmach = emax;
} else if (fla_lsame(cmach, "O", (ftnlen)1, (ftnlen)1)) {
rmach = rmax;
}
ret_val = rmach;
first = FALSE_;
return ret_val;
/* End of SLAMCH */
} /* fla_slamch_ */
| FLA_Error FLA_Sort_bsvd_ext | ( | FLA_Direct | direct, |
| FLA_Obj | s, | ||
| FLA_Bool | apply_U, | ||
| FLA_Obj | U, | ||
| FLA_Bool | apply_V, | ||
| FLA_Obj | V, | ||
| FLA_Bool | apply_C, | ||
| FLA_Obj | C | ||
| ) |
References FLA_Obj_length(), FLA_Obj_vector_dim(), FLA_Obj_vector_inc(), FLA_Obj_width(), FLA_Sort(), FLA_Sort_bsvd_ext_b_opc(), FLA_Sort_bsvd_ext_b_opd(), FLA_Sort_bsvd_ext_b_ops(), FLA_Sort_bsvd_ext_b_opz(), FLA_Sort_bsvd_ext_f_opc(), FLA_Sort_bsvd_ext_f_opd(), FLA_Sort_bsvd_ext_f_ops(), and FLA_Sort_bsvd_ext_f_opz().
{
FLA_Datatype datatype;
dim_t m_U, rs_U, cs_U;
dim_t m_V, rs_V, cs_V;
dim_t n_C, rs_C, cs_C;
dim_t m_s, inc_s;
//if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
// FLA_Sort_bsvd_check( direct, s,
// apply_U, U,
// apply_V, V,
// apply_C, C );
// Sort singular values only; quick sort
if ( apply_U == FALSE && apply_V == FALSE )
return FLA_Sort( direct, s );
// s dimensions must be provided.
m_s = FLA_Obj_vector_dim( s );
inc_s = FLA_Obj_vector_inc( s );
// Datatype of U, V and C must be consistent and must be defined from one of them.
FLA_SORT_BSVD_EXT_DEFINE_OBJ_VARIABLES( U, apply_U, datatype, m_U, FLA_Obj_length, rs_U, cs_U );
FLA_SORT_BSVD_EXT_DEFINE_OBJ_VARIABLES( V, apply_V, datatype, m_V, FLA_Obj_length, rs_V, cs_V );
FLA_SORT_BSVD_EXT_DEFINE_OBJ_VARIABLES( C, apply_C, datatype, n_C, FLA_Obj_width, rs_C, cs_C );
switch ( datatype )
{
case FLA_FLOAT:
{
float* s_p = ( float* ) FLA_FLOAT_PTR( s );
float* U_p = ( apply_U == TRUE ? ( float* ) FLA_FLOAT_PTR( U ) : NULL );
float* V_p = ( apply_V == TRUE ? ( float* ) FLA_FLOAT_PTR( V ) : NULL );
float* C_p = ( apply_C == TRUE ? ( float* ) FLA_FLOAT_PTR( C ) : NULL );
if ( direct == FLA_FORWARD )
FLA_Sort_bsvd_ext_f_ops( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_bsvd_ext_b_ops( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
break;
}
case FLA_DOUBLE:
{
double* s_p = ( double* ) FLA_DOUBLE_PTR( s );
double* U_p = ( apply_U == TRUE ? ( double* ) FLA_DOUBLE_PTR( U ) : NULL );
double* V_p = ( apply_V == TRUE ? ( double* ) FLA_DOUBLE_PTR( V ) : NULL );
double* C_p = ( apply_C == TRUE ? ( double* ) FLA_DOUBLE_PTR( C ) : NULL );
if ( direct == FLA_FORWARD )
FLA_Sort_bsvd_ext_f_opd( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_bsvd_ext_b_opd( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
break;
}
case FLA_COMPLEX:
{
float* s_p = ( float* ) FLA_FLOAT_PTR( s );
scomplex* U_p = ( apply_U == TRUE ? ( scomplex* ) FLA_COMPLEX_PTR( U ) : NULL );
scomplex* V_p = ( apply_V == TRUE ? ( scomplex* ) FLA_COMPLEX_PTR( V ) : NULL );
scomplex* C_p = ( apply_C == TRUE ? ( scomplex* ) FLA_COMPLEX_PTR( C ) : NULL );
if ( direct == FLA_FORWARD )
FLA_Sort_bsvd_ext_f_opc( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_bsvd_ext_b_opc( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
break;
}
case FLA_DOUBLE_COMPLEX:
{
double* s_p = ( double* ) FLA_DOUBLE_PTR( s );
dcomplex* U_p = ( apply_U == TRUE ? ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( U ) : NULL );
dcomplex* V_p = ( apply_V == TRUE ? ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( V ) : NULL );
dcomplex* C_p = ( apply_C == TRUE ? ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( C ) : NULL );
if ( direct == FLA_FORWARD )
FLA_Sort_bsvd_ext_f_opz( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_bsvd_ext_b_opz( m_s, s_p, inc_s,
m_U, U_p, rs_U, cs_U,
m_V, V_p, rs_V, cs_V,
n_C, C_p, rs_C, cs_C );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_b_opc | ( | int | m_s, |
| float * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| scomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| scomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| scomplex * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_cswapv().
Referenced by FLA_Bsvd_ext_opt_var1(), and FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
float p;
FLA_SORT_BSVD_EXT_BODY( BACKWARD, bl1_cswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_b_opd | ( | int | m_s, |
| double * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| double * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| double * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| double * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_dswapv().
Referenced by FLA_Bsvd_ext_opt_var1(), and FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
double p;
FLA_SORT_BSVD_EXT_BODY( BACKWARD, bl1_dswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_b_ops | ( | int | m_s, |
| float * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| float * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| float * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| float * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_sswapv().
Referenced by FLA_Bsvd_ext_opt_var1(), and FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
float p;
FLA_SORT_BSVD_EXT_BODY( BACKWARD, bl1_sswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_b_opz | ( | int | m_s, |
| double * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| dcomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| dcomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| dcomplex * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_zswapv().
Referenced by FLA_Bsvd_ext_opt_var1(), and FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
double p;
FLA_SORT_BSVD_EXT_BODY( BACKWARD, bl1_zswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_f_opc | ( | int | m_s, |
| float * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| scomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| scomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| scomplex * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_cswapv().
Referenced by FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
float p;
FLA_SORT_BSVD_EXT_BODY( FORWARD, bl1_cswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_f_opd | ( | int | m_s, |
| double * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| double * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| double * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| double * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_dswapv().
Referenced by FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
float p;
FLA_SORT_BSVD_EXT_BODY( FORWARD, bl1_dswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_f_ops | ( | int | m_s, |
| float * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| float * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| float * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| float * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_sswapv().
Referenced by FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
float p;
FLA_SORT_BSVD_EXT_BODY( FORWARD, bl1_sswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_bsvd_ext_f_opz | ( | int | m_s, |
| double * | s, | ||
| int | inc_s, | ||
| int | m_U, | ||
| dcomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| int | m_V, | ||
| dcomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V, | ||
| int | n_C, | ||
| dcomplex * | C, | ||
| int | rs_C, | ||
| int | cs_C | ||
| ) |
References bl1_zswapv().
Referenced by FLA_Sort_bsvd_ext().
{
int i, ii, j, k;
double p;
FLA_SORT_BSVD_EXT_BODY( FORWARD, bl1_zswapv );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd | ( | FLA_Direct | direct, |
| FLA_Obj | l, | ||
| FLA_Obj | V | ||
| ) |
References FLA_Check_error_level(), FLA_Obj_col_stride(), FLA_Obj_datatype(), FLA_Obj_length(), FLA_Obj_row_stride(), FLA_Obj_vector_inc(), FLA_Sort_evd_b_opc(), FLA_Sort_evd_b_opd(), FLA_Sort_evd_b_ops(), FLA_Sort_evd_b_opz(), FLA_Sort_evd_check(), FLA_Sort_evd_f_opc(), FLA_Sort_evd_f_opd(), FLA_Sort_evd_f_ops(), and FLA_Sort_evd_f_opz().
Referenced by FLA_Hevd_lv_unb_var1(), and FLA_Hevd_lv_unb_var2().
{
FLA_Datatype datatype;
dim_t m_A;
dim_t rs_V, cs_V;
dim_t inc_l;
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Sort_evd_check( direct, l, V );
datatype = FLA_Obj_datatype( V );
m_A = FLA_Obj_length( V );
rs_V = FLA_Obj_row_stride( V );
cs_V = FLA_Obj_col_stride( V );
inc_l = FLA_Obj_vector_inc( l );
switch ( datatype )
{
case FLA_FLOAT:
{
float* l_p = ( float* ) FLA_FLOAT_PTR( l );
float* V_p = ( float* ) FLA_FLOAT_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_evd_f_ops( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_evd_b_ops( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
break;
}
case FLA_DOUBLE:
{
double* l_p = ( double* ) FLA_DOUBLE_PTR( l );
double* V_p = ( double* ) FLA_DOUBLE_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_evd_f_opd( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_evd_b_opd( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
break;
}
case FLA_COMPLEX:
{
float* l_p = ( float* ) FLA_FLOAT_PTR( l );
scomplex* V_p = ( scomplex* ) FLA_COMPLEX_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_evd_f_opc( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_evd_b_opc( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
break;
}
case FLA_DOUBLE_COMPLEX:
{
double* l_p = ( double* ) FLA_DOUBLE_PTR( l );
dcomplex* V_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_evd_f_opz( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_evd_b_opz( m_A,
l_p, inc_l,
V_p, rs_V, cs_V );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_b_opc | ( | int | m_A, |
| float * | l, | ||
| int | inc_l, | ||
| scomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_evd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_b_opd | ( | int | m_A, |
| double * | l, | ||
| int | inc_l, | ||
| double * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_dswapv().
Referenced by FLA_Sort_evd().
{
int i, ii, j, k;
double p;
for ( ii = 1; ii < m_A; ++ii )
{
i = ii - 1;
k = i;
p = l[ i*inc_l ];
for ( j = ii; j < m_A; ++j )
{
if ( l[ j*inc_l ] > p )
{
k = j;
p = l[ j*inc_l ];
}
}
if ( k != i )
{
l[ k*inc_l ] = l[ i ];
l[ i ] = p;
bl1_dswapv( m_A,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_b_ops | ( | int | m_A, |
| float * | l, | ||
| int | inc_l, | ||
| float * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_evd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_b_opz | ( | int | m_A, |
| double * | l, | ||
| int | inc_l, | ||
| dcomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_zswapv().
Referenced by FLA_Sort_evd().
{
int i, ii, j, k;
double p;
for ( ii = 1; ii < m_A; ++ii )
{
i = ii - 1;
k = i;
p = l[ i*inc_l ];
for ( j = ii; j < m_A; ++j )
{
if ( l[ j*inc_l ] > p )
{
k = j;
p = l[ j*inc_l ];
}
}
if ( k != i )
{
l[ k*inc_l ] = l[ i ];
l[ i ] = p;
bl1_zswapv( m_A,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_check | ( | FLA_Direct | direct, |
| FLA_Obj | l, | ||
| FLA_Obj | V | ||
| ) |
References FLA_Check_floating_object(), FLA_Check_identical_object_precision(), FLA_Check_nonconstant_object(), FLA_Check_object_length_equals(), FLA_Check_real_object(), FLA_Check_valid_direct(), and FLA_Obj_vector_dim().
Referenced by FLA_Sort_evd().
{
FLA_Error e_val;
e_val = FLA_Check_valid_direct( direct );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( l );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( l );
FLA_Check_error_code( e_val );
e_val = FLA_Check_floating_object( V );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( V );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_precision( l, V );
FLA_Check_error_code( e_val );
e_val = FLA_Check_object_length_equals( V, FLA_Obj_vector_dim( l ) );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_f_opc | ( | int | m_A, |
| float * | l, | ||
| int | inc_l, | ||
| scomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_evd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_f_opd | ( | int | m_A, |
| double * | l, | ||
| int | inc_l, | ||
| double * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_dswapv().
Referenced by FLA_Sort_evd().
{
int i, ii, j, k;
double p;
for ( ii = 1; ii < m_A; ++ii )
{
i = ii - 1;
k = i;
p = l[ i*inc_l ];
for ( j = ii; j < m_A; ++j )
{
if ( l[ j*inc_l ] < p )
{
k = j;
p = l[ j*inc_l ];
}
}
if ( k != i )
{
l[ k*inc_l ] = l[ i ];
l[ i ] = p;
bl1_dswapv( m_A,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_f_ops | ( | int | m_A, |
| float * | l, | ||
| int | inc_l, | ||
| float * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_evd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_evd_f_opz | ( | int | m_A, |
| double * | l, | ||
| int | inc_l, | ||
| dcomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_zswapv().
Referenced by FLA_Sort_evd().
{
int i, ii, j, k;
double p;
for ( ii = 1; ii < m_A; ++ii )
{
i = ii - 1;
k = i;
p = l[ i*inc_l ];
for ( j = ii; j < m_A; ++j )
{
if ( l[ j*inc_l ] < p )
{
k = j;
p = l[ j*inc_l ];
}
}
if ( k != i )
{
l[ k*inc_l ] = l[ i ];
l[ i ] = p;
bl1_zswapv( m_A,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd | ( | FLA_Direct | direct, |
| FLA_Obj | s, | ||
| FLA_Obj | U, | ||
| FLA_Obj | V | ||
| ) |
References FLA_Check_error_level(), FLA_Obj_col_stride(), FLA_Obj_datatype(), FLA_Obj_length(), FLA_Obj_row_stride(), FLA_Obj_vector_inc(), FLA_Sort_svd_b_opc(), FLA_Sort_svd_b_opd(), FLA_Sort_svd_b_ops(), FLA_Sort_svd_b_opz(), FLA_Sort_svd_check(), FLA_Sort_svd_f_opc(), FLA_Sort_svd_f_opd(), FLA_Sort_svd_f_ops(), and FLA_Sort_svd_f_opz().
Referenced by FLA_Svd_uv_unb_var1(), and FLA_Svd_uv_unb_var2().
{
FLA_Datatype datatype;
dim_t m_U, n_V;
dim_t rs_U, cs_U;
dim_t rs_V, cs_V;
dim_t inc_s;
if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
FLA_Sort_svd_check( direct, s, U, V );
datatype = FLA_Obj_datatype( U );
m_U = FLA_Obj_length( U );
n_V = FLA_Obj_length( V );
rs_U = FLA_Obj_row_stride( U );
cs_U = FLA_Obj_col_stride( U );
rs_V = FLA_Obj_row_stride( V );
cs_V = FLA_Obj_col_stride( V );
inc_s = FLA_Obj_vector_inc( s );
switch ( datatype )
{
case FLA_FLOAT:
{
float* s_p = ( float* ) FLA_FLOAT_PTR( s );
float* U_p = ( float* ) FLA_FLOAT_PTR( U );
float* V_p = ( float* ) FLA_FLOAT_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_svd_f_ops( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_svd_b_ops( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
break;
}
case FLA_DOUBLE:
{
double* s_p = ( double* ) FLA_DOUBLE_PTR( s );
double* U_p = ( double* ) FLA_DOUBLE_PTR( U );
double* V_p = ( double* ) FLA_DOUBLE_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_svd_f_opd( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_svd_b_opd( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
break;
}
case FLA_COMPLEX:
{
float* s_p = ( float* ) FLA_FLOAT_PTR( s );
scomplex* U_p = ( scomplex* ) FLA_COMPLEX_PTR( U );
scomplex* V_p = ( scomplex* ) FLA_COMPLEX_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_svd_f_opc( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_svd_b_opc( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
break;
}
case FLA_DOUBLE_COMPLEX:
{
double* s_p = ( double* ) FLA_DOUBLE_PTR( s );
dcomplex* U_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( U );
dcomplex* V_p = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( V );
if ( direct == FLA_FORWARD )
FLA_Sort_svd_f_opz( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
else // if ( direct == FLA_BACKWARD )
FLA_Sort_svd_b_opz( m_U,
n_V,
s_p, inc_s,
U_p, rs_U, cs_U,
V_p, rs_V, cs_V );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_b_opc | ( | int | m_U, |
| int | n_V, | ||
| float * | s, | ||
| int | inc_s, | ||
| scomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| scomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_svd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_b_opd | ( | int | m_U, |
| int | n_V, | ||
| double * | s, | ||
| int | inc_s, | ||
| double * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| double * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_dswapv().
Referenced by FLA_Sort_svd().
{
int min_m_n = min( m_U, n_V );
int i, ii, j, k;
double p;
for ( ii = 1; ii < min_m_n; ++ii )
{
i = ii - 1;
k = i;
p = s[ i*inc_s ];
for ( j = ii; j < min_m_n; ++j )
{
if ( s[ j*inc_s ] > p )
{
k = j;
p = s[ j*inc_s ];
}
}
if ( k != i )
{
s[ k*inc_s ] = s[ i ];
s[ i ] = p;
bl1_dswapv( m_U,
U + i*cs_U, rs_U,
U + k*cs_U, rs_U );
bl1_dswapv( n_V,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_b_ops | ( | int | m_U, |
| int | n_V, | ||
| float * | s, | ||
| int | inc_s, | ||
| float * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| float * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_svd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_b_opz | ( | int | m_U, |
| int | n_V, | ||
| double * | s, | ||
| int | inc_s, | ||
| dcomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| dcomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_zswapv().
Referenced by FLA_Sort_svd().
{
int min_m_n = min( m_U, n_V );
int i, ii, j, k;
double p;
for ( ii = 1; ii < min_m_n; ++ii )
{
i = ii - 1;
k = i;
p = s[ i*inc_s ];
for ( j = ii; j < min_m_n; ++j )
{
if ( s[ j*inc_s ] > p )
{
k = j;
p = s[ j*inc_s ];
}
}
if ( k != i )
{
s[ k*inc_s ] = s[ i ];
s[ i ] = p;
bl1_zswapv( m_U,
U + i*cs_U, rs_U,
U + k*cs_U, rs_U );
bl1_zswapv( n_V,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_check | ( | FLA_Direct | direct, |
| FLA_Obj | s, | ||
| FLA_Obj | U, | ||
| FLA_Obj | V | ||
| ) |
References FLA_Check_floating_object(), FLA_Check_identical_object_datatype(), FLA_Check_identical_object_precision(), FLA_Check_nonconstant_object(), FLA_Check_real_object(), FLA_Check_valid_direct(), FLA_Check_vector_dim(), and FLA_Obj_length().
Referenced by FLA_Sort_svd().
{
FLA_Error e_val;
e_val = FLA_Check_valid_direct( direct );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( s );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( s );
FLA_Check_error_code( e_val );
e_val = FLA_Check_floating_object( U );
FLA_Check_error_code( e_val );
e_val = FLA_Check_nonconstant_object( U );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( U, V );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_precision( s, U );
FLA_Check_error_code( e_val );
//e_val = FLA_Check_square( U );
//FLA_Check_error_code( e_val );
//e_val = FLA_Check_square( V );
//FLA_Check_error_code( e_val );
e_val = FLA_Check_vector_dim( s, min( FLA_Obj_length( U ), FLA_Obj_length( V ) ) );
FLA_Check_error_code( e_val );
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_f_opc | ( | int | m_U, |
| int | n_V, | ||
| float * | s, | ||
| int | inc_s, | ||
| scomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| scomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_svd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_f_opd | ( | int | m_U, |
| int | n_V, | ||
| double * | s, | ||
| int | inc_s, | ||
| double * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| double * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_dswapv().
Referenced by FLA_Sort_svd().
{
int min_m_n = min( m_U, n_V );
int i, ii, j, k;
double p;
for ( ii = 1; ii < min_m_n; ++ii )
{
i = ii - 1;
k = i;
p = s[ i*inc_s ];
for ( j = ii; j < min_m_n; ++j )
{
if ( s[ j*inc_s ] < p )
{
k = j;
p = s[ j*inc_s ];
}
}
if ( k != i )
{
s[ k*inc_s ] = s[ i ];
s[ i ] = p;
bl1_dswapv( m_U,
U + i*cs_U, rs_U,
U + k*cs_U, rs_U );
bl1_dswapv( n_V,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_f_ops | ( | int | m_U, |
| int | n_V, | ||
| float * | s, | ||
| int | inc_s, | ||
| float * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| float * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
Referenced by FLA_Sort_svd().
{
return FLA_SUCCESS;
}
| FLA_Error FLA_Sort_svd_f_opz | ( | int | m_U, |
| int | n_V, | ||
| double * | s, | ||
| int | inc_s, | ||
| dcomplex * | U, | ||
| int | rs_U, | ||
| int | cs_U, | ||
| dcomplex * | V, | ||
| int | rs_V, | ||
| int | cs_V | ||
| ) |
References bl1_zswapv().
Referenced by FLA_Sort_svd().
{
int min_m_n = min( m_U, n_V );
int i, ii, j, k;
double p;
for ( ii = 1; ii < min_m_n; ++ii )
{
i = ii - 1;
k = i;
p = s[ i*inc_s ];
for ( j = ii; j < min_m_n; ++j )
{
if ( s[ j*inc_s ] < p )
{
k = j;
p = s[ j*inc_s ];
}
}
if ( k != i )
{
s[ k*inc_s ] = s[ i ];
s[ i ] = p;
bl1_zswapv( m_U,
U + i*cs_U, rs_U,
U + k*cs_U, rs_U );
bl1_zswapv( n_V,
V + i*cs_V, rs_V,
V + k*cs_V, rs_V );
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sv_2x2 | ( | FLA_Obj | alpha11, |
| FLA_Obj | alpha12, | ||
| FLA_Obj | alpha22, | ||
| FLA_Obj | sigma1, | ||
| FLA_Obj | sigma2 | ||
| ) |
References FLA_Obj_datatype(), FLA_Sv_2x2_opd(), and FLA_Sv_2x2_ops().
{
FLA_Datatype datatype;
datatype = FLA_Obj_datatype( alpha11 );
switch ( datatype )
{
case FLA_FLOAT:
{
float* buff_alpha11 = FLA_FLOAT_PTR( alpha11 );
float* buff_alpha12 = FLA_FLOAT_PTR( alpha12 );
float* buff_alpha22 = FLA_FLOAT_PTR( alpha22 );
float* buff_sigma1 = FLA_FLOAT_PTR( sigma1 );
float* buff_sigma2 = FLA_FLOAT_PTR( sigma2 );
FLA_Sv_2x2_ops( buff_alpha11,
buff_alpha12,
buff_alpha22,
buff_sigma1,
buff_sigma2 );
break;
}
case FLA_DOUBLE:
{
double* buff_alpha11 = FLA_DOUBLE_PTR( alpha11 );
double* buff_alpha12 = FLA_DOUBLE_PTR( alpha12 );
double* buff_alpha22 = FLA_DOUBLE_PTR( alpha22 );
double* buff_sigma1 = FLA_DOUBLE_PTR( sigma1 );
double* buff_sigma2 = FLA_DOUBLE_PTR( sigma2 );
FLA_Sv_2x2_opd( buff_alpha11,
buff_alpha12,
buff_alpha22,
buff_sigma1,
buff_sigma2 );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Sv_2x2_opd | ( | double * | alpha11, |
| double * | alpha12, | ||
| double * | alpha22, | ||
| double * | sigma1, | ||
| double * | sigma2 | ||
| ) |
Referenced by FLA_Bsvd_compute_shift_opd(), and FLA_Sv_2x2().
{
double zero = 0.0;
double one = 1.0;
double two = 2.0;
double f, g, h;
double as, at, au, c, fa, fhmin, fhmax, ga, ha;
double ssmin, ssmax;
double temp, temp2;
f = *alpha11;
g = *alpha12;
h = *alpha22;
fa = fabs( f );
ga = fabs( g );
ha = fabs( h );
fhmin = min( fa, ha );
fhmax = max( fa, ha );
if ( fhmin == zero )
{
ssmin = zero;
if ( fhmax == zero )
ssmax = ga;
else
{
temp = min( fhmax, ga ) / max( fhmax, ga );
ssmax = max( fhmax, ga ) * sqrt( one + temp * temp );
}
}
else
{
if ( ga < fhmax )
{
as = one + fhmin / fhmax;
at = ( fhmax - fhmin ) / fhmax;
au = ( ga / fhmax ) * ( ga / fhmax );
c = two / ( sqrt( as * as + au ) + sqrt( at * at + au ) );
ssmin = fhmin * c;
ssmax = fhmax / c;
}
else
{
au = fhmax / ga;
if ( au == zero )
{
ssmin = ( fhmin * fhmax ) / ga;
ssmax = ga;
}
else
{
as = one + fhmin / fhmax;
at = ( fhmax - fhmin ) / fhmax;
temp = as * au;
temp2 = at * au;
c = one / ( sqrt( one + temp * temp ) +
sqrt( one + temp2 * temp2 ) );
ssmin = ( fhmin * c ) * au;
ssmin = ssmin + ssmin;
ssmax = ga / ( c + c );
}
}
}
// Save the output values.
*sigma1 = ssmin;
*sigma2 = ssmax;
return FLA_SUCCESS;
}
| FLA_Error FLA_Sv_2x2_ops | ( | float * | alpha11, |
| float * | alpha12, | ||
| float * | alpha22, | ||
| float * | sigma1, | ||
| float * | sigma2 | ||
| ) |
Referenced by FLA_Bsvd_compute_shift_ops(), and FLA_Sv_2x2().
{
float zero = 0.0F;
float one = 1.0F;
float two = 2.0F;
float f, g, h;
float as, at, au, c, fa, fhmin, fhmax, ga, ha;
float ssmin, ssmax;
float temp, temp2;
f = *alpha11;
g = *alpha12;
h = *alpha22;
fa = fabsf( f );
ga = fabsf( g );
ha = fabsf( h );
fhmin = min( fa, ha );
fhmax = max( fa, ha );
if ( fhmin == zero )
{
ssmin = zero;
if ( fhmax == zero )
ssmax = ga;
else
{
temp = min( fhmax, ga ) / max( fhmax, ga );
ssmax = max( fhmax, ga ) * sqrtf( one + temp * temp );
}
}
else
{
if ( ga < fhmax )
{
as = one + fhmin / fhmax;
at = ( fhmax - fhmin ) / fhmax;
au = ( ga / fhmax ) * ( ga / fhmax );
c = two / ( sqrtf( as * as + au ) + sqrtf( at * at + au ) );
ssmin = fhmin * c;
ssmax = fhmax / c;
}
else
{
au = fhmax / ga;
if ( au == zero )
{
ssmin = ( fhmin * fhmax ) / ga;
ssmax = ga;
}
else
{
as = one + fhmin / fhmax;
at = ( fhmax - fhmin ) / fhmax;
temp = as * au;
temp2 = at * au;
c = one / ( sqrtf( one + temp * temp ) +
sqrtf( one + temp2 * temp2 ) );
ssmin = ( fhmin * c ) * au;
ssmin = ssmin + ssmin;
ssmax = ga / ( c + c );
}
}
}
// Save the output values.
*sigma1 = ssmin;
*sigma2 = ssmax;
return FLA_SUCCESS;
}
| FLA_Error FLA_Svv_2x2 | ( | FLA_Obj | alpha11, |
| FLA_Obj | alpha12, | ||
| FLA_Obj | alpha22, | ||
| FLA_Obj | sigma1, | ||
| FLA_Obj | sigma2, | ||
| FLA_Obj | gammaL, | ||
| FLA_Obj | sigmaL, | ||
| FLA_Obj | gammaR, | ||
| FLA_Obj | sigmaR | ||
| ) |
References FLA_Obj_datatype(), FLA_Svv_2x2_opd(), and FLA_Svv_2x2_ops().
{
FLA_Datatype datatype;
datatype = FLA_Obj_datatype( alpha11 );
switch ( datatype )
{
case FLA_FLOAT:
{
float* buff_alpha11 = FLA_FLOAT_PTR( alpha11 );
float* buff_alpha12 = FLA_FLOAT_PTR( alpha12 );
float* buff_alpha22 = FLA_FLOAT_PTR( alpha22 );
float* buff_sigma1 = FLA_FLOAT_PTR( sigma1 );
float* buff_sigma2 = FLA_FLOAT_PTR( sigma2 );
float* buff_gammaL = FLA_FLOAT_PTR( gammaL );
float* buff_sigmaL = FLA_FLOAT_PTR( sigmaL );
float* buff_gammaR = FLA_FLOAT_PTR( gammaR );
float* buff_sigmaR = FLA_FLOAT_PTR( sigmaR );
FLA_Svv_2x2_ops( buff_alpha11,
buff_alpha12,
buff_alpha22,
buff_sigma1,
buff_sigma2,
buff_gammaL,
buff_sigmaL,
buff_gammaR,
buff_sigmaR );
break;
}
case FLA_DOUBLE:
{
double* buff_alpha11 = FLA_DOUBLE_PTR( alpha11 );
double* buff_alpha12 = FLA_DOUBLE_PTR( alpha12 );
double* buff_alpha22 = FLA_DOUBLE_PTR( alpha22 );
double* buff_sigma1 = FLA_DOUBLE_PTR( sigma1 );
double* buff_sigma2 = FLA_DOUBLE_PTR( sigma2 );
double* buff_gammaL = FLA_DOUBLE_PTR( gammaL );
double* buff_sigmaL = FLA_DOUBLE_PTR( sigmaL );
double* buff_gammaR = FLA_DOUBLE_PTR( gammaR );
double* buff_sigmaR = FLA_DOUBLE_PTR( sigmaR );
FLA_Svv_2x2_opd( buff_alpha11,
buff_alpha12,
buff_alpha22,
buff_sigma1,
buff_sigma2,
buff_gammaL,
buff_sigmaL,
buff_gammaR,
buff_sigmaR );
break;
}
}
return FLA_SUCCESS;
}
| FLA_Error FLA_Svv_2x2_opd | ( | double * | alpha11, |
| double * | alpha12, | ||
| double * | alpha22, | ||
| double * | sigma1, | ||
| double * | sigma2, | ||
| double * | gammaL, | ||
| double * | sigmaL, | ||
| double * | gammaR, | ||
| double * | sigmaR | ||
| ) |
References FLA_Mach_params_opd().
Referenced by FLA_Bsvd_iteracc_v_opd_var1(), and FLA_Svv_2x2().
{
double zero = 0.0;
double half = 0.5;
double one = 1.0;
double two = 2.0;
double four = 4.0;
double eps;
double f, g, h;
double clt, crt, slt, srt;
double a, d, fa, ft, ga, gt, ha, ht, l;
double m, mm, r, s, t, temp, tsign, tt;
double ssmin, ssmax;
double csl, snl;
double csr, snr;
int gasmal, swap;
int pmax;
f = *alpha11;
g = *alpha12;
h = *alpha22;
eps = FLA_Mach_params_opd( FLA_MACH_EPS );
ft = f;
fa = fabs( f );
ht = h;
ha = fabs( h );
// pmax points to the maximum absolute element of matrix.
// pmax = 1 if f largest in absolute values.
// pmax = 2 if g largest in absolute values.
// pmax = 3 if h largest in absolute values.
pmax = 1;
swap = ( ha > fa );
if ( swap )
{
pmax = 3;
temp = ft;
ft = ht;
ht = temp;
temp = fa;
fa = ha;
ha = temp;
}
gt = g;
ga = fabs( g );
if ( ga == zero )
{
// Diagonal matrix case.
ssmin = ha;
ssmax = fa;
clt = one;
slt = zero;
crt = one;
srt = zero;
}
else
{
gasmal = TRUE;
if ( ga > fa )
{
pmax = 2;
if ( ( fa / ga ) < eps )
{
// Case of very large ga.
gasmal = FALSE;
ssmax = ga;
if ( ha > one ) ssmin = fa / ( ga / ha );
else ssmin = ( fa / ga ) * ha;
clt = one;
slt = ht / gt;
crt = ft / gt;
srt = one;
}
}
if ( gasmal )
{
// Normal case.
d = fa - ha;
if ( d == fa ) l = one;
else l = d / fa;
m = gt / ft;
t = two - l;
mm = m * m;
tt = t * t;
s = sqrt( tt + mm );
if ( l == zero ) r = fabs( m );
else r = sqrt( l * l + mm );
a = half * ( s + r );
ssmin = ha / a;
ssmax = fa * a;
if ( mm == zero )
{
// Here, m is tiny.
if ( l == zero ) t = signof( two, ft ) * signof( one, gt );
else t = gt / signof( d, ft ) + m / t;
}
else
{
t = ( m / ( s + t ) + m / ( r + l ) ) * ( one + a );
}
l = sqrt( t*t + four );
crt = two / l;
srt = t / l;
clt = ( crt + srt * m ) / a;
slt = ( ht / ft ) * srt / a;
}
}
if ( swap )
{
csl = srt;
snl = crt;
csr = slt;
snr = clt;
}
else
{
csl = clt;
snl = slt;
csr = crt;
snr = srt;
}
// Correct the signs of ssmax and ssmin.
if ( pmax == 1 )
tsign = signof( one, csr ) * signof( one, csl ) * signof( one, f );
else if ( pmax == 2 )
tsign = signof( one, snr ) * signof( one, csl ) * signof( one, g );
else // if ( pmax == 3 )
tsign = signof( one, snr ) * signof( one, snl ) * signof( one, h );
ssmax = signof( ssmax, tsign );
ssmin = signof( ssmin, tsign * signof( one, f ) * signof( one, h ) );
// Save the output values.
*sigma1 = ssmin;
*sigma2 = ssmax;
*gammaL = csl;
*sigmaL = snl;
*gammaR = csr;
*sigmaR = snr;
return FLA_SUCCESS;
}
| FLA_Error FLA_Svv_2x2_ops | ( | float * | alpha11, |
| float * | alpha12, | ||
| float * | alpha22, | ||
| float * | sigma1, | ||
| float * | sigma2, | ||
| float * | gammaL, | ||
| float * | sigmaL, | ||
| float * | gammaR, | ||
| float * | sigmaR | ||
| ) |
References FLA_Mach_params_ops().
Referenced by FLA_Bsvd_iteracc_v_ops_var1(), and FLA_Svv_2x2().
{
float zero = 0.0F;
float half = 0.5F;
float one = 1.0F;
float two = 2.0F;
float four = 4.0F;
float eps;
float f, g, h;
float clt, crt, slt, srt;
float a, d, fa, ft, ga, gt, ha, ht, l;
float m, mm, r, s, t, temp, tsign, tt;
float ssmin, ssmax;
float csl, snl;
float csr, snr;
int gasmal, swap;
int pmax;
f = *alpha11;
g = *alpha12;
h = *alpha22;
eps = FLA_Mach_params_ops( FLA_MACH_EPS );
ft = f;
fa = fabsf( f );
ht = h;
ha = fabsf( h );
// pmax points to the maximum absolute element of matrix.
// pmax = 1 if f largest in absolute values.
// pmax = 2 if g largest in absolute values.
// pmax = 3 if h largest in absolute values.
pmax = 1;
swap = ( ha > fa );
if ( swap )
{
pmax = 3;
temp = ft;
ft = ht;
ht = temp;
temp = fa;
fa = ha;
ha = temp;
}
gt = g;
ga = fabsf( g );
if ( ga == zero )
{
// Diagonal matrix case.
ssmin = ha;
ssmax = fa;
clt = one;
slt = zero;
crt = one;
srt = zero;
}
else
{
gasmal = TRUE;
if ( ga > fa )
{
pmax = 2;
if ( ( fa / ga ) < eps )
{
// Case of very large ga.
gasmal = FALSE;
ssmax = ga;
if ( ha > one ) ssmin = fa / ( ga / ha );
else ssmin = ( fa / ga ) * ha;
clt = one;
slt = ht / gt;
crt = ft / gt;
srt = one;
}
}
if ( gasmal )
{
// Normal case.
d = fa - ha;
if ( d == fa ) l = one;
else l = d / fa;
m = gt / ft;
t = two - l;
mm = m * m;
tt = t * t;
s = sqrtf( tt + mm );
if ( l == zero ) r = fabsf( m );
else r = sqrtf( l * l + mm );
a = half * ( s + r );
ssmin = ha / a;
ssmax = fa * a;
if ( mm == zero )
{
// Here, m is tiny.
if ( l == zero ) t = signof( two, ft ) * signof( one, gt );
else t = gt / signof( d, ft ) + m / t;
}
else
{
t = ( m / ( s + t ) + m / ( r + l ) ) * ( one + a );
}
l = sqrtf( t*t + four );
crt = two / l;
srt = t / l;
clt = ( crt + srt * m ) / a;
slt = ( ht / ft ) * srt / a;
}
}
if ( swap )
{
csl = srt;
snl = crt;
csr = slt;
snr = clt;
}
else
{
csl = clt;
snl = slt;
csr = crt;
snr = srt;
}
// Correct the signs of ssmax and ssmin.
if ( pmax == 1 )
tsign = signof( one, csr ) * signof( one, csl ) * signof( one, f );
else if ( pmax == 2 )
tsign = signof( one, snr ) * signof( one, csl ) * signof( one, g );
else // if ( pmax == 3 )
tsign = signof( one, snr ) * signof( one, snl ) * signof( one, h );
ssmax = signof( ssmax, tsign );
ssmin = signof( ssmin, tsign * signof( one, f ) * signof( one, h ) );
// Save the output values.
*sigma1 = ssmin;
*sigma2 = ssmax;
*gammaL = csl;
*sigmaL = snl;
*gammaR = csr;
*sigmaR = snr;
return FLA_SUCCESS;
}
| FLA_Error FLA_Wilkshift_bidiag_check | ( | FLA_Obj | epsilon1, |
| FLA_Obj | delta1, | ||
| FLA_Obj | epsilon2, | ||
| 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_check | ( | FLA_Obj | delta1, |
| FLA_Obj | epsilon, | ||
| FLA_Obj | delta2, | ||
| FLA_Obj | kappa | ||
| ) |
References FLA_Check_identical_object_datatype(), FLA_Check_if_scalar(), FLA_Check_nonconstant_object(), and FLA_Check_real_object().
Referenced by FLA_Wilkshift_tridiag().
{
FLA_Error e_val;
e_val = FLA_Check_nonconstant_object( delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_real_object( delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, epsilon );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, delta2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_identical_object_datatype( delta1, kappa );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( delta1 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( epsilon );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( delta2 );
FLA_Check_error_code( e_val );
e_val = FLA_Check_if_scalar( kappa );
FLA_Check_error_code( e_val );
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;
}
1.7.6.1