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