QDP++
generic_blas_vcaxmby3.h
Go to the documentation of this file.
1// $Id: generic_blas_vcaxmby3.h,v 1.4 2009-09-15 20:48:42 bjoo Exp $
2
7
8#ifndef QDP_GENERIC_BLAS_VCAXMBY3
9#define QDP_GENERIC_BLAS_VCAXMBY3
10
11namespace QDP {
12// (Vector) out = (Complex) (*scalep) * (Vector) InScale - (Vector) Add
13inline
14void vcaxmby3(REAL* Out, REAL* ap, REAL* xp, REAL* bp, REAL* yp, int n_3vec)
15{
16 double a_r;
17 double a_i;
18 double b_r;
19 double b_i;
20
21 double x0r;
22 double x0i;
23
24 double x1r;
25 double x1i;
26
27 double x2r;
28 double x2i;
29
30 double y0r;
31 double y0i;
32
33 double y1r;
34 double y1i;
35
36 double y2r;
37 double y2i;
38
39 double z0r;
40 double z0i;
41
42 double z1r;
43 double z1i;
44
45 double z2r;
46 double z2i;
47
48 a_r =(double)(*ap);
49 a_i =(double)*(ap+1);
50
51 b_r =(double)(*bp);
52 b_i =(double)*(bp+1);
53
54 int index_x = 0;
55 int index_y = 0;
56 int index_z = 0;
57
58 int counter;
59
60 if( n_3vec > 0 ) {
61 // Prefetch whole vectors
62 x0r = (double)xp[index_x++];
63 x0i = (double)xp[index_x++];
64 y0r = (double)yp[index_y++];
65 y0i = (double)yp[index_y++];
66
67 x1r = (double)xp[index_x++];
68 x1i = (double)xp[index_x++];
69 y1r = (double)yp[index_y++];
70 y1i = (double)yp[index_y++];
71
72 x2r = (double)xp[index_x++];
73 x2i = (double)xp[index_x++];
74 y2r = (double)yp[index_y++];
75 y2i = (double)yp[index_y++];
76
77 int len = 4*n_3vec;
78
79 for( counter = 0; counter < len-1; counter++) {
80
81 z0r = a_r * x0r;
82 z0i = a_i * x0r;
83 x0r = (double)xp[index_x++];
84 z0r -= b_r * y0r;
85 z0i -= b_i * y0r;
86 y0r = (double)yp[index_y++];
87
88 z0r -= a_i * x0i;
89 z0i += a_r * x0i;
90 x0i = (double)xp[index_x++];
91 z0r += b_i * y0i;
92 z0i -= b_r * y0i;
93 y0i = (double)yp[index_y++];
94 Out[index_z++] = (REAL)z0r;
95 Out[index_z++] = (REAL)z0i;
96
97 z1r = a_r * x1r;
98 z1i = a_i * x1r;
99 x1r = (double)xp[index_x++];
100 z1r -= b_r * y1r;
101 z1i -= b_i * y1r;
102 y1r = (double)yp[index_y++];
103
104 z1r -= a_i * x1i;
105 z1i += a_r * x1i;
106 x1i = (double)xp[index_x++];
107 z1r += b_i * y1i;
108 z1i -= b_r * y1i;
109 y1i = (double)yp[index_y++];
110 Out[index_z++] = (REAL)z1r;
111 Out[index_z++] = (REAL)z1i;
112
113 z2r = a_r * x2r;
114 z2i = a_i * x2r;
115 x2r = (double)xp[index_x++];
116 z2r -= b_r * y2r;
117 z2i -= b_i * y2r;
118 y2r = (double)yp[index_y++];
119
120 z2r -= a_i * x2i;
121 z2i += a_r * x2i;
122 x2i = (double)xp[index_x++];
123 z2r += b_i * y2i;
124 z2i -= b_r * y2i;
125 y2i = (double)yp[index_y++];
126 Out[index_z++] = (REAL)z2r;
127 Out[index_z++] = (REAL)z2i;
128 }
129
130 z0r = a_r * x0r;
131 z0i = a_i * x0r;
132 z0r -= b_r * y0r;
133 z0i -= b_i * y0r;
134
135 z0r -= a_i * x0i;
136 z0i += a_r * x0i;
137 z0r += b_i * y0i;
138 Out[index_z++] = (REAL)z0r;
139 z0i -= b_r * y0i;
140 Out[index_z++] = (REAL)z0i;
141
142 z1r = a_r * x1r;
143 z1i = a_i * x1r;
144 z1r -= b_r * y1r;
145 z1i -= b_i * y1r;
146
147 z1r -= a_i * x1i;
148 z1i += a_r * x1i;
149 z1r += b_i * y1i;
150 Out[index_z++] = (REAL)z1r;
151 z1i -= b_r * y1i;
152 Out[index_z++] = (REAL)z1i;
153
154 z2r = a_r * x2r;
155 z2i = a_i * x2r;
156 z2r -= b_r * y2r;
157 z2i -= b_i * y2r;
158
159 z2r -= a_i * x2i;
160 z2i += a_r * x2i;
161 z2r += b_i * y2i;
162 Out[index_z++] = (REAL)z2r;
163 z2i -= b_r * y2i;
164 Out[index_z++] = (REAL)z2i;
165 }
166}
167
168
169} // namespace QDP;
170
171#endif // guard
REAL32 REAL
Yet another random number generator.
void vcaxmby3(REAL *Out, REAL *ap, REAL *xp, REAL *bp, REAL *yp, int n_3vec)