libflame  revision_anchor
Functions
FLA_Svd_uv.h File Reference

(r)

Go to the source code of this file.

Functions

FLA_Error FLA_Svd_uv_unb_var1 (dim_t n_iter_max, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V, dim_t k_accum, dim_t b_alg)
 
FLA_Error FLA_Svd_uv_unb_var2 (dim_t n_iter_max, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V, dim_t k_accum, dim_t b_alg)
 

Function Documentation

◆ FLA_Svd_uv_unb_var1()

FLA_Error FLA_Svd_uv_unb_var1 ( dim_t  n_iter_max,
FLA_Obj  A,
FLA_Obj  s,
FLA_Obj  U,
FLA_Obj  V,
dim_t  k_accum,
dim_t  b_alg 
)
14 {
15  FLA_Error r_val = FLA_SUCCESS;
16  FLA_Datatype dt;
17  FLA_Datatype dt_real;
18  FLA_Datatype dt_comp;
19  FLA_Obj scale, T, S, rL, rR, d, e, G, H;
20  dim_t m_A, n_A;
21  dim_t min_m_n;
22  dim_t n_GH;
23  double crossover_ratio = 17.0 / 9.0;
24 
25  n_GH = k_accum;
26 
27  m_A = FLA_Obj_length( A );
28  n_A = FLA_Obj_width( A );
29  min_m_n = FLA_Obj_min_dim( A );
30  dt = FLA_Obj_datatype( A );
31  dt_real = FLA_Obj_datatype_proj_to_real( A );
32  dt_comp = FLA_Obj_datatype_proj_to_complex( A );
33 
34  // Create matrices to hold block Householder transformations.
35  FLA_Bidiag_UT_create_T( A, &T, &S );
36 
37  // Create vectors to hold the realifying scalars.
38  FLA_Obj_create( dt, min_m_n, 1, 0, 0, &rL );
39  FLA_Obj_create( dt, min_m_n, 1, 0, 0, &rR );
40 
41  // Create vectors to hold the diagonal and sub-diagonal.
42  FLA_Obj_create( dt_real, min_m_n, 1, 0, 0, &d );
43  FLA_Obj_create( dt_real, min_m_n-1, 1, 0, 0, &e );
44 
45  // Create matrices to hold the left and right Givens scalars.
46  FLA_Obj_create( dt_comp, min_m_n-1, n_GH, 0, 0, &G );
47  FLA_Obj_create( dt_comp, min_m_n-1, n_GH, 0, 0, &H );
48 
49  // Create a real scaling factor.
50  FLA_Obj_create( dt_real, 1, 1, 0, 0, &scale );
51 
52  // Compute a scaling factor; If none is needed, sigma will be set to one.
53  FLA_Svd_compute_scaling( A, scale );
54 
55  // Scale the matrix if scale is non-unit.
56  if ( !FLA_Obj_equals( scale, FLA_ONE ) )
57  FLA_Scal( scale, A );
58 
59  if ( m_A < crossover_ratio * n_A )
60  {
61  // Reduce the matrix to bidiagonal form.
62  // Apply scalars to rotate elements on the superdiagonal to the real domain.
63  // Extract the diagonal and superdiagonal from A.
64  FLA_Bidiag_UT( A, T, S );
65  FLA_Bidiag_UT_realify( A, rL, rR );
67 
68  // Form U and V.
69  FLA_Bidiag_UT_form_U( A, T, U );
70  FLA_Bidiag_UT_form_V( A, S, V );
71 
72  // Apply the realifying scalars in rL and rR to U and V, respectively.
73  {
74  FLA_Obj UL, UR;
75  FLA_Obj VL, VR;
76 
77  FLA_Part_1x2( U, &UL, &UR, min_m_n, FLA_LEFT );
78  FLA_Part_1x2( V, &VL, &VR, min_m_n, FLA_LEFT );
79 
80  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_CONJUGATE, rL, UL );
81  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_NO_CONJUGATE, rR, VL );
82  }
83 
84  // Perform a singular value decomposition on the bidiagonal matrix.
85  r_val = FLA_Bsvd_v_opt_var1( n_iter_max, d, e, G, H, U, V, b_alg );
86  }
87  else // if ( crossover_ratio * n_A <= m_A )
88  {
89  FLA_Obj TQ, R;
90  FLA_Obj AT,
91  AB;
92  FLA_Obj UL, UR;
93 
94  // Perform a QR factorization on A and form Q in U.
95  FLA_QR_UT_create_T( A, &TQ );
96  FLA_QR_UT( A, TQ );
97  FLA_QR_UT_form_Q( A, TQ, U );
98  FLA_Obj_free( &TQ );
99 
100  // Set the lower triangle of R to zero and then copy the upper
101  // triangle of A to R.
102  FLA_Part_2x1( A, &AT,
103  &AB, n_A, FLA_TOP );
104  FLA_Obj_create( dt, n_A, n_A, 0, 0, &R );
105  FLA_Setr( FLA_LOWER_TRIANGULAR, FLA_ZERO, R );
106  FLA_Copyr( FLA_UPPER_TRIANGULAR, AT, R );
107 
108  // Reduce the matrix to bidiagonal form.
109  // Apply scalars to rotate elements on the superdiagonal to the real domain.
110  // Extract the diagonal and superdiagonal from A.
111  FLA_Bidiag_UT( R, T, S );
112  FLA_Bidiag_UT_realify( R, rL, rR );
114 
115  // Form V from right Householder vectors in upper triangle of R.
116  FLA_Bidiag_UT_form_V( R, S, V );
117 
118  // Form U in R.
119  FLA_Bidiag_UT_form_U( R, T, R );
120 
121  // Apply the realifying scalars in rL and rR to U and V, respectively.
122  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_CONJUGATE, rL, R );
123  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_NO_CONJUGATE, rR, V );
124 
125  // Perform a singular value decomposition on the bidiagonal matrix.
126  r_val = FLA_Bsvd_v_opt_var1( n_iter_max, d, e, G, H, R, V, b_alg );
127 
128  // Multiply R into U, storing the result in A and then copying back
129  // to U.
130  FLA_Part_1x2( U, &UL, &UR, n_A, FLA_LEFT );
131  FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
132  FLA_ONE, UL, R, FLA_ZERO, A );
133  FLA_Copy( A, UL );
134 
135  FLA_Obj_free( &R );
136  }
137 
138  // Copy the converged eigenvalues to the output vector.
139  FLA_Copy( d, s );
140 
141  // Sort the singular values and singular vectors in descending order.
142  FLA_Sort_svd( FLA_BACKWARD, s, U, V );
143 
144  // If the matrix was scaled, rescale the singular values.
145  if ( !FLA_Obj_equals( scale, FLA_ONE ) )
146  FLA_Inv_scal( scale, s );
147 
148  FLA_Obj_free( &scale );
149  FLA_Obj_free( &T );
150  FLA_Obj_free( &S );
151  FLA_Obj_free( &rL );
152  FLA_Obj_free( &rR );
153  FLA_Obj_free( &d );
154  FLA_Obj_free( &e );
155  FLA_Obj_free( &G );
156  FLA_Obj_free( &H );
157 
158  return r_val;
159 }
FLA_Error FLA_Bidiag_UT(FLA_Obj A, FLA_Obj TU, FLA_Obj TV)
Definition: FLA_Bidiag_UT.c:17
FLA_Error FLA_Bidiag_UT_form_V(FLA_Obj A, FLA_Obj S, FLA_Obj V)
Definition: FLA_Bidiag_UT_form_V.c:13
FLA_Error FLA_Bidiag_UT_extract_real_diagonals(FLA_Obj A, FLA_Obj d, FLA_Obj e)
Definition: FLA_Bidiag_UT_extract_real_diagonals.c:13
FLA_Error FLA_Bidiag_UT_form_U(FLA_Obj A, FLA_Obj T, FLA_Obj U)
Definition: FLA_Bidiag_UT_form_U.c:13
FLA_Error FLA_Bidiag_UT_create_T(FLA_Obj A, FLA_Obj *TU, FLA_Obj *TV)
Definition: FLA_Bidiag_UT_create_T.c:13
FLA_Error FLA_Bidiag_UT_realify(FLA_Obj A, FLA_Obj d, FLA_Obj e)
Definition: FLA_Bidiag_UT_realify.c:13
FLA_Error FLA_Bsvd_v_opt_var1(dim_t n_iter_max, FLA_Obj d, FLA_Obj e, FLA_Obj G, FLA_Obj H, FLA_Obj U, FLA_Obj V, dim_t b_alg)
Definition: FLA_Bsvd_v_opt_var1.c:15
FLA_Error FLA_QR_UT_form_Q(FLA_Obj A, FLA_Obj T, FLA_Obj Q)
Definition: FLA_QR_UT_form_Q.c:13
FLA_Error FLA_QR_UT_create_T(FLA_Obj A, FLA_Obj *T)
Definition: FLA_QR_UT_create_T.c:13
FLA_Error FLA_Svd_compute_scaling(FLA_Obj A, FLA_Obj sigma)
Definition: FLA_Svd_compute_scaling.c:13
FLA_Error FLA_Copy(FLA_Obj A, FLA_Obj B)
Definition: FLA_Copy.c:15
FLA_Error FLA_Scal(FLA_Obj alpha, FLA_Obj A)
Definition: FLA_Scal.c:15
FLA_Error FLA_Inv_scal(FLA_Obj alpha, FLA_Obj A)
Definition: FLA_Inv_scal.c:13
FLA_Error FLA_Copyr(FLA_Uplo uplo, FLA_Obj A, FLA_Obj B)
Definition: FLA_Copyr.c:15
FLA_Error FLA_Gemm(FLA_Trans transa, FLA_Trans transb, FLA_Obj alpha, FLA_Obj A, FLA_Obj B, FLA_Obj beta, FLA_Obj C)
Definition: FLA_Gemm.c:15
FLA_Obj FLA_ZERO
Definition: FLA_Init.c:20
FLA_Obj FLA_ONE
Definition: FLA_Init.c:18
FLA_Error FLA_QR_UT(FLA_Obj A, FLA_Obj T)
Definition: FLA_QR_UT.c:15
FLA_Datatype FLA_Obj_datatype_proj_to_complex(FLA_Obj A)
Definition: FLA_Query.c:37
dim_t FLA_Obj_width(FLA_Obj obj)
Definition: FLA_Query.c:123
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_Part_1x2(FLA_Obj A, FLA_Obj *A1, FLA_Obj *A2, dim_t nb, FLA_Side side)
Definition: FLA_View.c:110
FLA_Error FLA_Part_2x1(FLA_Obj A, FLA_Obj *A1, FLA_Obj *A2, dim_t mb, FLA_Side side)
Definition: FLA_View.c:76
FLA_Datatype FLA_Obj_datatype_proj_to_real(FLA_Obj A)
Definition: FLA_Query.c:23
dim_t FLA_Obj_length(FLA_Obj obj)
Definition: FLA_Query.c:116
FLA_Bool FLA_Obj_equals(FLA_Obj A, FLA_Obj B)
Definition: FLA_Query.c:507
dim_t FLA_Obj_min_dim(FLA_Obj obj)
Definition: FLA_Query.c:153
FLA_Error FLA_Obj_free(FLA_Obj *obj)
Definition: FLA_Obj.c:588
FLA_Datatype FLA_Obj_datatype(FLA_Obj obj)
Definition: FLA_Query.c:13
int FLA_Error
Definition: FLA_type_defs.h:47
int FLA_Datatype
Definition: FLA_type_defs.h:49
unsigned long dim_t
Definition: FLA_type_defs.h:71
FLA_Error FLA_Setr(FLA_Uplo uplo, FLA_Obj alpha, FLA_Obj A)
Definition: FLA_Setr.c:13
FLA_Error FLA_Apply_diag_matrix(FLA_Side side, FLA_Conj conj, FLA_Obj x, FLA_Obj A)
Definition: FLA_Apply_diag_matrix.c:13
FLA_Error FLA_Sort_svd(FLA_Direct direct, FLA_Obj s, FLA_Obj U, FLA_Obj V)
Definition: FLA_Sort_svd.c:13
Definition: FLA_type_defs.h:159

