1 #include "blas_extended.h"
2 #include "blas_extended_private.h"
BLAS_cdot_c_s(enum blas_conj_type conj,int n,const void * alpha,const void * x,int incx,const void * beta,const float * y,int incy,void * r)3 void BLAS_cdot_c_s(enum blas_conj_type conj, int n, const void *alpha,
4 const void *x, int incx, const void *beta,
5 const float *y, int incy, void *r)
6
7 /*
8 * Purpose
9 * =======
10 *
11 * This routine computes the inner product:
12 *
13 * r <- beta * r + alpha * SUM_{i=0, n-1} x[i] * y[i].
14 *
15 * Arguments
16 * =========
17 *
18 * conj (input) enum blas_conj_type
19 * When x and y are complex vectors, specifies whether vector
20 * components x[i] are used unconjugated or conjugated.
21 *
22 * n (input) int
23 * The length of vectors x and y.
24 *
25 * alpha (input) const void*
26 *
27 * x (input) const void*
28 * Array of length n.
29 *
30 * incx (input) int
31 * The stride used to access components x[i].
32 *
33 * beta (input) const void*
34 *
35 * y (input) const float*
36 * Array of length n.
37 *
38 * incy (input) int
39 * The stride used to access components y[i].
40 *
41 * r (input/output) void*
42 *
43 */
44 {
45 static const char routine_name[] = "BLAS_cdot_c_s";
46
47 int i, ix = 0, iy = 0;
48 float *r_i = (float *) r;
49 const float *x_i = (float *) x;
50 const float *y_i = y;
51 float *alpha_i = (float *) alpha;
52 float *beta_i = (float *) beta;
53 float x_ii[2];
54 float y_ii;
55 float r_v[2];
56 float prod[2];
57 float sum[2];
58 float tmp1[2];
59 float tmp2[2];
60
61
62 /* Test the input parameters. */
63 if (n < 0)
64 BLAS_error(routine_name, -2, n, NULL);
65 else if (incx == 0)
66 BLAS_error(routine_name, -5, incx, NULL);
67 else if (incy == 0)
68 BLAS_error(routine_name, -8, incy, NULL);
69
70 /* Immediate return. */
71 if (((beta_i[0] == 1.0 && beta_i[1] == 0.0))
72 && (n == 0 || (alpha_i[0] == 0.0 && alpha_i[1] == 0.0)))
73 return;
74
75
76
77 r_v[0] = r_i[0];
78 r_v[1] = r_i[0 + 1];
79 sum[0] = sum[1] = 0.0;
80 incx *= 2;
81
82 if (incx < 0)
83 ix = (-n + 1) * incx;
84 if (incy < 0)
85 iy = (-n + 1) * incy;
86
87 if (conj == blas_conj) {
88 for (i = 0; i < n; ++i) {
89 x_ii[0] = x_i[ix];
90 x_ii[1] = x_i[ix + 1];
91 y_ii = y_i[iy];
92 x_ii[1] = -x_ii[1];
93 {
94 prod[0] = x_ii[0] * y_ii;
95 prod[1] = x_ii[1] * y_ii;
96 } /* prod = x[i]*y[i] */
97 sum[0] = sum[0] + prod[0];
98 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
99 ix += incx;
100 iy += incy;
101 } /* endfor */
102 } else {
103 /* do not conjugate */
104
105 for (i = 0; i < n; ++i) {
106 x_ii[0] = x_i[ix];
107 x_ii[1] = x_i[ix + 1];
108 y_ii = y_i[iy];
109
110 {
111 prod[0] = x_ii[0] * y_ii;
112 prod[1] = x_ii[1] * y_ii;
113 } /* prod = x[i]*y[i] */
114 sum[0] = sum[0] + prod[0];
115 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
116 ix += incx;
117 iy += incy;
118 } /* endfor */
119 }
120
121 {
122 tmp1[0] = sum[0] * alpha_i[0] - sum[1] * alpha_i[1];
123 tmp1[1] = sum[0] * alpha_i[1] + sum[1] * alpha_i[0];
124 }
125 /* tmp1 = sum*alpha */
126 {
127 tmp2[0] = r_v[0] * beta_i[0] - r_v[1] * beta_i[1];
128 tmp2[1] = r_v[0] * beta_i[1] + r_v[1] * beta_i[0];
129 }
130 /* tmp2 = r*beta */
131 tmp1[0] = tmp1[0] + tmp2[0];
132 tmp1[1] = tmp1[1] + tmp2[1]; /* tmp1 = tmp1+tmp2 */
133 ((float *) r)[0] = tmp1[0];
134 ((float *) r)[1] = tmp1[1]; /* r = tmp1 */
135
136
137
138 }
139