libflame  revision_anchor
Functions
FLA_util_lapack_prototypes.h File Reference

(r)

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_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)

Function Documentation

References bli_capdiagmv(), bli_csapdiagmv(), bli_dapdiagmv(), bli_sapdiagmv(), bli_zapdiagmv(), bli_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_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;
  side_t       blis_side; 
  conj_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 );

      bli_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 );

      bli_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 );

        bli_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 );

        bli_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 );

        bli_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 );

        bli_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;
}
FLA_Error FLA_Apply_G_1x2_check ( FLA_Obj  gamma,
FLA_Obj  sigma,
FLA_Obj  beta,
FLA_Obj  epsilon 
)

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;
}
FLA_Error FLA_Apply_G_mx2_check ( FLA_Obj  gamma,
FLA_Obj  sigma,
FLA_Obj  a1,
FLA_Obj  a2 
)

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_ */

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;
}
FLA_Error FLA_Househ2_UT ( FLA_Side  side,
FLA_Obj  chi_1,
FLA_Obj  x2,
FLA_Obj  tau 
)

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_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;
}
FLA_Error FLA_Househ2_UT_check ( FLA_Side  side,
FLA_Obj  chi_1,
FLA_Obj  x2,
FLA_Obj  tau 
)

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 bli_cinvscalv(), bli_cnrm2(), BLIS_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_cnrm2( i_two,
             y, i_one,
             &norm_x );

  //
  // Compute alpha:
  //
  //   alpha := - || x ||_2 * chi_1 / | chi_1 |
  //

  alpha.real = -chi_1->real * norm_x / abs_chi_1;
  alpha.imag = -chi_1->imag * norm_x / 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;

  bli_cinvscalv( BLIS_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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = chi_1_minus_alpha.real * chi_1_minus_alpha.real +
                             chi_1_minus_alpha.imag * chi_1_minus_alpha.imag;

  tau->real = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
              ( 2.0F * abs_sq_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 bli_dinvscalv(), bli_dnrm2(), BLIS_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_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;

  bli_dinvscalv( BLIS_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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = chi_1_minus_alpha * chi_1_minus_alpha;

  *tau = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
         ( 2.0 * abs_sq_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 bli_sinvscalv(), bli_snrm2(), BLIS_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_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;

  bli_sinvscalv( BLIS_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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = chi_1_minus_alpha * chi_1_minus_alpha;

  *tau = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
         ( 2.0F * abs_sq_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 bli_zinvscalv(), bli_znrm2(), BLIS_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_znrm2( i_two,
             y, i_one,
             &norm_x );

  //
  // Compute alpha:
  //
  //   alpha := - || x ||_2 * chi_1 / | chi_1 |
  //

  alpha.real = -chi_1->real * norm_x / abs_chi_1;
  alpha.imag = -chi_1->imag * norm_x / 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;

  bli_zinvscalv( BLIS_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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = chi_1_minus_alpha.real * chi_1_minus_alpha.real +
                             chi_1_minus_alpha.imag * chi_1_minus_alpha.imag;

  tau->real = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
              ( 2.0 * abs_sq_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 
)
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_ops ( int  m_x2,
float *  chi_1,
float *  x2,
int  inc_x2,
float *  tau 
)
FLA_Error FLA_Househ2_UT_r_opz ( int  m_x2,
dcomplex chi_1,
dcomplex x2,
int  inc_x2,
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 
)

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 bli_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_cnrm2( i_two,
             y, i_one,
             &norm_x );

  //
  // Compute alpha:
  //
  //   alpha := - || x ||_2 * chi_1 / | chi_1 |
  //

  alpha->real = -chi_1->real * norm_x / abs_chi_1;
  alpha->imag = -chi_1->imag * norm_x / 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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = chi_1_minus_alpha->real * chi_1_minus_alpha->real +
                             chi_1_minus_alpha->imag * chi_1_minus_alpha->imag;

  tau->real = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
              ( 2.0F * abs_sq_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 bli_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = (*chi_1_minus_alpha) * (*chi_1_minus_alpha);

  *tau = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
         ( 2.0 * abs_sq_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 bli_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = (*chi_1_minus_alpha) * (*chi_1_minus_alpha);

  *tau = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
         ( 2.0F * abs_sq_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 bli_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_sq_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
  //

  bli_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
  //

  bli_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;

  bli_znrm2( i_two,
             y, i_one,
             &norm_x );

  //
  // Compute alpha:
  //
  //   alpha := - || x ||_2 * chi_1 / | chi_1 |
  //

  alpha->real = -chi_1->real * norm_x / abs_chi_1;
  alpha->imag = -chi_1->imag * norm_x / 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 ) )
  //        = ( | chi_1 - alpha |^2 + || x_2 ||_2^2 ) /
  //          ( 2 * | chi_1 - alpha |^2 )
  //

  abs_sq_chi_1_minus_alpha = chi_1_minus_alpha->real * chi_1_minus_alpha->real +
                             chi_1_minus_alpha->imag * chi_1_minus_alpha->imag;

  tau->real = ( abs_sq_chi_1_minus_alpha + norm_x_2 * norm_x_2 ) /
              ( 2.0 * abs_sq_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;
}
FLA_Error FLA_Househ3UD_UT ( FLA_Obj  chi_1,
FLA_Obj  x2,
FLA_Obj  y2,
FLA_Obj  tau 
)

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;
}
FLA_Error FLA_Househ3UD_UT_check ( FLA_Obj  chi_1,
FLA_Obj  x2,
FLA_Obj  y2,
FLA_Obj  tau 
)

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 bli_cinvscalv(), bli_cnrm2(), BLIS_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
  //

  bli_cnrm2( m_x1,
             x1, inc_x1,
             &norm_x_1 );

  bli_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
  //

  bli_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;

  bli_cinvscalv( BLIS_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;

  bli_cinvscalv( BLIS_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 bli_dinvscalv(), bli_dnrm2(), BLIS_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
  //

  bli_dnrm2( m_x1,
             x1, inc_x1,
             &norm_x_1 );

  bli_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
  //

  bli_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;

  bli_dinvscalv( BLIS_NO_CONJUGATE,
                 m_x1,
                 &chi_0_minus_alpha,
                 x1, inc_x1 );

  neg_chi_0_minus_alpha = -chi_0_minus_alpha;

  bli_dinvscalv( BLIS_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 bli_sinvscalv(), bli_snrm2(), BLIS_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
  //

  bli_snrm2( m_x1,
             x1, inc_x1,
             &norm_x_1 );

  bli_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
  //

  bli_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;

  bli_sinvscalv( BLIS_NO_CONJUGATE,
                 m_x1,
                 &chi_0_minus_alpha,
                 x1, inc_x1 );

  neg_chi_0_minus_alpha = -chi_0_minus_alpha;

  bli_sinvscalv( BLIS_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 bli_zinvscalv(), bli_znrm2(), BLIS_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
  //

  bli_znrm2( m_x1,
             x1, inc_x1,
             &norm_x_1 );

  bli_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
  //

  bli_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;

  bli_zinvscalv( BLIS_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;

  bli_zinvscalv( BLIS_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;
}
logical fla_lsame ( char *  ca,
char *  cb,
ftnlen  ca_len,
ftnlen  cb_len 
)

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(), FLA_LU_piv(), 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;
}

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_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_Mach_params(), 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 bli_d0(), and bli_d1().

Referenced by FLA_Pythag2().

{
    double zero = bli_d0();
    double one  = bli_d1();

    double xabs, yabs;
    double w, z;
    double zdivw;

    xabs = fabs( *chi );
    yabs = fabs( *psi );
    w    = max( xabs, yabs );
    z    = min( xabs, yabs );

    if ( z == zero )
    {
        *rho = w;
    }
    else
    {
        zdivw = z / w;

        *rho = w * sqrt( one + zdivw * zdivw );
    }

    return FLA_SUCCESS;
}
FLA_Error FLA_Pythag2_ops ( float *  chi,
float *  psi,
float *  rho 
)

References bli_s0(), and bli_s1().

Referenced by FLA_Pythag2().

{
    float  zero = bli_s0();
    float  one  = bli_s1();

    float  xabs, yabs;
    float  w, z;
    float  zdivw;

    xabs = fabsf( *chi );
    yabs = fabsf( *psi );
    w    = max( xabs, yabs );
    z    = min( xabs, yabs );

    if ( z == zero )
    {
        *rho = w;
    }
    else
    {
        zdivw = z / w;

        *rho = w * sqrt( one + zdivw * zdivw );
    }

    return FLA_SUCCESS;
}
FLA_Error FLA_Pythag3 ( FLA_Obj  chi,
FLA_Obj  psi,
FLA_Obj  zeta,
FLA_Obj  rho 
)

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 bli_d0().

Referenced by FLA_Pythag3().

{
    double zero = bli_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 bli_s0().

Referenced by FLA_Pythag3().

{
    float  zero = bli_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;
}

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;
}

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_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 bli_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;
            bli_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 bli_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;
            bli_zswapv( m_A,
                        V + i*cs_V, rs_V,
                        V + k*cs_V, rs_V );
        }
    }

    return FLA_SUCCESS;
}

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 bli_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;
            bli_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 bli_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;
            bli_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 bli_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;
            bli_dswapv( m_U,
                        U + i*cs_U, rs_U,
                        U + k*cs_U, rs_U );
            bli_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 bli_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;
            bli_zswapv( m_U,
                        U + i*cs_U, rs_U,
                        U + k*cs_U, rs_U );
            bli_zswapv( n_V,
                        V + i*cs_V, rs_V,
                        V + k*cs_V, rs_V );
        }
    }

    return FLA_SUCCESS;
}

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_square(), 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 bli_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;
            bli_dswapv( m_U,
                        U + i*cs_U, rs_U,
                        U + k*cs_U, rs_U );
            bli_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 bli_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;
            bli_zswapv( m_U,
                        U + i*cs_U, rs_U,
                        U + k*cs_U, rs_U );
            bli_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_Sv_2x2().

{
    FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );

    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 
)

Referenced by FLA_Svv_2x2().

{
    FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );

    return FLA_SUCCESS;
}
FLA_Error FLA_Wilkshift_bidiag_check ( FLA_Obj  epsilon1,
FLA_Obj  delta1,
FLA_Obj  epsilon2,
FLA_Obj  delta2,
FLA_Obj  kappa 
)
FLA_Error FLA_Wilkshift_tridiag ( FLA_Obj  delta1,
FLA_Obj  epsilon,
FLA_Obj  delta2,
FLA_Obj  kappa 
)

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

{
    FLA_Datatype datatype;

    datatype = FLA_Obj_datatype( delta1 );

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

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

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

            break;
        }

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

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

            break;
        }
    }

    return FLA_SUCCESS;
}
FLA_Error FLA_Wilkshift_tridiag_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;
}