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