1 /*
2
3 Copyright (C) 2014, The University of Texas at Austin
4
5 This file is part of libflame and is available under the 3-Clause
6 BSD license, which can be found in the LICENSE file at the top-level
7 directory, or at http://opensource.org/licenses/BSD-3-Clause
8
9 */
10
11 #include "FLAME.h"
12
FLA_Bidiag_UT_extract_diagonals(FLA_Obj A,FLA_Obj d,FLA_Obj e)13 FLA_Error FLA_Bidiag_UT_extract_diagonals( FLA_Obj A, FLA_Obj d, FLA_Obj e )
14 {
15 FLA_Error r_val = FLA_SUCCESS;
16
17 if ( FLA_Check_error_level() >= FLA_MIN_ERROR_CHECKING )
18 FLA_Bidiag_UT_extract_diagonals_check( A, d, e );
19
20 if ( FLA_Obj_length( A ) >= FLA_Obj_width( A ) )
21 r_val = FLA_Bidiag_UT_u_extract_diagonals( A, d, e );
22 else
23 r_val = FLA_Bidiag_UT_l_extract_diagonals( A, d, e );
24
25 return r_val;
26 }
27
28
29
FLA_Bidiag_UT_u_extract_diagonals(FLA_Obj A,FLA_Obj d,FLA_Obj e)30 FLA_Error FLA_Bidiag_UT_u_extract_diagonals( FLA_Obj A, FLA_Obj d, FLA_Obj e )
31 {
32 FLA_Datatype datatype;
33 int n_A;
34 int rs_A, cs_A;
35 int inc_d;
36 int inc_e;
37 int i;
38
39 datatype = FLA_Obj_datatype( A );
40
41 n_A = FLA_Obj_width( A );
42
43 rs_A = FLA_Obj_row_stride( A );
44 cs_A = FLA_Obj_col_stride( A );
45
46 inc_d = FLA_Obj_vector_inc( d );
47
48 if ( n_A != 1 )
49 inc_e = FLA_Obj_vector_inc( e );
50 else
51 inc_e = 0;
52
53 switch ( datatype )
54 {
55 case FLA_FLOAT:
56 {
57 float* buff_A = FLA_FLOAT_PTR( A );
58 float* buff_d = FLA_FLOAT_PTR( d );
59 float* buff_e = ( n_A != 1 ? FLA_FLOAT_PTR( e ) : NULL );
60
61 for ( i = 0; i < n_A; ++i )
62 {
63 float* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
64 float* a12t_l = buff_A + (i+1)*cs_A + (i )*rs_A;
65 float* delta1 = buff_d + (i )*inc_d;
66 float* epsilon1 = buff_e + (i )*inc_e;
67
68 int n_ahead = n_A - i - 1;
69
70 // delta1 = alpha11;
71 *delta1 = *alpha11;
72
73 // epsilon1 = a12t_l;
74 if ( n_ahead > 0 )
75 *epsilon1 = *a12t_l;
76 }
77
78 break;
79 }
80
81 case FLA_DOUBLE:
82 {
83 double* buff_A = FLA_DOUBLE_PTR( A );
84 double* buff_d = FLA_DOUBLE_PTR( d );
85 double* buff_e = ( n_A != 1 ? FLA_DOUBLE_PTR( e ) : NULL );
86
87 for ( i = 0; i < n_A; ++i )
88 {
89 double* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
90 double* a12t_l = buff_A + (i+1)*cs_A + (i )*rs_A;
91 double* delta1 = buff_d + (i )*inc_d;
92 double* epsilon1 = buff_e + (i )*inc_e;
93
94 int n_ahead = n_A - i - 1;
95
96 // delta1 = alpha11;
97 *delta1 = *alpha11;
98
99 // epsilon1 = a12t_l;
100 if ( n_ahead > 0 )
101 *epsilon1 = *a12t_l;
102 }
103
104 break;
105 }
106
107 case FLA_COMPLEX:
108 {
109 scomplex* buff_A = FLA_COMPLEX_PTR( A );
110 scomplex* buff_d = FLA_COMPLEX_PTR( d );
111 scomplex* buff_e = ( n_A != 1 ? FLA_COMPLEX_PTR( e ) : NULL );
112
113 for ( i = 0; i < n_A; ++i )
114 {
115 scomplex* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
116 scomplex* a12t_l = buff_A + (i+1)*cs_A + (i )*rs_A;
117 scomplex* delta1 = buff_d + (i )*inc_d;
118 scomplex* epsilon1 = buff_e + (i )*inc_e;
119
120 int n_ahead = n_A - i - 1;
121
122 // delta1 = alpha11;
123 *delta1 = *alpha11;
124
125 // epsilon1 = a12t_l;
126 if ( n_ahead > 0 )
127 *epsilon1 = *a12t_l;
128 }
129
130 break;
131 }
132
133 case FLA_DOUBLE_COMPLEX:
134 {
135 dcomplex* buff_A = FLA_DOUBLE_COMPLEX_PTR( A );
136 dcomplex* buff_d = FLA_DOUBLE_COMPLEX_PTR( d );
137 dcomplex* buff_e = ( n_A != 1 ? FLA_DOUBLE_COMPLEX_PTR( e ) : NULL );
138
139 for ( i = 0; i < n_A; ++i )
140 {
141 dcomplex* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
142 dcomplex* a12t_l = buff_A + (i+1)*cs_A + (i )*rs_A;
143 dcomplex* delta1 = buff_d + (i )*inc_d;
144 dcomplex* epsilon1 = buff_e + (i )*inc_e;
145
146 int n_ahead = n_A - i - 1;
147
148 // delta1 = alpha11;
149 *delta1 = *alpha11;
150
151 // epsilon1 = a12t_l;
152 if ( n_ahead > 0 )
153 *epsilon1 = *a12t_l;
154 }
155
156 break;
157 }
158 }
159
160 return FLA_SUCCESS;
161 }
162
163
164
FLA_Bidiag_UT_l_extract_diagonals(FLA_Obj A,FLA_Obj d,FLA_Obj e)165 FLA_Error FLA_Bidiag_UT_l_extract_diagonals( FLA_Obj A, FLA_Obj d, FLA_Obj e )
166 {
167 FLA_Datatype datatype;
168 int m_A;
169 int rs_A, cs_A;
170 int inc_d;
171 int inc_e;
172 int i;
173
174 datatype = FLA_Obj_datatype( A );
175
176 m_A = FLA_Obj_length( A );
177
178 rs_A = FLA_Obj_row_stride( A );
179 cs_A = FLA_Obj_col_stride( A );
180
181 inc_d = FLA_Obj_vector_inc( d );
182
183 if ( m_A != 1 )
184 inc_e = FLA_Obj_vector_inc( e );
185 else
186 inc_e = 0;
187
188 switch ( datatype )
189 {
190 case FLA_FLOAT:
191 {
192 float* buff_A = FLA_FLOAT_PTR( A );
193 float* buff_d = FLA_FLOAT_PTR( d );
194 float* buff_e = ( m_A != 1 ? FLA_FLOAT_PTR( e ) : NULL );
195
196 for ( i = 0; i < m_A; ++i )
197 {
198 float* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
199 float* a21_t = buff_A + (i )*cs_A + (i+1)*rs_A;
200 float* delta1 = buff_d + (i )*inc_d;
201 float* epsilon1 = buff_e + (i )*inc_e;
202
203 int m_ahead = m_A - i - 1;
204
205 // delta1 = alpha11;
206 *delta1 = *alpha11;
207
208 // epsilon1 = a21_t;
209 if ( m_ahead > 0 )
210 *epsilon1 = *a21_t;
211 }
212
213 break;
214 }
215
216 case FLA_DOUBLE:
217 {
218 double* buff_A = FLA_DOUBLE_PTR( A );
219 double* buff_d = FLA_DOUBLE_PTR( d );
220 double* buff_e = ( m_A != 1 ? FLA_DOUBLE_PTR( e ) : NULL );
221
222 for ( i = 0; i < m_A; ++i )
223 {
224 double* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
225 double* a21_t = buff_A + (i )*cs_A + (i+1)*rs_A;
226 double* delta1 = buff_d + (i )*inc_d;
227 double* epsilon1 = buff_e + (i )*inc_e;
228
229 int m_ahead = m_A - i - 1;
230
231 // delta1 = alpha11;
232 *delta1 = *alpha11;
233
234 // epsilon1 = a21_t;
235 if ( m_ahead > 0 )
236 *epsilon1 = *a21_t;
237 }
238
239 break;
240 }
241
242 case FLA_COMPLEX:
243 {
244 scomplex* buff_A = FLA_COMPLEX_PTR( A );
245 scomplex* buff_d = FLA_COMPLEX_PTR( d );
246 scomplex* buff_e = ( m_A != 1 ? FLA_COMPLEX_PTR( e ) : NULL );
247
248 for ( i = 0; i < m_A; ++i )
249 {
250 scomplex* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
251 scomplex* a21_t = buff_A + (i )*cs_A + (i+1)*rs_A;
252 scomplex* delta1 = buff_d + (i )*inc_d;
253 scomplex* epsilon1 = buff_e + (i )*inc_e;
254
255 int m_ahead = m_A - i - 1;
256
257 // delta1 = alpha11;
258 *delta1 = *alpha11;
259
260 // epsilon1 = a21_t;
261 if ( m_ahead > 0 )
262 *epsilon1 = *a21_t;
263 }
264
265 break;
266 }
267
268 case FLA_DOUBLE_COMPLEX:
269 {
270 dcomplex* buff_A = FLA_DOUBLE_COMPLEX_PTR( A );
271 dcomplex* buff_d = FLA_DOUBLE_COMPLEX_PTR( d );
272 dcomplex* buff_e = ( m_A != 1 ? FLA_DOUBLE_COMPLEX_PTR( e ) : NULL );
273
274 for ( i = 0; i < m_A; ++i )
275 {
276 dcomplex* alpha11 = buff_A + (i )*cs_A + (i )*rs_A;
277 dcomplex* a21_t = buff_A + (i )*cs_A + (i+1)*rs_A;
278 dcomplex* delta1 = buff_d + (i )*inc_d;
279 dcomplex* epsilon1 = buff_e + (i )*inc_e;
280
281 int m_ahead = m_A - i - 1;
282
283 // delta1 = alpha11;
284 *delta1 = *alpha11;
285
286 // epsilon1 = a21_t;
287 if ( m_ahead > 0 )
288 *epsilon1 = *a21_t;
289 }
290
291 break;
292 }
293 }
294
295 return FLA_SUCCESS;
296 }
297
298