libflame
revision_anchor
|
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) |
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; }