1 #include "blas_extended.h"
2 #include "blas_extended_private.h"
BLAS_zdot_z_c(enum blas_conj_type conj,int n,const void * alpha,const void * x,int incx,const void * beta,const void * y,int incy,void * r)3 void BLAS_zdot_z_c(enum blas_conj_type conj, int n, const void *alpha,
4 const void *x, int incx, const void *beta,
5 const void *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 void*
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_zdot_z_c";
46
47 int i, ix = 0, iy = 0;
48 double *r_i = (double *) r;
49 const double *x_i = (double *) x;
50 const float *y_i = (float *) y;
51 double *alpha_i = (double *) alpha;
52 double *beta_i = (double *) beta;
53 double x_ii[2];
54 float y_ii[2];
55 double r_v[2];
56 double prod[2];
57 double sum[2];
58 double tmp1[2];
59 double 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 incy *= 2;
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[0] = y_i[iy];
92 y_ii[1] = y_i[iy + 1];
93 x_ii[1] = -x_ii[1];
94 {
95 prod[0] = (double) x_ii[0] * y_ii[0] - (double) x_ii[1] * y_ii[1];
96 prod[1] = (double) x_ii[0] * y_ii[1] + (double) x_ii[1] * y_ii[0];
97 } /* prod = x[i]*y[i] */
98 sum[0] = sum[0] + prod[0];
99 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
100 ix += incx;
101 iy += incy;
102 } /* endfor */
103 } else {
104 /* do not conjugate */
105
106 for (i = 0; i < n; ++i) {
107 x_ii[0] = x_i[ix];
108 x_ii[1] = x_i[ix + 1];
109 y_ii[0] = y_i[iy];
110 y_ii[1] = y_i[iy + 1];
111
112 {
113 prod[0] = (double) x_ii[0] * y_ii[0] - (double) x_ii[1] * y_ii[1];
114 prod[1] = (double) x_ii[0] * y_ii[1] + (double) x_ii[1] * y_ii[0];
115 } /* prod = x[i]*y[i] */
116 sum[0] = sum[0] + prod[0];
117 sum[1] = sum[1] + prod[1]; /* sum = sum+prod */
118 ix += incx;
119 iy += incy;
120 } /* endfor */
121 }
122
123 {
124 tmp1[0] = (double) sum[0] * alpha_i[0] - (double) sum[1] * alpha_i[1];
125 tmp1[1] = (double) sum[0] * alpha_i[1] + (double) sum[1] * alpha_i[0];
126 } /* tmp1 = sum*alpha */
127 {
128 tmp2[0] = (double) r_v[0] * beta_i[0] - (double) r_v[1] * beta_i[1];
129 tmp2[1] = (double) r_v[0] * beta_i[1] + (double) r_v[1] * beta_i[0];
130 } /* tmp2 = r*beta */
131 tmp1[0] = tmp1[0] + tmp2[0];
132 tmp1[1] = tmp1[1] + tmp2[1]; /* tmp1 = tmp1+tmp2 */
133 ((double *) r)[0] = tmp1[0];
134 ((double *) r)[1] = tmp1[1]; /* r = tmp1 */
135
136
137
138 }
139