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

FLA_Error FLA_Apply_diag_matrix ( FLA_Side  side,
FLA_Conj  conj,
FLA_Obj  x,
FLA_Obj  A 
)

References bli_capdiagmv(), bli_dapdiagmv(), bli_sapdiagmv(), bli_zapdiagmv(), 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().

{
  FLA_Datatype datatype;
  int          m_A, n_A;
  int          rs_A, cs_A;
  int          inc_x;
  char         blis_side; 
  char         blis_conj;

  if ( FLA_Check_error_level() == FLA_MIN_ERROR_CHECKING ) 
    FLA_Apply_diag_matrix_check( side, conj, x, A );

  datatype = 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 ( datatype )
  {
    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:
    {
      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:
    {
      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;
}
FLA_Error FLA_Apply_diag_matrix_check ( FLA_Side  side,
FLA_Conj  conj,
FLA_Obj  x,
FLA_Obj  A 
)

References FLA_Check_floating_object(), FLA_Check_identical_object_datatype(), FLA_Check_nonconstant_object(), FLA_Check_valid_conj(), and FLA_Check_valid_leftright_side().

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_datatype( A, x );
  FLA_Check_error_code( e_val );

  return FLA_SUCCESS;
}
FLA_Error FLA_Form_perm_matrix ( FLA_Obj  p,
FLA_Obj  A 
)

References FLA_Apply_pivots(), FLA_Check_error_level(), FLA_Form_perm_matrix_check(), and FLA_Obj_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_Obj_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;
}
FLA_Error FLA_Form_perm_matrix_check ( FLA_Obj  p,
FLA_Obj  A 
)

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_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_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(), 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_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(), 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_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(), 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_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(), 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_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(), 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(), 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(), 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(), 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_LU_find_zero_on_diagonal ( FLA_Obj  A)

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;
}
FLA_Error FLA_LU_find_zero_on_diagonal_check ( FLA_Obj  A)

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_Shift_pivots_to ( FLA_Pivot_type  ptype,
FLA_Obj  p 
)

References FLA_Check_error_level(), FLA_Obj_length(), FLA_Obj_width(), and FLA_Shift_pivots_to_check().

Referenced by FLA_LU_piv_blk_external(), and FLA_LU_piv_unb_external().

{
  int  m_p, n_p;
  int* buff_p;
  int  i;

  if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
    FLA_Shift_pivots_to_check( ptype, p );

  m_p    = FLA_Obj_length( p );
  n_p    = FLA_Obj_width( p );
  buff_p = FLA_INT_PTR( p );

  if ( m_p < 1 || n_p < 1 ) return FLA_SUCCESS;

  if ( ptype == FLA_LAPACK_PIVOTS )
  {
    // Shift FLAME pivots to LAPACK pivots.
    for ( i = 0; i < m_p; i++ )
      buff_p[ i ] += i + 1;
  }
  else
  {
    // Otherwise, shift LAPACK pivots back to FLAME.
    for ( i = 0; i < m_p; i++ )
      buff_p[ i ] -= i + 1;
  }

  return FLA_SUCCESS;
}
FLA_Error FLA_Shift_pivots_to_check ( FLA_Pivot_type  ptype,
FLA_Obj  p 
)

References FLA_Check_col_vector(), FLA_Check_int_object(), FLA_Check_nonconstant_object(), and FLA_Check_valid_pivot_type().

Referenced by FLA_Shift_pivots_to().

{
  FLA_Error e_val;

  e_val = FLA_Check_valid_pivot_type( ptype );
  FLA_Check_error_code( e_val );

  e_val = FLA_Check_int_object( p );
  FLA_Check_error_code( e_val );

  e_val = FLA_Check_nonconstant_object( p );
  FLA_Check_error_code( e_val );

  e_val = FLA_Check_col_vector( p );
  FLA_Check_error_code( e_val );

  return FLA_SUCCESS;
}