1 #include <stdio.h>
2 #include "blas_extended.h"
3
che_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,enum blas_side_type side,int n,void * a,int lda,void * a_vec,int row)4 void che_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
5 enum blas_side_type side, int n, void *a, int lda,
6 void *a_vec, int row)
7
8 /* Copies the given row of matrix a into the supplied vector. */
9 {
10 int conj_flag;
11 int i;
12 int ai, incai1, incai2;
13 int vi, incvi = 1;
14
15 float a_elem[2];
16 const float *a_i = (float *) a;
17 float *a_vec_i = (float *) a_vec;
18
19 if ((side == blas_left_side && uplo == blas_upper) ||
20 (side == blas_right_side && uplo == blas_lower))
21 conj_flag = 1;
22 else
23 conj_flag = 0;
24
25 if ((order == blas_colmajor && uplo == blas_upper) ||
26 (order == blas_rowmajor && uplo == blas_lower)) {
27 incai1 = 1;
28 incai2 = lda;
29 ai = lda * row;
30 } else {
31 incai1 = lda;
32 incai2 = 1;
33 ai = row;
34 }
35
36 incai1 *= 2;
37 incai2 *= 2;
38 ai *= 2;
39 incvi *= 2;
40
41 for (i = 0, vi = 0; i < row; i++, ai += incai1, vi += incvi) {
42 a_elem[0] = a_i[ai];
43 a_elem[1] = a_i[ai + 1];
44 if (conj_flag == 1) {
45 a_elem[1] = -a_elem[1];
46 }
47 a_vec_i[vi] = a_elem[0];
48 a_vec_i[vi + 1] = a_elem[1];
49 }
50 for (; i < n; i++, ai += incai2, vi += incvi) {
51 a_elem[0] = a_i[ai];
52 a_elem[1] = a_i[ai + 1];
53 if (conj_flag == 0) {
54 a_elem[1] = -a_elem[1];
55 }
56 a_vec_i[vi] = a_elem[0];
57 a_vec_i[vi + 1] = a_elem[1];
58 }
59 }
zhe_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,enum blas_side_type side,int n,void * a,int lda,void * a_vec,int row)60 void zhe_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
61 enum blas_side_type side, int n, void *a, int lda,
62 void *a_vec, int row)
63
64 /* Copies the given row of matrix a into the supplied vector. */
65 {
66 int conj_flag;
67 int i;
68 int ai, incai1, incai2;
69 int vi, incvi = 1;
70
71 double a_elem[2];
72 const double *a_i = (double *) a;
73 double *a_vec_i = (double *) a_vec;
74
75 if ((side == blas_left_side && uplo == blas_upper) ||
76 (side == blas_right_side && uplo == blas_lower))
77 conj_flag = 1;
78 else
79 conj_flag = 0;
80
81 if ((order == blas_colmajor && uplo == blas_upper) ||
82 (order == blas_rowmajor && uplo == blas_lower)) {
83 incai1 = 1;
84 incai2 = lda;
85 ai = lda * row;
86 } else {
87 incai1 = lda;
88 incai2 = 1;
89 ai = row;
90 }
91
92 incai1 *= 2;
93 incai2 *= 2;
94 ai *= 2;
95 incvi *= 2;
96
97 for (i = 0, vi = 0; i < row; i++, ai += incai1, vi += incvi) {
98 a_elem[0] = a_i[ai];
99 a_elem[1] = a_i[ai + 1];
100 if (conj_flag == 1) {
101 a_elem[1] = -a_elem[1];
102 }
103 a_vec_i[vi] = a_elem[0];
104 a_vec_i[vi + 1] = a_elem[1];
105 }
106 for (; i < n; i++, ai += incai2, vi += incvi) {
107 a_elem[0] = a_i[ai];
108 a_elem[1] = a_i[ai + 1];
109 if (conj_flag == 0) {
110 a_elem[1] = -a_elem[1];
111 }
112 a_vec_i[vi] = a_elem[0];
113 a_vec_i[vi + 1] = a_elem[1];
114 }
115 }
116
che_print_matrix(void * a,int n,int lda,enum blas_order_type order,enum blas_uplo_type uplo)117 void che_print_matrix(void *a, int n, int lda,
118 enum blas_order_type order, enum blas_uplo_type uplo)
119 {
120 int ai, aij;
121 int incai, incaij1, incaij2;
122 int i, j;
123
124 float a_elem[2];
125 const float *a_i = (float *) a;
126
127 if ((order == blas_colmajor && uplo == blas_upper) ||
128 (order == blas_rowmajor && uplo == blas_lower)) {
129 incai = lda;
130 incaij1 = 1;
131 incaij2 = lda;
132 } else {
133 incai = 1;
134 incaij1 = lda;
135 incaij2 = 1;
136 }
137
138 incai *= 2;
139 incaij1 *= 2;
140 incaij2 *= 2;
141
142 for (i = 0, ai = 0; i < n; i++, ai += incai) {
143
144 for (j = 0, aij = ai; j < i; j++, aij += incaij1) {
145 a_elem[0] = a_i[aij];
146 a_elem[1] = a_i[aij + 1];
147 if (uplo == blas_upper) {
148 a_elem[1] = -a_elem[1];
149 }
150 printf("(%16.8e, %16.8e)", a_elem[0], a_elem[1]);
151 }
152 a_elem[0] = a_i[ai];
153 a_elem[1] = 0.0;
154 printf("(%16.8e, %16.8e)", a_elem[0], a_elem[1]);
155 j++;
156 aij += incaij2;
157 for (; j < n; j++, aij += incaij2) {
158 a_elem[0] = a_i[aij];
159 a_elem[1] = a_i[aij + 1];
160 if (uplo == blas_lower) {
161 a_elem[1] = -a_elem[1];
162 }
163 printf("(%16.8e, %16.8e)", a_elem[0], a_elem[1]);
164 }
165 printf("\n");
166 }
167
168 }
zhe_print_matrix(void * a,int n,int lda,enum blas_order_type order,enum blas_uplo_type uplo)169 void zhe_print_matrix(void *a, int n, int lda,
170 enum blas_order_type order, enum blas_uplo_type uplo)
171 {
172 int ai, aij;
173 int incai, incaij1, incaij2;
174 int i, j;
175
176 double a_elem[2];
177 const double *a_i = (double *) a;
178
179 if ((order == blas_colmajor && uplo == blas_upper) ||
180 (order == blas_rowmajor && uplo == blas_lower)) {
181 incai = lda;
182 incaij1 = 1;
183 incaij2 = lda;
184 } else {
185 incai = 1;
186 incaij1 = lda;
187 incaij2 = 1;
188 }
189
190 incai *= 2;
191 incaij1 *= 2;
192 incaij2 *= 2;
193
194 for (i = 0, ai = 0; i < n; i++, ai += incai) {
195
196 for (j = 0, aij = ai; j < i; j++, aij += incaij1) {
197 a_elem[0] = a_i[aij];
198 a_elem[1] = a_i[aij + 1];
199 if (uplo == blas_upper) {
200 a_elem[1] = -a_elem[1];
201 }
202 printf("(%24.16e, %24.16e)", a_elem[0], a_elem[1]);
203 }
204 a_elem[0] = a_i[ai];
205 a_elem[1] = 0.0;
206 printf("(%24.16e, %24.16e)", a_elem[0], a_elem[1]);
207 j++;
208 aij += incaij2;
209 for (; j < n; j++, aij += incaij2) {
210 a_elem[0] = a_i[aij];
211 a_elem[1] = a_i[aij + 1];
212 if (uplo == blas_lower) {
213 a_elem[1] = -a_elem[1];
214 }
215 printf("(%24.16e, %24.16e)", a_elem[0], a_elem[1]);
216 }
217 printf("\n");
218 }
219
220 }
221
sskew_commit_row(enum blas_order_type order,enum blas_uplo_type uplo,enum blas_side_type side,int n,float * a,int lda,float * a_vec,int row)222 void sskew_commit_row(enum blas_order_type order, enum blas_uplo_type uplo,
223 enum blas_side_type side, int n, float *a, int lda,
224 float *a_vec, int row)
225
226 /*
227 * Copies the given vector into the given row of symmetric matrix a.
228 */
229 {
230 int conj_flag;
231 int i;
232 int ai, incai1, incai2;
233 int vi, incvi = 1;
234
235 float a_elem;
236 float *a_i = a;
237 const float *a_vec_i = a_vec;
238
239 if ((side == blas_left_side && uplo == blas_upper) ||
240 (side == blas_right_side && uplo == blas_lower))
241 conj_flag = 1;
242 else
243 conj_flag = 0;
244
245 if ((order == blas_colmajor && uplo == blas_upper) ||
246 (order == blas_rowmajor && uplo == blas_lower)) {
247 incai1 = 1;
248 incai2 = lda;
249 ai = lda * row;
250 } else {
251 incai1 = lda;
252 incai2 = 1;
253 ai = row;
254 }
255
256
257
258
259
260
261 for (i = 0, vi = 0; i < row; i++, ai += incai1, vi += incvi) {
262 a_elem = a_vec_i[vi];
263 if (conj_flag == 1) {
264 a_elem = -a_elem;
265 }
266 a_i[ai] = a_elem;
267 }
268
269 for (; i < n; i++, ai += incai2, vi += incvi) {
270 a_elem = a_vec_i[vi];
271 if (conj_flag == 0) {
272 a_elem = -a_elem;
273 }
274 a_i[ai] = a_elem;
275 }
276 }
dskew_commit_row(enum blas_order_type order,enum blas_uplo_type uplo,enum blas_side_type side,int n,double * a,int lda,double * a_vec,int row)277 void dskew_commit_row(enum blas_order_type order, enum blas_uplo_type uplo,
278 enum blas_side_type side, int n, double *a, int lda,
279 double *a_vec, int row)
280
281 /*
282 * Copies the given vector into the given row of symmetric matrix a.
283 */
284 {
285 int conj_flag;
286 int i;
287 int ai, incai1, incai2;
288 int vi, incvi = 1;
289
290 double a_elem;
291 double *a_i = a;
292 const double *a_vec_i = a_vec;
293
294 if ((side == blas_left_side && uplo == blas_upper) ||
295 (side == blas_right_side && uplo == blas_lower))
296 conj_flag = 1;
297 else
298 conj_flag = 0;
299
300 if ((order == blas_colmajor && uplo == blas_upper) ||
301 (order == blas_rowmajor && uplo == blas_lower)) {
302 incai1 = 1;
303 incai2 = lda;
304 ai = lda * row;
305 } else {
306 incai1 = lda;
307 incai2 = 1;
308 ai = row;
309 }
310
311
312
313
314
315
316 for (i = 0, vi = 0; i < row; i++, ai += incai1, vi += incvi) {
317 a_elem = a_vec_i[vi];
318 if (conj_flag == 1) {
319 a_elem = -a_elem;
320 }
321 a_i[ai] = a_elem;
322 }
323
324 for (; i < n; i++, ai += incai2, vi += incvi) {
325 a_elem = a_vec_i[vi];
326 if (conj_flag == 0) {
327 a_elem = -a_elem;
328 }
329 a_i[ai] = a_elem;
330 }
331 }
332
sskew_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,enum blas_side_type side,int n,float * a,int lda,float * a_vec,int row)333 void sskew_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
334 enum blas_side_type side, int n, float *a, int lda,
335 float *a_vec, int row)
336
337 /*
338 * Copies the given row of skew matrix a into the supplied vector.
339 */
340 {
341 int conj_flag;
342 int i;
343 int ai, incai1, incai2;
344 int vi, incvi = 1;
345
346 float a_elem;
347 const float *a_i = a;
348 float *a_vec_i = a_vec;
349
350 if ((side == blas_left_side && uplo == blas_upper) ||
351 (side == blas_right_side && uplo == blas_lower))
352 conj_flag = 1;
353 else
354 conj_flag = 0;
355
356 if ((order == blas_colmajor && uplo == blas_upper) ||
357 (order == blas_rowmajor && uplo == blas_lower)) {
358 incai1 = 1;
359 incai2 = lda;
360 ai = lda * row;
361 } else {
362 incai1 = lda;
363 incai2 = 1;
364 ai = row;
365 }
366
367
368
369
370
371
372 for (i = 0, vi = 0; i < row; i++, ai += incai1, vi += incvi) {
373 a_elem = a_i[ai];
374 if (conj_flag == 1) {
375 a_elem = -a_elem;
376 }
377 a_vec_i[vi] = a_elem;
378 }
379
380 for (; i < n; i++, ai += incai2, vi += incvi) {
381 a_elem = a_i[ai];
382 if (conj_flag == 0) {
383 a_elem = -a_elem;
384 }
385 a_vec_i[vi] = a_elem;
386 }
387 }
dskew_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,enum blas_side_type side,int n,double * a,int lda,double * a_vec,int row)388 void dskew_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
389 enum blas_side_type side, int n, double *a, int lda,
390 double *a_vec, int row)
391
392 /*
393 * Copies the given row of skew matrix a into the supplied vector.
394 */
395 {
396 int conj_flag;
397 int i;
398 int ai, incai1, incai2;
399 int vi, incvi = 1;
400
401 double a_elem;
402 const double *a_i = a;
403 double *a_vec_i = a_vec;
404
405 if ((side == blas_left_side && uplo == blas_upper) ||
406 (side == blas_right_side && uplo == blas_lower))
407 conj_flag = 1;
408 else
409 conj_flag = 0;
410
411 if ((order == blas_colmajor && uplo == blas_upper) ||
412 (order == blas_rowmajor && uplo == blas_lower)) {
413 incai1 = 1;
414 incai2 = lda;
415 ai = lda * row;
416 } else {
417 incai1 = lda;
418 incai2 = 1;
419 ai = row;
420 }
421
422
423
424
425
426
427 for (i = 0, vi = 0; i < row; i++, ai += incai1, vi += incvi) {
428 a_elem = a_i[ai];
429 if (conj_flag == 1) {
430 a_elem = -a_elem;
431 }
432 a_vec_i[vi] = a_elem;
433 }
434
435 for (; i < n; i++, ai += incai2, vi += incvi) {
436 a_elem = a_i[ai];
437 if (conj_flag == 0) {
438 a_elem = -a_elem;
439 }
440 a_vec_i[vi] = a_elem;
441 }
442 }
443