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