1
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include "blas_extended.h"
5 #include "blas_extended_private.h"
6 #include "blas_extended_test.h"
7
sskew_commit_row_hbmv(enum blas_order_type order,enum blas_uplo_type uplo,int n,float * a,int k,int lda,float * a_vec,int row)8 void sskew_commit_row_hbmv(enum blas_order_type order,
9 enum blas_uplo_type uplo, int n, float *a, int k,
10 int lda, float *a_vec, int row)
11
12 /*
13 * Copies the given vector into the given row of symmetric matrix a.
14 */
15 {
16
17 int ai;
18 int i;
19 int incai1;
20 int incai2;
21 int vi;
22 int incvi = 1;
23 int loopmax;
24
25 float a_elem;
26 float *a_i = a;
27 const float *a_vec_i = a_vec;
28
29 if (order == blas_colmajor) {
30 if (uplo == blas_upper) {
31 incai1 = 1;
32 incai2 = lda - 1;
33 } else {
34 incai1 = lda - 1;
35 incai2 = 1;
36 }
37 } else {
38 if (uplo == blas_upper) {
39 incai1 = lda - 1;
40 incai2 = 1;
41 } else {
42 incai1 = 1;
43 incai2 = lda - 1;
44 }
45 }
46
47 ai = 0;
48 if ((uplo == blas_upper && order == blas_colmajor) ||
49 (uplo == blas_lower && order == blas_rowmajor)) {
50 /* starting place */
51 ai = lda * row + ((row < k) ? (k - row) : 0);
52 } else {
53 /* starting place */
54 ai = (row > k) ? (k + lda * (row - k)) : row;
55 }
56
57
58
59
60
61
62 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
63 /* this is a wasteful loop but important */
64 }
65
66 for (; i < row; i++, ai += incai1, vi += incvi) {
67 a_elem = a_vec_i[vi];
68 if (uplo == blas_upper) {
69 a_elem = -a_elem;;
70 }
71 a_i[ai] = a_elem;
72 }
73
74 loopmax = MIN(row + k + 1, n);
75 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
76 a_elem = a_vec_i[vi];
77 if (uplo == blas_lower) {
78 a_elem = -a_elem;;
79 }
80 a_i[ai] = a_elem;
81 }
82 }
dskew_commit_row_hbmv(enum blas_order_type order,enum blas_uplo_type uplo,int n,double * a,int k,int lda,double * a_vec,int row)83 void dskew_commit_row_hbmv(enum blas_order_type order,
84 enum blas_uplo_type uplo, int n, double *a, int k,
85 int lda, double *a_vec, int row)
86
87 /*
88 * Copies the given vector into the given row of symmetric matrix a.
89 */
90 {
91
92 int ai;
93 int i;
94 int incai1;
95 int incai2;
96 int vi;
97 int incvi = 1;
98 int loopmax;
99
100 double a_elem;
101 double *a_i = a;
102 const double *a_vec_i = a_vec;
103
104 if (order == blas_colmajor) {
105 if (uplo == blas_upper) {
106 incai1 = 1;
107 incai2 = lda - 1;
108 } else {
109 incai1 = lda - 1;
110 incai2 = 1;
111 }
112 } else {
113 if (uplo == blas_upper) {
114 incai1 = lda - 1;
115 incai2 = 1;
116 } else {
117 incai1 = 1;
118 incai2 = lda - 1;
119 }
120 }
121
122 ai = 0;
123 if ((uplo == blas_upper && order == blas_colmajor) ||
124 (uplo == blas_lower && order == blas_rowmajor)) {
125 /* starting place */
126 ai = lda * row + ((row < k) ? (k - row) : 0);
127 } else {
128 /* starting place */
129 ai = (row > k) ? (k + lda * (row - k)) : row;
130 }
131
132
133
134
135
136
137 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
138 /* this is a wasteful loop but important */
139 }
140
141 for (; i < row; i++, ai += incai1, vi += incvi) {
142 a_elem = a_vec_i[vi];
143 if (uplo == blas_upper) {
144 a_elem = -a_elem;;
145 }
146 a_i[ai] = a_elem;
147 }
148
149 loopmax = MIN(row + k + 1, n);
150 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
151 a_elem = a_vec_i[vi];
152 if (uplo == blas_lower) {
153 a_elem = -a_elem;;
154 }
155 a_i[ai] = a_elem;
156 }
157 }
158
chbmv_commit_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,int k,int lda,void * a_vec,int row)159 void chbmv_commit_row(enum blas_order_type order, enum blas_uplo_type uplo,
160 int n, void *a, int k, int lda, void *a_vec, int row)
161
162 /*
163 * Copies the given vector into the given row of symmetric matrix a.
164 */
165 {
166
167 int ai;
168 int i;
169 int incai1;
170 int incai2;
171 int vi;
172 int incvi = 1;
173 int loopmax;
174
175 float a_elem[2];
176 float *a_i = (float *) a;
177 const float *a_vec_i = (float *) a_vec;
178
179 if (order == blas_colmajor) {
180 if (uplo == blas_upper) {
181 incai1 = 1;
182 incai2 = lda - 1;
183 } else {
184 incai1 = lda - 1;
185 incai2 = 1;
186 }
187 } else {
188 if (uplo == blas_upper) {
189 incai1 = lda - 1;
190 incai2 = 1;
191 } else {
192 incai1 = 1;
193 incai2 = lda - 1;
194 }
195 }
196
197 ai = 0;
198 if ((uplo == blas_upper && order == blas_colmajor) ||
199 (uplo == blas_lower && order == blas_rowmajor)) {
200 /* starting place */
201 ai = lda * row + ((row < k) ? (k - row) : 0);
202 } else {
203 /* starting place */
204 ai = (row > k) ? (k + lda * (row - k)) : row;
205 }
206
207 incai1 *= 2;
208 incai2 *= 2;
209 ai *= 2;
210 incvi *= 2;
211
212 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
213 /* this is a wasteful loop but important */
214 }
215
216 for (; i < row; i++, ai += incai1, vi += incvi) {
217 a_elem[0] = a_vec_i[vi];
218 a_elem[1] = a_vec_i[vi + 1];
219 if (uplo == blas_upper) {
220 a_elem[1] = -a_elem[1];;
221 }
222 a_i[ai] = a_elem[0];
223 a_i[ai + 1] = a_elem[1];
224 }
225
226 loopmax = MIN(row + k + 1, n);
227 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
228 a_elem[0] = a_vec_i[vi];
229 a_elem[1] = a_vec_i[vi + 1];
230 if (uplo == blas_lower) {
231 a_elem[1] = -a_elem[1];;
232 }
233 a_i[ai] = a_elem[0];
234 a_i[ai + 1] = a_elem[1];
235 }
236 }
zhbmv_commit_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,int k,int lda,void * a_vec,int row)237 void zhbmv_commit_row(enum blas_order_type order, enum blas_uplo_type uplo,
238 int n, void *a, int k, int lda, void *a_vec, int row)
239
240 /*
241 * Copies the given vector into the given row of symmetric matrix a.
242 */
243 {
244
245 int ai;
246 int i;
247 int incai1;
248 int incai2;
249 int vi;
250 int incvi = 1;
251 int loopmax;
252
253 double a_elem[2];
254 double *a_i = (double *) a;
255 const double *a_vec_i = (double *) a_vec;
256
257 if (order == blas_colmajor) {
258 if (uplo == blas_upper) {
259 incai1 = 1;
260 incai2 = lda - 1;
261 } else {
262 incai1 = lda - 1;
263 incai2 = 1;
264 }
265 } else {
266 if (uplo == blas_upper) {
267 incai1 = lda - 1;
268 incai2 = 1;
269 } else {
270 incai1 = 1;
271 incai2 = lda - 1;
272 }
273 }
274
275 ai = 0;
276 if ((uplo == blas_upper && order == blas_colmajor) ||
277 (uplo == blas_lower && order == blas_rowmajor)) {
278 /* starting place */
279 ai = lda * row + ((row < k) ? (k - row) : 0);
280 } else {
281 /* starting place */
282 ai = (row > k) ? (k + lda * (row - k)) : row;
283 }
284
285 incai1 *= 2;
286 incai2 *= 2;
287 ai *= 2;
288 incvi *= 2;
289
290 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
291 /* this is a wasteful loop but important */
292 }
293
294 for (; i < row; i++, ai += incai1, vi += incvi) {
295 a_elem[0] = a_vec_i[vi];
296 a_elem[1] = a_vec_i[vi + 1];
297 if (uplo == blas_upper) {
298 a_elem[1] = -a_elem[1];;
299 }
300 a_i[ai] = a_elem[0];
301 a_i[ai + 1] = a_elem[1];
302 }
303
304 loopmax = MIN(row + k + 1, n);
305 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
306 a_elem[0] = a_vec_i[vi];
307 a_elem[1] = a_vec_i[vi + 1];
308 if (uplo == blas_lower) {
309 a_elem[1] = -a_elem[1];;
310 }
311 a_i[ai] = a_elem[0];
312 a_i[ai + 1] = a_elem[1];
313 }
314 }
315
sskew_copy_row_hbmv(enum blas_order_type order,enum blas_uplo_type uplo,int n,float * a,int k,int lda,float * a_vec,int row)316 void sskew_copy_row_hbmv(enum blas_order_type order, enum blas_uplo_type uplo,
317 int n, float *a, int k, int lda,
318 float *a_vec, int row)
319
320 /*
321 * Copies the given row of matrix a into the supplied vector.
322 */
323 {
324
325 int ai;
326 int i;
327 int incai1;
328 int incai2;
329 int vi;
330 int incvi = 1;
331 int loopmax;
332
333 float a_elem;
334 const float *a_i = a;
335 float *a_vec_i = a_vec;
336
337
338 if (order == blas_colmajor) {
339 if (uplo == blas_upper) {
340 incai1 = 1;
341 incai2 = lda - 1;
342 } else {
343 incai1 = lda - 1;
344 incai2 = 1;
345 }
346 } else {
347 if (uplo == blas_upper) {
348 incai1 = lda - 1;
349 incai2 = 1;
350 } else {
351 incai1 = 1;
352 incai2 = lda - 1;
353 }
354 }
355
356 ai = 0;
357 if ((uplo == blas_upper && order == blas_colmajor) ||
358 (uplo == blas_lower && order == blas_rowmajor)) {
359 /* starting place */
360 ai = lda * row + ((row < k) ? (k - row) : 0);
361 } else {
362 /* starting place */
363 ai = (row > k) ? (k + lda * (row - k)) : row;
364 }
365
366
367
368
369
370
371 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
372 a_vec_i[vi] = 0.0;
373 }
374
375 for (; i < row; i++, ai += incai1, vi += incvi) {
376 a_elem = a_i[ai];
377 if (uplo == blas_upper) {
378 a_elem = -a_elem;;
379 }
380 a_vec_i[vi] = a_elem;
381 }
382
383 loopmax = MIN(row + k + 1, n);
384 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
385 a_elem = a_i[ai];
386 if (i == row) {
387
388 }
389 if (uplo == blas_lower) {
390 a_elem = -a_elem;;
391 }
392 a_vec_i[vi] = a_elem;
393 }
394 for (; i < n; i++, vi += incvi) {
395 a_vec_i[vi] = 0.0;
396 }
397
398
399 }
dskew_copy_row_hbmv(enum blas_order_type order,enum blas_uplo_type uplo,int n,double * a,int k,int lda,double * a_vec,int row)400 void dskew_copy_row_hbmv(enum blas_order_type order, enum blas_uplo_type uplo,
401 int n, double *a, int k, int lda,
402 double *a_vec, int row)
403
404 /*
405 * Copies the given row of matrix a into the supplied vector.
406 */
407 {
408
409 int ai;
410 int i;
411 int incai1;
412 int incai2;
413 int vi;
414 int incvi = 1;
415 int loopmax;
416
417 double a_elem;
418 const double *a_i = a;
419 double *a_vec_i = a_vec;
420
421
422 if (order == blas_colmajor) {
423 if (uplo == blas_upper) {
424 incai1 = 1;
425 incai2 = lda - 1;
426 } else {
427 incai1 = lda - 1;
428 incai2 = 1;
429 }
430 } else {
431 if (uplo == blas_upper) {
432 incai1 = lda - 1;
433 incai2 = 1;
434 } else {
435 incai1 = 1;
436 incai2 = lda - 1;
437 }
438 }
439
440 ai = 0;
441 if ((uplo == blas_upper && order == blas_colmajor) ||
442 (uplo == blas_lower && order == blas_rowmajor)) {
443 /* starting place */
444 ai = lda * row + ((row < k) ? (k - row) : 0);
445 } else {
446 /* starting place */
447 ai = (row > k) ? (k + lda * (row - k)) : row;
448 }
449
450
451
452
453
454
455 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
456 a_vec_i[vi] = 0.0;
457 }
458
459 for (; i < row; i++, ai += incai1, vi += incvi) {
460 a_elem = a_i[ai];
461 if (uplo == blas_upper) {
462 a_elem = -a_elem;;
463 }
464 a_vec_i[vi] = a_elem;
465 }
466
467 loopmax = MIN(row + k + 1, n);
468 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
469 a_elem = a_i[ai];
470 if (i == row) {
471
472 }
473 if (uplo == blas_lower) {
474 a_elem = -a_elem;;
475 }
476 a_vec_i[vi] = a_elem;
477 }
478 for (; i < n; i++, vi += incvi) {
479 a_vec_i[vi] = 0.0;
480 }
481
482
483 }
484
chbmv_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,int k,int lda,void * a_vec,int row)485 void chbmv_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
486 int n, void *a, int k, int lda, void *a_vec, int row)
487
488 /*
489 * Copies the given row of matrix a into the supplied vector.
490 */
491 {
492
493 int ai;
494 int i;
495 int incai1;
496 int incai2;
497 int vi;
498 int incvi = 1;
499 int loopmax;
500
501 float a_elem[2];
502 const float *a_i = (float *) a;
503 float *a_vec_i = (float *) a_vec;
504
505
506 if (order == blas_colmajor) {
507 if (uplo == blas_upper) {
508 incai1 = 1;
509 incai2 = lda - 1;
510 } else {
511 incai1 = lda - 1;
512 incai2 = 1;
513 }
514 } else {
515 if (uplo == blas_upper) {
516 incai1 = lda - 1;
517 incai2 = 1;
518 } else {
519 incai1 = 1;
520 incai2 = lda - 1;
521 }
522 }
523
524 ai = 0;
525 if ((uplo == blas_upper && order == blas_colmajor) ||
526 (uplo == blas_lower && order == blas_rowmajor)) {
527 /* starting place */
528 ai = lda * row + ((row < k) ? (k - row) : 0);
529 } else {
530 /* starting place */
531 ai = (row > k) ? (k + lda * (row - k)) : row;
532 }
533
534 incai1 *= 2;
535 incai2 *= 2;
536 ai *= 2;
537 incvi *= 2;
538
539 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
540 a_vec_i[vi] = 0.0;
541 a_vec_i[vi + 1] = 0.0;
542 }
543
544 for (; i < row; i++, ai += incai1, vi += incvi) {
545 a_elem[0] = a_i[ai];
546 a_elem[1] = a_i[ai + 1];
547 if (uplo == blas_upper) {
548 a_elem[1] = -a_elem[1];;
549 }
550 a_vec_i[vi] = a_elem[0];
551 a_vec_i[vi + 1] = a_elem[1];
552 }
553
554 loopmax = MIN(row + k + 1, n);
555 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
556 a_elem[0] = a_i[ai];
557 a_elem[1] = a_i[ai + 1];
558 if (i == row) {
559 a_elem[1] = 0.0;
560 }
561 if (uplo == blas_lower) {
562 a_elem[1] = -a_elem[1];;
563 }
564 a_vec_i[vi] = a_elem[0];
565 a_vec_i[vi + 1] = a_elem[1];
566 }
567 for (; i < n; i++, vi += incvi) {
568 a_vec_i[vi] = 0.0;
569 a_vec_i[vi + 1] = 0.0;
570 }
571 }
zhbmv_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,int k,int lda,void * a_vec,int row)572 void zhbmv_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
573 int n, void *a, int k, int lda, void *a_vec, int row)
574
575 /*
576 * Copies the given row of matrix a into the supplied vector.
577 */
578 {
579
580 int ai;
581 int i;
582 int incai1;
583 int incai2;
584 int vi;
585 int incvi = 1;
586 int loopmax;
587
588 double a_elem[2];
589 const double *a_i = (double *) a;
590 double *a_vec_i = (double *) a_vec;
591
592
593 if (order == blas_colmajor) {
594 if (uplo == blas_upper) {
595 incai1 = 1;
596 incai2 = lda - 1;
597 } else {
598 incai1 = lda - 1;
599 incai2 = 1;
600 }
601 } else {
602 if (uplo == blas_upper) {
603 incai1 = lda - 1;
604 incai2 = 1;
605 } else {
606 incai1 = 1;
607 incai2 = lda - 1;
608 }
609 }
610
611 ai = 0;
612 if ((uplo == blas_upper && order == blas_colmajor) ||
613 (uplo == blas_lower && order == blas_rowmajor)) {
614 /* starting place */
615 ai = lda * row + ((row < k) ? (k - row) : 0);
616 } else {
617 /* starting place */
618 ai = (row > k) ? (k + lda * (row - k)) : row;
619 }
620
621 incai1 *= 2;
622 incai2 *= 2;
623 ai *= 2;
624 incvi *= 2;
625
626 for (i = 0, vi = 0; i < row - k; i++, vi += incvi) {
627 a_vec_i[vi] = 0.0;
628 a_vec_i[vi + 1] = 0.0;
629 }
630
631 for (; i < row; i++, ai += incai1, vi += incvi) {
632 a_elem[0] = a_i[ai];
633 a_elem[1] = a_i[ai + 1];
634 if (uplo == blas_upper) {
635 a_elem[1] = -a_elem[1];;
636 }
637 a_vec_i[vi] = a_elem[0];
638 a_vec_i[vi + 1] = a_elem[1];
639 }
640
641 loopmax = MIN(row + k + 1, n);
642 for (; i < loopmax; i++, ai += incai2, vi += incvi) {
643 a_elem[0] = a_i[ai];
644 a_elem[1] = a_i[ai + 1];
645 if (i == row) {
646 a_elem[1] = 0.0;
647 }
648 if (uplo == blas_lower) {
649 a_elem[1] = -a_elem[1];;
650 }
651 a_vec_i[vi] = a_elem[0];
652 a_vec_i[vi + 1] = a_elem[1];
653 }
654 for (; i < n; i++, vi += incvi) {
655 a_vec_i[vi] = 0.0;
656 a_vec_i[vi + 1] = 0.0;
657 }
658 }
659
cprint_hbmv_matrix(void * a,int n,int k,int lda,enum blas_order_type order,enum blas_uplo_type uplo)660 void cprint_hbmv_matrix(void *a, int n, int k, int lda,
661 enum blas_order_type order, enum blas_uplo_type uplo)
662 {
663
664 int row;
665 float *x;
666 x = (float *) blas_malloc(n * sizeof(float) * 2);
667 if (n > 0 && x == NULL) {
668 BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
669 }
670
671 for (row = 0; row < n; row++) {
672 chbmv_copy_row(order, uplo, n, a, k, lda, x, row);
673 cprint_vector(x, n, 1, NULL);
674 }
675 printf("\n");
676 blas_free(x);
677
678 }
zprint_hbmv_matrix(void * a,int n,int k,int lda,enum blas_order_type order,enum blas_uplo_type uplo)679 void zprint_hbmv_matrix(void *a, int n, int k, int lda,
680 enum blas_order_type order, enum blas_uplo_type uplo)
681 {
682
683 int row;
684 double *x;
685 x = (double *) blas_malloc(n * sizeof(double) * 2);
686 if (n > 0 && x == NULL) {
687 BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
688 }
689
690 for (row = 0; row < n; row++) {
691 zhbmv_copy_row(order, uplo, n, a, k, lda, x, row);
692 zprint_vector(x, n, 1, NULL);
693 }
694 printf("\n");
695 blas_free(x);
696
697 }
698