References FLA_Apply_diag_matrix(), FLA_Bidiag_UT(), FLA_Bidiag_UT_create_T(), FLA_Bidiag_UT_extract_real_diagonals(), FLA_Bidiag_UT_form_U(), FLA_Bidiag_UT_form_V(), FLA_Bidiag_UT_realify(), FLA_Bsvd_v_opt_var1(), FLA_Copy(), FLA_Copyr(), FLA_Gemm(), FLA_Inv_scal(), FLA_Obj_create(), FLA_Obj_datatype(), FLA_Obj_datatype_proj_to_complex(), FLA_Obj_datatype_proj_to_real(), FLA_Obj_equals(), FLA_Obj_free(), FLA_Obj_length(), FLA_Obj_min_dim(), FLA_Obj_width(), FLA_ONE, FLA_Part_1x2(), FLA_Part_2x1(), FLA_QR_UT(), FLA_QR_UT_create_T(), FLA_QR_UT_form_Q(), FLA_Scal(), FLA_Setr(), FLA_Sort_svd(), FLA_Svd_compute_scaling(), and FLA_ZERO.

◆ FLA_Svd_uv_unb_var2()

FLA_Error FLA_Svd_uv_unb_var2 ( dim_t  n_iter_max,
FLA_Obj  A,
FLA_Obj  s,
FLA_Obj  U,
FLA_Obj  V,
dim_t  k_accum,
dim_t  b_alg 
)
14 {
15  FLA_Error r_val = FLA_SUCCESS;
16  FLA_Datatype dt;
17  FLA_Datatype dt_real;
18  FLA_Datatype dt_comp;
19  FLA_Obj scale, T, S, rL, rR, d, e, G, H, RG, RH, W;
20  dim_t m_A, n_A;
21  dim_t min_m_n;
22  dim_t n_GH;
23  double crossover_ratio = 17.0 / 9.0;
24 
25  n_GH = k_accum;
26 
27  m_A = FLA_Obj_length( A );
28  n_A = FLA_Obj_width( A );
29  min_m_n = FLA_Obj_min_dim( A );
30  dt = FLA_Obj_datatype( A );
31  dt_real = FLA_Obj_datatype_proj_to_real( A );
32  dt_comp = FLA_Obj_datatype_proj_to_complex( A );
33 
34  // If the matrix is a scalar, then the SVD is easy.
35  if ( min_m_n == 1 )
36  {
37  FLA_Copy( A, s );
40 
41  return FLA_SUCCESS;
42  }
43 
44  // Create matrices to hold block Householder transformations.
45  FLA_Bidiag_UT_create_T( A, &T, &S );
46 
47  // Create vectors to hold the realifying scalars.
48  FLA_Obj_create( dt, min_m_n, 1, 0, 0, &rL );
49  FLA_Obj_create( dt, min_m_n, 1, 0, 0, &rR );
50 
51  // Create vectors to hold the diagonal and sub-diagonal.
52  FLA_Obj_create( dt_real, min_m_n, 1, 0, 0, &d );
53  FLA_Obj_create( dt_real, min_m_n-1, 1, 0, 0, &e );
54 
55  // Create matrices to hold the left and right Givens scalars.
56  FLA_Obj_create( dt_comp, min_m_n-1, n_GH, 0, 0, &G );
57  FLA_Obj_create( dt_comp, min_m_n-1, n_GH, 0, 0, &H );
58 
59  // Create matrices to hold the left and right Givens matrices.
60  FLA_Obj_create( dt_real, min_m_n, min_m_n, 0, 0, &RG );
61  FLA_Obj_create( dt_real, min_m_n, min_m_n, 0, 0, &RH );
62  FLA_Obj_create( dt, m_A, n_A, 0, 0, &W );
63 
64  // Create a real scaling factor.
65  FLA_Obj_create( dt_real, 1, 1, 0, 0, &scale );
66 
67  // Compute a scaling factor; If none is needed, sigma will be set to one.
68  FLA_Svd_compute_scaling( A, scale );
69 
70  // Scale the matrix if scale is non-unit.
71  if ( !FLA_Obj_equals( scale, FLA_ONE ) )
72  FLA_Scal( scale, A );
73 
74  if ( m_A >= n_A )
75  {
76  if ( m_A < crossover_ratio * n_A )
77  {
78  // Reduce the matrix to bidiagonal form.
79  // Apply scalars to rotate elements on the sub-diagonal to the real domain.
80  // Extract the diagonal and sub-diagonal from A.
81  FLA_Bidiag_UT( A, T, S );
82  FLA_Bidiag_UT_realify( A, rL, rR );
84 
85  // Form U and V.
86  FLA_Bidiag_UT_form_U( A, T, U );
87  FLA_Bidiag_UT_form_V( A, S, V );
88 
89  // Apply the realifying scalars in rL and rR to U and V, respectively.
90  {
91  FLA_Obj UL, UR;
92  FLA_Obj VL, VR;
93 
94  FLA_Part_1x2( U, &UL, &UR, min_m_n, FLA_LEFT );
95  FLA_Part_1x2( V, &VL, &VR, min_m_n, FLA_LEFT );
96 
97  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_CONJUGATE, rL, UL );
98  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_NO_CONJUGATE, rR, VL );
99  }
100 
101  // Perform a singular value decomposition on the bidiagonal matrix.
102  r_val = FLA_Bsvd_v_opt_var2( n_iter_max, d, e, G, H, RG, RH, W, U, V, b_alg );
103  }
104  else // if ( crossover_ratio * n_A <= m_A )
105  {
106  FLA_Obj TQ, R;
107  FLA_Obj AT,
108  AB;
109  FLA_Obj UL, UR;
110 
111  // Perform a QR factorization on A and form Q in U.
112  FLA_QR_UT_create_T( A, &TQ );
113  FLA_QR_UT( A, TQ );
114  FLA_QR_UT_form_Q( A, TQ, U );
115  FLA_Obj_free( &TQ );
116 
117  // Set the lower triangle of R to zero and then copy the upper
118  // triangle of A to R.
119  FLA_Part_2x1( A, &AT,
120  &AB, n_A, FLA_TOP );
121  FLA_Obj_create( dt, n_A, n_A, 0, 0, &R );
122  FLA_Setr( FLA_LOWER_TRIANGULAR, FLA_ZERO, R );
123  FLA_Copyr( FLA_UPPER_TRIANGULAR, AT, R );
124 
125  // Reduce the matrix to bidiagonal form.
126  // Apply scalars to rotate elements on the superdiagonal to the real domain.
127  // Extract the diagonal and superdiagonal from A.
128  FLA_Bidiag_UT( R, T, S );
129  FLA_Bidiag_UT_realify( R, rL, rR );
131 
132  // Form V from right Householder vectors in upper triangle of R.
133  FLA_Bidiag_UT_form_V( R, S, V );
134 
135  // Form U in R.
136  FLA_Bidiag_UT_form_U( R, T, R );
137 
138  // Apply the realifying scalars in rL and rR to U and V, respectively.
139  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_CONJUGATE, rL, R );
140  FLA_Apply_diag_matrix( FLA_RIGHT, FLA_NO_CONJUGATE, rR, V );
141 
142  // Perform a singular value decomposition on the bidiagonal matrix.
143  r_val = FLA_Bsvd_v_opt_var2( n_iter_max, d, e, G, H, RG, RH, W, R, V, b_alg );
144 
145  // Multiply R into U, storing the result in A and then copying back
146  // to U.
147  FLA_Part_1x2( U, &UL, &UR, n_A, FLA_LEFT );
148  FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
149  FLA_ONE, UL, R, FLA_ZERO, A );
150  FLA_Copy( A, UL );
151 
152  FLA_Obj_free( &R );
153  }
154  }
155  else // if ( m_A < n_A )
156  {
157  FLA_Check_error_code( FLA_NOT_YET_IMPLEMENTED );
158  }
159 
160  // Copy the converged eigenvalues to the output vector.
161  FLA_Copy( d, s );
162 
163  // Sort the singular values and singular vectors in descending order.
164  FLA_Sort_svd( FLA_BACKWARD, s, U, V );
165 
166  // If the matrix was scaled, rescale the singular values.
167  if ( !FLA_Obj_equals( scale, FLA_ONE ) )
168  FLA_Inv_scal( scale, s );
169 
170  FLA_Obj_free( &scale );
171  FLA_Obj_free( &T );
172  FLA_Obj_free( &S );
173  FLA_Obj_free( &rL );
174  FLA_Obj_free( &rR );
175  FLA_Obj_free( &d );
176  FLA_Obj_free( &e );
177  FLA_Obj_free( &G );
178  FLA_Obj_free( &H );
179  FLA_Obj_free( &RG );
180  FLA_Obj_free( &RH );
181  FLA_Obj_free( &W );
182 
183  return r_val;
184 }
FLA_Error FLA_Bsvd_v_opt_var2(dim_t n_iter_max, FLA_Obj d, FLA_Obj e, FLA_Obj G, FLA_Obj H, FLA_Obj RG, FLA_Obj RH, FLA_Obj W, FLA_Obj U, FLA_Obj V, dim_t b_alg)
Definition: FLA_Bsvd_v_opt_var2.c:13
FLA_Error FLA_Set_to_identity(FLA_Obj A)
Definition: FLA_Set_to_identity.c:13

References FLA_Apply_diag_matrix(), FLA_Bidiag_UT(), FLA_Bidiag_UT_create_T(), FLA_Bidiag_UT_extract_real_diagonals(), FLA_Bidiag_UT_form_U(), FLA_Bidiag_UT_form_V(), FLA_Bidiag_UT_realify(), FLA_Bsvd_v_opt_var2(), FLA_Copy(), FLA_Copyr(), FLA_Gemm(), FLA_Inv_scal(), FLA_Obj_create(), FLA_Obj_datatype(), FLA_Obj_datatype_proj_to_complex(), FLA_Obj_datatype_proj_to_real(), FLA_Obj_equals(), FLA_Obj_free(), FLA_Obj_length(), FLA_Obj_min_dim(), FLA_Obj_width(), FLA_ONE, FLA_Part_1x2(), FLA_Part_2x1(), FLA_QR_UT(), FLA_QR_UT_create_T(), FLA_QR_UT_form_Q(), FLA_Scal(), FLA_Set_to_identity(), FLA_Setr(), FLA_Sort_svd(), FLA_Svd_compute_scaling(), and FLA_ZERO.