FLA_Svd.h File Reference

(r)

Go to the source code of this file.

Functions
FLA_Error	FLA_Svd_compute_scaling (FLA_Obj A, FLA_Obj sigma)
FLA_Error	FLA_Svd (FLA_Svd_type jobu, FLA_Svd_type jobv, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V)
FLA_Error	FLA_Svd_ext (FLA_Svd_type jobu, FLA_Trans transu, FLA_Svd_type jobv, FLA_Trans transv, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V)

---

Function Documentation

FLA_Error FLA_Svd	(	FLA_Svd_type	jobu,
		FLA_Svd_type	jobv,
		FLA_Obj	A,
		FLA_Obj	s,
		FLA_Obj	U,
		FLA_Obj	V
	)

References FLA_Check_error_level(), FLA_Conjugate(), FLA_Obj_flip_base(), FLA_Obj_flip_view(), FLA_Obj_is_complex(), FLA_Obj_length(), FLA_Obj_min_dim(), FLA_Obj_width(), FLA_Part_1x2(), FLA_Svd_check(), and FLA_Svd_ext_u_unb_var1().

{
  FLA_Error r_val      = FLA_SUCCESS;
  dim_t     n_iter_max = 30;
  dim_t     k_accum    = 32;
  dim_t     b_alg      = 512;
  dim_t     min_m_n    = FLA_Obj_min_dim( A );
  dim_t     m_A        = FLA_Obj_length( A );
  dim_t     n_A        = FLA_Obj_width( A );
  FLA_Obj   W;         // Dummy variable for partitioning of matrices.

  // Check parameters.
  if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
    FLA_Svd_check( jobu, jobv, A, s, U, V );

  // Partition U and V if necessary.
  if ( jobu == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( U, &U, &W, min_m_n, FLA_LEFT );
  if ( jobv == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( V, &V, &W, min_m_n, FLA_LEFT );

  // Use extension version
  if ( m_A >= n_A )
  {
    r_val = FLA_Svd_ext_u_unb_var1( jobu, jobv, 
                                    n_iter_max, 
                                    A, s, U, V,
                                    k_accum, b_alg );
  }
  else
  {
    // Flip A and change U and V
    FLA_Obj_flip_base( &A );
    FLA_Obj_flip_view( &A );
    
    r_val = FLA_Svd_ext_u_unb_var1( jobu, jobv, 
                                    n_iter_max, 
                                    A, s, V, U,
                                    k_accum, b_alg );
    
    // Recover A and conjugate U and V for complex cases
    FLA_Obj_flip_base( &A );
    
    if ( FLA_Obj_is_complex( A ) )
    {
      if ( jobu != FLA_SVD_VECTORS_NONE ) FLA_Conjugate( U );
      if ( jobv != FLA_SVD_VECTORS_NONE ) FLA_Conjugate( V );
    }
  }

  return r_val;
}

FLA_Error FLA_Svd_compute_scaling	(	FLA_Obj	A,
		FLA_Obj	sigma
	)

References FLA_Check_error_level(), FLA_Copy(), FLA_Inv_scal(), FLA_Invert(), FLA_Mach_params(), FLA_Max_abs_value(), FLA_Obj_create(), FLA_Obj_datatype_proj_to_real(), FLA_Obj_free(), FLA_Obj_gt(), FLA_Obj_lt(), FLA_ONE, FLA_Sqrt(), FLA_Svd_compute_scaling_check(), and FLA_ZERO.

Referenced by FLA_Svd_uv_unb_var1(), and FLA_Svd_uv_unb_var2().

{
    FLA_Datatype dt_real;
    FLA_Obj      norm;
    FLA_Obj      safmin;
    FLA_Obj      prec;
    FLA_Obj      rmin;
    FLA_Obj      rmax;

    if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
        FLA_Svd_compute_scaling_check( A, sigma );

    dt_real = FLA_Obj_datatype_proj_to_real( A );

    FLA_Obj_create( dt_real, 1, 1, 0, 0, &norm );
    FLA_Obj_create( dt_real, 1, 1, 0, 0, &prec );
    FLA_Obj_create( dt_real, 1, 1, 0, 0, &safmin );
    FLA_Obj_create( dt_real, 1, 1, 0, 0, &rmin );
    FLA_Obj_create( dt_real, 1, 1, 0, 0, &rmax );

    // Query safmin, precision.
    FLA_Mach_params( FLA_MACH_PREC,  prec );
    FLA_Mach_params( FLA_MACH_SFMIN, safmin );

//FLA_Obj_show( "safmin", safmin, "%20.12e", "" );
//FLA_Obj_show( "prec", prec, "%20.12e", "" );

    // rmin = sqrt( safmin ) / prec;
    FLA_Copy( safmin, rmin );
    FLA_Sqrt( rmin );
    FLA_Inv_scal( prec, rmin );

    // rmax = 1 / rmin;
    FLA_Copy( rmin, rmax );
    FLA_Invert( FLA_NO_CONJUGATE, rmax );

//FLA_Obj_show( "rmin", rmin, "%20.12e", "" );
//FLA_Obj_show( "rmax", rmax, "%20.12e", "" );

    // Find the maximum absolute value of A.
    FLA_Max_abs_value( A, norm );

    if ( FLA_Obj_gt( norm, FLA_ZERO ) && FLA_Obj_lt( norm, rmin ) )
    {
        // sigma = rmin / norm;
        FLA_Copy( rmin, sigma );
        FLA_Inv_scal( norm, sigma );
    }
    else if ( FLA_Obj_gt( norm, rmax ) )
    {
        // sigma = rmax / norm;
        FLA_Copy( rmax, sigma );
        FLA_Inv_scal( norm, sigma );
    }
    else
    {
        // sigma = 1.0;
        FLA_Copy( FLA_ONE, sigma );
    }

    FLA_Obj_free( &norm );
    FLA_Obj_free( &prec );
    FLA_Obj_free( &safmin );
    FLA_Obj_free( &rmin );
    FLA_Obj_free( &rmax );

    return FLA_SUCCESS;
}

FLA_Error FLA_Svd_ext	(	FLA_Svd_type	jobu,
		FLA_Trans	transu,
		FLA_Svd_type	jobv,
		FLA_Trans	transv,
		FLA_Obj	A,
		FLA_Obj	s,
		FLA_Obj	U,
		FLA_Obj	V
	)

References FLA_Check_error_level(), FLA_Conjugate(), FLA_Obj_flip_base(), FLA_Obj_flip_view(), FLA_Obj_is_complex(), FLA_Obj_length(), FLA_Obj_min_dim(), FLA_Obj_width(), FLA_Part_1x2(), FLA_Svd_ext_check(), and FLA_Svd_ext_u_unb_var1().

{
  FLA_Error r_val      = FLA_SUCCESS;
  dim_t     n_iter_max = 30;
  dim_t     k_accum    = 32;
  dim_t     b_alg      = 512;
  dim_t     min_m_n    = FLA_Obj_min_dim( A );
  dim_t     m_A        = FLA_Obj_length( A );
  dim_t     n_A        = FLA_Obj_width( A );
  FLA_Bool  u_flipped  = FALSE;
  FLA_Bool  v_flipped  = FALSE;
  FLA_Obj   W;           // Dummy variable for partitioning of matrices.

  // Check parameters.
  if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
    FLA_Svd_ext_check( jobu, transu, jobv, transv, A, s, U, V );

  // Transpose U and V to match dimensions used in SVD. 
  if ( ( transu == FLA_TRANSPOSE        || transu == FLA_CONJ_TRANSPOSE            ) &&
       ( jobu   != FLA_SVD_VECTORS_NONE && jobu   != FLA_SVD_VECTORS_MIN_OVERWRITE ) )
  {
    FLA_Obj_flip_base( &U );
    FLA_Obj_flip_view( &U );
    u_flipped = TRUE;
  }
  if ( ( transv == FLA_TRANSPOSE        || transv == FLA_CONJ_TRANSPOSE            ) &&
       ( jobv   != FLA_SVD_VECTORS_NONE && jobv   != FLA_SVD_VECTORS_MIN_OVERWRITE ) )
  {
    FLA_Obj_flip_base( &V );
    FLA_Obj_flip_view( &V );
    v_flipped = TRUE;
  }

  // Partition U and V if necessary.
  if ( jobu == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( U, &U, &W, min_m_n, FLA_LEFT );
  if ( jobv == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( V, &V, &W, min_m_n, FLA_LEFT );

  if ( m_A >= n_A )
  {
    r_val = FLA_Svd_ext_u_unb_var1( jobu, jobv, 
                                    n_iter_max, 
                                    A, s, U, V, 
                                    k_accum, b_alg );

    // Recover U and V.
    if ( u_flipped == TRUE ) 
    {
      if ( FLA_Obj_is_complex( U ) )
        FLA_Conjugate( U );
      FLA_Obj_flip_base( &U );
    }
    if ( v_flipped == TRUE ) 
    {
      if ( FLA_Obj_is_complex( V ) )
        FLA_Conjugate( V );
      FLA_Obj_flip_base( &V );
    }
  }
  else
  {
    // Flip A and exchange U and V parameters.
    FLA_Obj_flip_base( &A );
    FLA_Obj_flip_view( &A );

    // Note that U and V are also swapped.
    r_val = FLA_Svd_ext_u_unb_var1( jobv, jobu, 
                                    n_iter_max, 
                                    A, s, V, U, 
                                    k_accum, b_alg );
    
    // Recover A.
    FLA_Obj_flip_base( &A );

    // Recover U and V. Consider a case that U and V are not created. 
    if ( u_flipped == TRUE ) 
      FLA_Obj_flip_base( &U );
    else if ( jobu != FLA_SVD_VECTORS_NONE && 
              jobu != FLA_SVD_VECTORS_MIN_OVERWRITE ) 
      if ( FLA_Obj_is_complex( U ) )
        FLA_Conjugate( U );
    
    if ( v_flipped == TRUE ) 
      FLA_Obj_flip_base( &V );
    else if ( jobv != FLA_SVD_VECTORS_NONE && 
              jobv != FLA_SVD_VECTORS_MIN_OVERWRITE )
      if ( FLA_Obj_is_complex( V ) )
        FLA_Conjugate( V );
  }

  return r_val;
}