libflame  revision_anchor
Functions
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_Svd()

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

14 {
15  FLA_Error r_val = FLA_SUCCESS;
16  dim_t n_iter_max = 30;
17  dim_t k_accum = 32;
18  dim_t b_alg = 512;
19  dim_t min_m_n = FLA_Obj_min_dim( A );
20  dim_t m_A = FLA_Obj_length( A );
21  dim_t n_A = FLA_Obj_width( A );
22  FLA_Obj W; // Dummy variable for partitioning of matrices.
23 
24  // Check parameters.
25  if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
26  FLA_Svd_check( jobu, jobv, A, s, U, V );
27 
28  // Partition U and V if necessary.
29  if ( jobu == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( U, &U, &W, min_m_n, FLA_LEFT );
30  if ( jobv == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( V, &V, &W, min_m_n, FLA_LEFT );
31 
32  // Use extension version
33  if ( m_A >= n_A )
34  {
35  r_val = FLA_Svd_ext_u_unb_var1( jobu, jobv,
36  n_iter_max,
37  A, s, U, V,
38  k_accum, b_alg );
39  }
40  else
41  {
42  // Flip A and change U and V
43  FLA_Obj_flip_base( &A );
44  FLA_Obj_flip_view( &A );
45 
46  r_val = FLA_Svd_ext_u_unb_var1( jobu, jobv,
47  n_iter_max,
48  A, s, V, U,
49  k_accum, b_alg );
50 
51  // Recover A and conjugate U and V for complex cases
52  FLA_Obj_flip_base( &A );
53 
54  if ( FLA_Obj_is_complex( A ) )
55  {
56  if ( jobu != FLA_SVD_VECTORS_NONE ) FLA_Conjugate( U );
57  if ( jobv != FLA_SVD_VECTORS_NONE ) FLA_Conjugate( V );
58  }
59  }
60 
61  return r_val;
62 }
unsigned long dim_t
Definition: FLA_type_defs.h:71
FLA_Error FLA_Obj_flip_base(FLA_Obj *obj)
Definition: FLA_Obj.c:647
FLA_Error FLA_Conjugate(FLA_Obj A)
Definition: FLA_Conjugate.c:13
int FLA_Error
Definition: FLA_type_defs.h:47
Definition: FLA_type_defs.h:158
dim_t FLA_Obj_width(FLA_Obj obj)
Definition: FLA_Query.c:123
FLA_Error FLA_Svd_ext_u_unb_var1(FLA_Svd_type jobu, FLA_Svd_type jobv, dim_t n_iter_max, FLA_Obj A, FLA_Obj s, FLA_Obj V, FLA_Obj U, dim_t k_accum, dim_t b_alg)
Definition: FLA_Svd_ext_u_unb_var1.c:14
FLA_Error FLA_Svd_check(FLA_Svd_type jobu, FLA_Svd_type jobv, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V)
Definition: FLA_Svd_check.c:13
unsigned int FLA_Check_error_level(void)
Definition: FLA_Check.c:18
FLA_Error FLA_Part_1x2(FLA_Obj A, FLA_Obj *A1, FLA_Obj *A2, dim_t nb, FLA_Side side)
Definition: FLA_View.c:110
FLA_Bool FLA_Obj_is_complex(FLA_Obj A)
Definition: FLA_Query.c:324
dim_t FLA_Obj_length(FLA_Obj obj)
Definition: FLA_Query.c:116
FLA_Error FLA_Obj_flip_view(FLA_Obj *obj)
Definition: FLA_Obj.c:669
dim_t FLA_Obj_min_dim(FLA_Obj obj)
Definition: FLA_Query.c:153

◆ FLA_Svd_compute_scaling()

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

14 {
15  FLA_Datatype dt_real;
16  FLA_Obj norm;
17  FLA_Obj safmin;
18  FLA_Obj prec;
19  FLA_Obj rmin;
20  FLA_Obj rmax;
21 
22  if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
24 
25  dt_real = FLA_Obj_datatype_proj_to_real( A );
26 
27  FLA_Obj_create( dt_real, 1, 1, 0, 0, &norm );
28  FLA_Obj_create( dt_real, 1, 1, 0, 0, &prec );
29  FLA_Obj_create( dt_real, 1, 1, 0, 0, &safmin );
30  FLA_Obj_create( dt_real, 1, 1, 0, 0, &rmin );
31  FLA_Obj_create( dt_real, 1, 1, 0, 0, &rmax );
32 
33  // Query safmin, precision.
34  FLA_Mach_params( FLA_MACH_PREC, prec );
35  FLA_Mach_params( FLA_MACH_SFMIN, safmin );
36 
37 //FLA_Obj_show( "safmin", safmin, "%20.12e", "" );
38 //FLA_Obj_show( "prec", prec, "%20.12e", "" );
39 
40  // rmin = sqrt( safmin ) / prec;
41  FLA_Copy( safmin, rmin );
42  FLA_Sqrt( rmin );
43  FLA_Inv_scal( prec, rmin );
44 
45  // rmax = 1 / rmin;
46  FLA_Copy( rmin, rmax );
47  FLA_Invert( FLA_NO_CONJUGATE, rmax );
48 
49 //FLA_Obj_show( "rmin", rmin, "%20.12e", "" );
50 //FLA_Obj_show( "rmax", rmax, "%20.12e", "" );
51 
52  // Find the maximum absolute value of A.
53  FLA_Max_abs_value( A, norm );
54 
55  if ( FLA_Obj_gt( norm, FLA_ZERO ) && FLA_Obj_lt( norm, rmin ) )
56  {
57  // sigma = rmin / norm;
58  FLA_Copy( rmin, sigma );
59  FLA_Inv_scal( norm, sigma );
60  }
61  else if ( FLA_Obj_gt( norm, rmax ) )
62  {
63  // sigma = rmax / norm;
64  FLA_Copy( rmax, sigma );
65  FLA_Inv_scal( norm, sigma );
66  }
67  else
68  {
69  // sigma = 1.0;
70  FLA_Copy( FLA_ONE, sigma );
71  }
72 
73  FLA_Obj_free( &norm );
74  FLA_Obj_free( &prec );
75  FLA_Obj_free( &safmin );
76  FLA_Obj_free( &rmin );
77  FLA_Obj_free( &rmax );
78 
79  return FLA_SUCCESS;
80 }
FLA_Error FLA_Obj_create(FLA_Datatype datatype, dim_t m, dim_t n, dim_t rs, dim_t cs, FLA_Obj *obj)
Definition: FLA_Obj.c:55
FLA_Error FLA_Obj_free(FLA_Obj *obj)
Definition: FLA_Obj.c:588
FLA_Error FLA_Copy(FLA_Obj A, FLA_Obj B)
Definition: FLA_Copy.c:15
FLA_Bool FLA_Obj_gt(FLA_Obj A, FLA_Obj B)
Definition: FLA_Query.c:658
FLA_Error FLA_Max_abs_value(FLA_Obj A, FLA_Obj amax)
Definition: FLA_Max_abs_value.c:13
FLA_Obj FLA_ONE
Definition: FLA_Init.c:18
FLA_Bool FLA_Obj_lt(FLA_Obj A, FLA_Obj B)
Definition: FLA_Query.c:813
FLA_Error FLA_Mach_params(FLA_Machval machval, FLA_Obj val)
Definition: FLA_Mach_params.c:13
Definition: FLA_type_defs.h:158
FLA_Error FLA_Svd_compute_scaling_check(FLA_Obj A, FLA_Obj sigma)
Definition: FLA_Svd_compute_scaling_check.c:13
FLA_Error FLA_Invert(FLA_Conj conj, FLA_Obj x)
Definition: FLA_Invert.c:13
FLA_Datatype FLA_Obj_datatype_proj_to_real(FLA_Obj A)
Definition: FLA_Query.c:23
FLA_Error FLA_Inv_scal(FLA_Obj alpha, FLA_Obj A)
Definition: FLA_Inv_scal.c:13
unsigned int FLA_Check_error_level(void)
Definition: FLA_Check.c:18
int FLA_Datatype
Definition: FLA_type_defs.h:49
FLA_Error FLA_Sqrt(FLA_Obj alpha)
Definition: FLA_Sqrt.c:13
FLA_Obj FLA_ZERO
Definition: FLA_Init.c:20

◆ FLA_Svd_ext()

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

28 {
29  FLA_Error r_val = FLA_SUCCESS;
30  dim_t n_iter_max = 30;
31  dim_t k_accum = 32;
32  dim_t b_alg = 512;
33  dim_t min_m_n = FLA_Obj_min_dim( A );
34  dim_t m_A = FLA_Obj_length( A );
35  dim_t n_A = FLA_Obj_width( A );
36  FLA_Bool u_flipped = FALSE;
37  FLA_Bool v_flipped = FALSE;
38  FLA_Obj W; // Dummy variable for partitioning of matrices.
39 
40  // Check parameters.
41  if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
42  FLA_Svd_ext_check( jobu, transu, jobv, transv, A, s, U, V );
43 
44  // Transpose U and V to match dimensions used in SVD.
45  if ( ( transu == FLA_TRANSPOSE || transu == FLA_CONJ_TRANSPOSE ) &&
46  ( jobu != FLA_SVD_VECTORS_NONE && jobu != FLA_SVD_VECTORS_MIN_OVERWRITE ) )
47  {
48  FLA_Obj_flip_base( &U );
49  FLA_Obj_flip_view( &U );
50  u_flipped = TRUE;
51  }
52  if ( ( transv == FLA_TRANSPOSE || transv == FLA_CONJ_TRANSPOSE ) &&
53  ( jobv != FLA_SVD_VECTORS_NONE && jobv != FLA_SVD_VECTORS_MIN_OVERWRITE ) )
54  {
55  FLA_Obj_flip_base( &V );
56  FLA_Obj_flip_view( &V );
57  v_flipped = TRUE;
58  }
59 
60  // Partition U and V if necessary.
61  if ( jobu == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( U, &U, &W, min_m_n, FLA_LEFT );
62  if ( jobv == FLA_SVD_VECTORS_MIN_COPY ) FLA_Part_1x2( V, &V, &W, min_m_n, FLA_LEFT );
63 
64  if ( m_A >= n_A )
65  {
66  r_val = FLA_Svd_ext_u_unb_var1( jobu, jobv,
67  n_iter_max,
68  A, s, U, V,
69  k_accum, b_alg );
70 
71  // Recover U and V.
72  if ( u_flipped == TRUE )
73  {
74  if ( FLA_Obj_is_complex( U ) )
75  FLA_Conjugate( U );
76  FLA_Obj_flip_base( &U );
77  }
78  if ( v_flipped == TRUE )
79  {
80  if ( FLA_Obj_is_complex( V ) )
81  FLA_Conjugate( V );
82  FLA_Obj_flip_base( &V );
83  }
84  }
85  else
86  {
87  // Flip A and exchange U and V parameters.
88  FLA_Obj_flip_base( &A );
89  FLA_Obj_flip_view( &A );
90 
91  // Note that U and V are also swapped.
92  r_val = FLA_Svd_ext_u_unb_var1( jobv, jobu,
93  n_iter_max,
94  A, s, V, U,
95  k_accum, b_alg );
96 
97  // Recover A.
98  FLA_Obj_flip_base( &A );
99 
100  // Recover U and V. Consider a case that U and V are not created.
101  if ( u_flipped == TRUE )
102  FLA_Obj_flip_base( &U );
103  else if ( jobu != FLA_SVD_VECTORS_NONE &&
104  jobu != FLA_SVD_VECTORS_MIN_OVERWRITE )
105  if ( FLA_Obj_is_complex( U ) )
106  FLA_Conjugate( U );
107 
108  if ( v_flipped == TRUE )
109  FLA_Obj_flip_base( &V );
110  else if ( jobv != FLA_SVD_VECTORS_NONE &&
111  jobv != FLA_SVD_VECTORS_MIN_OVERWRITE )
112  if ( FLA_Obj_is_complex( V ) )
113  FLA_Conjugate( V );
114  }
115 
116  return r_val;
117 }
unsigned long dim_t
Definition: FLA_type_defs.h:71
FLA_Error FLA_Obj_flip_base(FLA_Obj *obj)
Definition: FLA_Obj.c:647
FLA_Error FLA_Conjugate(FLA_Obj A)
Definition: FLA_Conjugate.c:13
int FLA_Error
Definition: FLA_type_defs.h:47
Definition: FLA_type_defs.h:158
dim_t FLA_Obj_width(FLA_Obj obj)
Definition: FLA_Query.c:123
FLA_Error FLA_Svd_ext_u_unb_var1(FLA_Svd_type jobu, FLA_Svd_type jobv, dim_t n_iter_max, FLA_Obj A, FLA_Obj s, FLA_Obj V, FLA_Obj U, dim_t k_accum, dim_t b_alg)
Definition: FLA_Svd_ext_u_unb_var1.c:14
FLA_Error FLA_Svd_ext_check(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)
Definition: FLA_Svd_ext_check.c:13
unsigned int FLA_Check_error_level(void)
Definition: FLA_Check.c:18
FLA_Error FLA_Part_1x2(FLA_Obj A, FLA_Obj *A1, FLA_Obj *A2, dim_t nb, FLA_Side side)
Definition: FLA_View.c:110
int FLA_Bool
Definition: FLA_type_defs.h:46
FLA_Bool FLA_Obj_is_complex(FLA_Obj A)
Definition: FLA_Query.c:324
dim_t FLA_Obj_length(FLA_Obj obj)
Definition: FLA_Query.c:116
FLA_Error FLA_Obj_flip_view(FLA_Obj *obj)
Definition: FLA_Obj.c:669
dim_t FLA_Obj_min_dim(FLA_Obj obj)
Definition: FLA_Query.c:153