← Back to branch summary

~ubuntu-branches/ubuntu/vivid/atlas/vivid

« back to all changes in this revision

Viewing changes to tune/blas/gemv/CASES/ATL_gemvN_16x2_1.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot, Sylvestre Ledru, Sébastien Villemot
  • Date: 2013-06-11 15:58:16 UTC
  • mfrom: (1.1.4) (25 sid)
  • mto: This revision was merged to the branch mainline in revision 26.
  • Revision ID: package-import@ubuntu.com-20130611155816-8xeeiziu1iml040c
Tags: 3.10.1-1
http://bugs.debian.org/609287
http://bugs.debian.org/602524
http://bugs.debian.org/701921
http://bugs.debian.org/666203
http://bugs.debian.org/701068
http://bugs.debian.org/697431
[ Sylvestre Ledru ]
* New upstream release (Closes: #609287)

[ Sébastien Villemot ]
* Provide architectural defaults (i.e. precomputed timings) for all
  release archs (except armel and mips for the time being, due to slow
  porterboxes). This will make the package build much faster and should
  eliminate transient build failures due to excessive variance in the
  timings.
* Move symlinks for lib{cblas,f77blas,atlas,lapack_atlas} out of the
  libblas.so.3 alternative and make them always present, so that
  software relying on these libs do not break when another alternative
  is selected for BLAS
* ATLAS now has improved ARM support with native asm constructs. This required
  the following tunes:
  + armel-is-v4t.diff: new patch, prevents FTBFS on armel; otherwise,
    ATLAS uses asm constructs too recent for the platform (armel is only v4t)
  + debian/rules: on armhf, define the ATL_ARM_HARDFP flag; otherwise the asm
    constructs use the soft-float ABI for passing floating points
  + on armhf, ensure that -mfloat-abi=softfp and -mcpu=vfpv3 flags are never
    used; this is implemented via a patch (armhf.diff) and by the use of fixed
    archdefs
* The generic package is now built without multi-threading, because otherwise
  the package fails to build on some single-processor machines (this required
  the introduction of a patch: fix-non-threaded-build.diff). As a side effect,
  the build of the custom package gracefully handles non-threaded
  builds. (Closes: #602524)
* Add libblas.a as slave in the libblas.so alternative (Closes: #701921)
* Add symlinks for lib{f77blas,atlas}.a in /usr/lib (Closes: #666203)
* Modify shlibs file of libatlas3-base, such that packages using
  libblas/liblapack depend on any BLAS/LAPACK alternative, while packages
  depending on ATLAS-specific libraries (e.g. libatlas.so) depend specifically
  on libatlas3-base.
* corei1.diff: remove patch, applied upstream
* Use my @debian.org email address
* Remove obsolete DM-Upload-Allowed flag
* Switch VCS to git
* Remove Conflicts/Replaces against pre-squeeze packages
* libatlas-base-dev now provides libblas.so, as libblas-dev
* No longer use -Wa,--noexecstack in CFLAGS, it makes the package FTBFS
* Do not use POWER3 arch for powerpcspe port (Closes: #701068)
* Bump to debhelper compat level 9
* README.Debian: mention that devscripts is needed to compile the custom
  package (Closes: #697431)
* Bump Standards-Version to 3.9.4. As a consequence, add Built-Using
  fields because the package embeds stuff from liblapack-pic

Show diffs side-by-side

added added

removed removed

Lines of Context:
tune/blas/gemv/CASES/ATL_gemvN_16x2_1.c
1
 
/*
2
 
 *             Automatically Tuned Linear Algebra Software v3.8.4
3
 
 *                    (C) Copyright 1999 R. Clint Whaley
4
 
 *
5
 
 * Redistribution and use in source and binary forms, with or without
6
 
 * modification, are permitted provided that the following conditions
7
 
 * are met:
8
 
 *   1. Redistributions of source code must retain the above copyright
9
 
 *      notice, this list of conditions and the following disclaimer.
10
 
 *   2. Redistributions in binary form must reproduce the above copyright
11
 
 *      notice, this list of conditions, and the following disclaimer in the
12
 
 *      documentation and/or other materials provided with the distribution.
13
 
 *   3. The name of the ATLAS group or the names of its contributers may
14
 
 *      not be used to endorse or promote products derived from this
15
 
 *      software without specific written permission.
16
 
 *
17
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
18
 
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
19
 
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
20
 
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
21
 
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22
 
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23
 
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24
 
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25
 
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26
 
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27
 
 * POSSIBILITY OF SUCH DAMAGE.
28
 
 *
29
 
 */
30
 
 
31
 
#include "atlas_misc.h"
32
 
#include "atlas_level1.h"
33
 
#include "atlas_level2.h"
34
 
#include "atlas_lvl2.h"
35
 
 
36
 
 
37
 
static void gemvMlt8(const int M, const int N, const TYPE *A, const int lda,
38
 
                     const TYPE *X, const SCALAR beta, TYPE *Y);
39
 
static void gemvNle4(const int M, const int N, const TYPE *A, const int lda,
40
 
                     const TYPE *X, const SCALAR beta, TYPE *Y);
41
 
#ifdef BETA0
42
 
   #define Yget(y_, yp_, bet_) (y_) = ATL_rzero
43
 
#elif defined BETAX
44
 
   #define Yget(y_, yp_, bet_) (y_) = (yp_) * (bet_)
45
 
#else
46
 
   #define Yget(y_, yp_, bet_) (y_) = (yp_)
47
 
#endif
48
 
static void gemvN32x4(const int M, const int N, const TYPE *A, const int lda,
49
 
                      const TYPE *x, const SCALAR beta0, TYPE *y)
50
 
/*
51
 
 * rank-4 daxpy based NoTrans gemv
52
 
 */
53
 
{
54
 
   const int M16 = (M>>4)<<4;
55
 
   TYPE *stY = y + M16 - 32;
56
 
   const TYPE *A0 = A, *A1 = A+lda, *A2 = A1 + lda, *A3 = A2 + lda;
57
 
   register TYPE z0, z1, z2, z3, z4, z5, z6, z7;
58
 
   register TYPE y0, y1, y2, y3, y4, y5, y6, y7;
59
 
   const register TYPE x0 = *x, x1 = x[1], x2 = x[2], x3 = x[3];
60
 
   #ifdef BETAX
61
 
      const register TYPE beta = beta0;
62
 
   #else
63
 
      #define beta beta0
64
 
   #endif
65
 
 
66
 
   ATL_assert(N == 4);
67
 
   if (M16 >= 32)
68
 
   {
69
 
      #ifdef BETA0
70
 
         y0 = y1 = y2 = y3 = y4 = y5 = y6 = y7 = ATL_rzero;
71
 
      #else
72
 
         y0 = *y;   y1 = y[1]; y2 = y[2]; y3 = y[3];
73
 
         y4 = y[4]; y5 = y[5]; y6 = y[6]; y7 = y[7];
74
 
         #ifdef BETAX
75
 
            y0 *= beta; y1 *= beta; y2 *= beta; y3 *= beta;
76
 
            y4 *= beta; y5 *= beta; y6 *= beta; y7 *= beta;
77
 
         #endif
78
 
      #endif
79
 
      y0 += x0 * *A0;   Yget(z0, y[8], beta);
80
 
      y1 += x1 * A1[1];
81
 
      y2 += x2 * A2[2]; Yget(z1, y[9], beta);
82
 
      y3 += x3 * A3[3];
83
 
      y4 += x0 * A0[4]; Yget(z2, y[10], beta);
84
 
      y5 += x1 * A1[5];
85
 
      y6 += x2 * A2[6]; Yget(z3, y[11], beta);
86
 
      y7 += x3 * A3[7];
87
 
 
88
 
      y0 += x1 * *A1;   Yget(z4, y[12], beta);
89
 
      y1 += x2 * A2[1];
90
 
      y2 += x3 * A3[2]; Yget(z5, y[13], beta);
91
 
      y3 += x0 * A0[3];
92
 
      y4 += x1 * A1[4]; Yget(z6, y[14], beta);
93
 
      y5 += x2 * A2[5];
94
 
      y6 += x3 * A3[6]; Yget(z7, y[15], beta);
95
 
      y7 += x0 * A0[7];
96
 
 
97
 
      y0 += x2 * *A2;
98
 
      y1 += x3 * A3[1];
99
 
      y2 += x0 * A0[2];
100
 
      y3 += x1 * A1[3];
101
 
      y4 += x2 * A2[4];
102
 
      y5 += x3 * A3[5];
103
 
      y6 += x0 * A0[6];
104
 
      y7 += x1 * A1[7];
105
 
 
106
 
      y0 += x3 * *A3;
107
 
      y1 += x0 * A0[1];
108
 
      y2 += x1 * A1[2];
109
 
      y3 += x2 * A2[3];
110
 
      y4 += x3 * A3[4];
111
 
      y5 += x0 * A0[5];
112
 
      y6 += x1 * A1[6];
113
 
      y7 += x2 * A2[7];
114
 
 
115
 
      z0 += x0 * A0[8]; *y = y0;
116
 
      z1 += x1 * A1[9];
117
 
      z2 += x2 * A2[10]; y[1] = y1;
118
 
      z3 += x3 * A3[11];
119
 
      z4 += x0 * A0[12]; y[2] = y2;
120
 
      z5 += x1 * A1[13];
121
 
      z6 += x2 * A2[14]; y[3] = y3;
122
 
      z7 += x3 * A3[15];
123
 
 
124
 
      z0 += x1 * A1[8]; y[4] = y4;
125
 
      z1 += x2 * A2[9];
126
 
      z2 += x3 * A3[10]; y[5] = y5;
127
 
      z3 += x0 * A0[11];
128
 
      z4 += x1 * A1[12]; y[6] = y6;
129
 
      z5 += x2 * A2[13];
130
 
      z6 += x3 * A3[14]; y[7] = y7;
131
 
      z7 += x0 * A0[15];
132
 
 
133
 
      z0 += x2 * A2[8];  Yget(y0, y[16], beta);
134
 
      z1 += x3 * A3[ 9];
135
 
      z2 += x0 * A0[10]; Yget(y1, y[17], beta);
136
 
      z3 += x1 * A1[11];
137
 
      z4 += x2 * A2[12]; Yget(y2, y[18], beta);
138
 
      z5 += x3 * A3[13];
139
 
      z6 += x0 * A0[14]; Yget(y3, y[19], beta);
140
 
      z7 += x1 * A1[15];
141
 
 
142
 
      z0 += x3 * A3[8];  Yget(y4, y[20], beta);
143
 
      z1 += x0 * A0[9];
144
 
      z2 += x1 * A1[10]; Yget(y5, y[21], beta);
145
 
      z3 += x2 * A2[11];
146
 
      z4 += x3 * A3[12]; Yget(y6, y[22], beta); A3 += 16;
147
 
      z5 += x0 * A0[13];                        A0 += 16;
148
 
      z6 += x1 * A1[14]; Yget(y7, y[23], beta); A1 += 16;
149
 
      z7 += x2 * A2[15];                        A2 += 16;
150
 
      if (M16 != 32)
151
 
      {
152
 
         do
153
 
         {
154
 
            y0 += x0 * *A0;   y[8] = z0;
155
 
            y1 += x1 * A1[1];
156
 
            y2 += x2 * A2[2]; y[9] = z1;
157
 
            y3 += x3 * A3[3];
158
 
            y4 += x0 * A0[4]; y[10] = z2;
159
 
            y5 += x1 * A1[5];
160
 
            y6 += x2 * A2[6]; y[11] = z3;
161
 
            y7 += x3 * A3[7];
162
 
 
163
 
            y0 += x1 * *A1;   y[12] = z4;
164
 
            y1 += x2 * A2[1];
165
 
            y2 += x3 * A3[2]; y[13] = z5;
166
 
            y3 += x0 * A0[3];
167
 
            y4 += x1 * A1[4]; y[14] = z6;
168
 
            y5 += x2 * A2[5];
169
 
            y6 += x3 * A3[6]; y[15] = z7; y += 16;
170
 
            y7 += x0 * A0[7];
171
 
 
172
 
            y0 += x2 * *A2;   Yget(z0, y[8], beta);
173
 
            y1 += x3 * A3[1];
174
 
            y2 += x0 * A0[2]; Yget(z1, y[9], beta);
175
 
            y3 += x1 * A1[3];
176
 
            y4 += x2 * A2[4]; Yget(z2, y[10], beta);
177
 
            y5 += x3 * A3[5];
178
 
            y6 += x0 * A0[6]; Yget(z3, y[11], beta);
179
 
            y7 += x1 * A1[7];
180
 
 
181
 
            y0 += x3 * *A3;   Yget(z4, y[12], beta);
182
 
            y1 += x0 * A0[1];
183
 
            y2 += x1 * A1[2]; Yget(z5, y[13], beta);
184
 
            y3 += x2 * A2[3];
185
 
            y4 += x3 * A3[4]; Yget(z6, y[14], beta);
186
 
            y5 += x0 * A0[5];
187
 
            y6 += x1 * A1[6]; Yget(z7, y[15], beta);
188
 
            y7 += x2 * A2[7];
189
 
 
190
 
            z0 += x0 * A0[8];  *y = y0;
191
 
            z1 += x1 * A1[9];
192
 
            z2 += x2 * A2[10]; y[1] = y1;
193
 
            z3 += x3 * A3[11];
194
 
            z4 += x0 * A0[12]; y[2] = y2;
195
 
            z5 += x1 * A1[13];
196
 
            z6 += x2 * A2[14]; y[3] = y3;
197
 
            z7 += x3 * A3[15];
198
 
 
199
 
            z0 += x1 * A1[8]; y[4] = y4;
200
 
            z1 += x2 * A2[9];
201
 
            z2 += x3 * A3[10]; y[5] = y5;
202
 
            z3 += x0 * A0[11];
203
 
            z4 += x1 * A1[12]; y[6] = y6;
204
 
            z5 += x2 * A2[13];
205
 
            z6 += x3 * A3[14]; y[7] = y7;
206
 
            z7 += x0 * A0[15];
207
 
 
208
 
            z0 += x2 * A2[8];  Yget(y0, y[16], beta);
209
 
            z1 += x3 * A3[ 9];
210
 
            z2 += x0 * A0[10]; Yget(y1, y[17], beta);
211
 
            z3 += x1 * A1[11];
212
 
            z4 += x2 * A2[12]; Yget(y2, y[18], beta);
213
 
            z5 += x3 * A3[13];
214
 
            z6 += x0 * A0[14]; Yget(y3, y[19], beta);
215
 
            z7 += x1 * A1[15];
216
 
 
217
 
            z0 += x3 * A3[8];  Yget(y4, y[20], beta);
218
 
            z1 += x0 * A0[9];
219
 
            z2 += x1 * A1[10]; Yget(y5, y[21], beta);
220
 
            z3 += x2 * A2[11];
221
 
            z4 += x3 * A3[12]; Yget(y6, y[22], beta);  A3 += 16;
222
 
            z5 += x0 * A0[13];                         A0 += 16;
223
 
            z6 += x1 * A1[14]; Yget(y7, y[23], beta);  A1 += 16;
224
 
            z7 += x2 * A2[15];                         A2 += 16;
225
 
         }
226
 
         while (y != stY);
227
 
      }
228
 
      y0 += x0 * *A0;   y[8] = z0;
229
 
      y1 += x1 * A1[1];
230
 
      y2 += x2 * A2[2]; y[9] = z1;
231
 
      y3 += x3 * A3[3];
232
 
      y4 += x0 * A0[4]; y[10] = z2;
233
 
      y5 += x1 * A1[5];
234
 
      y6 += x2 * A2[6]; y[11] = z3;
235
 
      y7 += x3 * A3[7];
236
 
 
237
 
      y0 += x1 * *A1;   y[12] = z4;
238
 
      y1 += x2 * A2[1];
239
 
      y2 += x3 * A3[2]; y[13] = z5;
240
 
      y3 += x0 * A0[3];
241
 
      y4 += x1 * A1[4]; y[14] = z6;
242
 
      y5 += x2 * A2[5];
243
 
      y6 += x3 * A3[6]; y[15] = z7; y += 16;
244
 
      y7 += x0 * A0[7];
245
 
 
246
 
      y0 += x2 * *A2;   Yget(z0, y[8], beta);
247
 
      y1 += x3 * A3[1];
248
 
      y2 += x0 * A0[2]; Yget(z1, y[9], beta);
249
 
      y3 += x1 * A1[3];
250
 
      y4 += x2 * A2[4]; Yget(z2, y[10], beta);
251
 
      y5 += x3 * A3[5];
252
 
      y6 += x0 * A0[6]; Yget(z3, y[11], beta);
253
 
      y7 += x1 * A1[7];
254
 
 
255
 
      y0 += x3 * *A3;   Yget(z4, y[12], beta);
256
 
      y1 += x0 * A0[1];
257
 
      y2 += x1 * A1[2]; Yget(z5, y[13], beta);
258
 
      y3 += x2 * A2[3];
259
 
      y4 += x3 * A3[4]; Yget(z6, y[14], beta);
260
 
      y5 += x0 * A0[5];
261
 
      y6 += x1 * A1[6]; Yget(z7, y[15], beta);
262
 
      y7 += x2 * A2[7];
263
 
 
264
 
      z0 += x0 * A0[8];  *y = y0;
265
 
      z1 += x1 * A1[9];
266
 
      z2 += x2 * A2[10]; y[1] = y1;
267
 
      z3 += x3 * A3[11];
268
 
      z4 += x0 * A0[12]; y[2] = y2;
269
 
      z5 += x1 * A1[13];
270
 
      z6 += x2 * A2[14]; y[3] = y3;
271
 
      z7 += x3 * A3[15];
272
 
 
273
 
      z0 += x1 * A1[8]; y[4] = y4;
274
 
      z1 += x2 * A2[9];
275
 
      z2 += x3 * A3[10]; y[5] = y5;
276
 
      z3 += x0 * A0[11];
277
 
      z4 += x1 * A1[12]; y[6] = y6;
278
 
      z5 += x2 * A2[13];
279
 
      z6 += x3 * A3[14]; y[7] = y7;
280
 
      z7 += x0 * A0[15];
281
 
 
282
 
      z0 += x2 * A2[8];
283
 
      z1 += x3 * A3[ 9];
284
 
      z2 += x0 * A0[10];
285
 
      z3 += x1 * A1[11];
286
 
      z4 += x2 * A2[12];
287
 
      z5 += x3 * A3[13];
288
 
      z6 += x0 * A0[14];
289
 
      z7 += x1 * A1[15];
290
 
 
291
 
      z0 += x3 * A3[8];
292
 
      z1 += x0 * A0[9];
293
 
      z2 += x1 * A1[10];
294
 
      z3 += x2 * A2[11];
295
 
      z4 += x3 * A3[12];
296
 
      z5 += x0 * A0[13];
297
 
      z6 += x1 * A1[14];
298
 
      z7 += x2 * A2[15];
299
 
      y[8] = z0;
300
 
      y[9] = z1;
301
 
      y[10] = z2;
302
 
      y[11] = z3;
303
 
      y[12] = z4;
304
 
      y[13] = z5;
305
 
      y[14] = z6;
306
 
      y[15] = z7;
307
 
      if (M-M16) gemvMlt8(M-M16, N, A0+16, lda, x, beta, y+16);
308
 
   }
309
 
   else if (N) gemvMlt8(M, N, A, lda, x, beta, y);
310
 
}
311
 
 
312
 
static void gemv32x4(const int M, const int N, const TYPE *A, const int lda,
313
 
                     const TYPE *X, const SCALAR beta, TYPE *Y)
314
 
{
315
 
   #ifdef BETA1
316
 
      int j;
317
 
   #endif
318
 
   const int incA = lda<<2;
319
 
 
320
 
   if (N >= 4)
321
 
   {
322
 
      if (M >= 32)
323
 
      {
324
 
         #ifdef BETA1
325
 
            for (j=(N>>2); j; j--, A += incA, X += 4)
326
 
               gemvN32x4(M, 4, A, lda, X, ATL_rone, Y);
327
 
            if ( (j = N-((N>>2)<<2)) ) gemvNle4(M, j, A, lda, X, ATL_rone, Y);
328
 
         #else
329
 
            gemvN32x4(M, 4, A, lda, X, beta, Y);
330
 
            if (N != 4)
331
 
               Mjoin(PATL,gemvN_a1_x1_b1_y1)
332
 
                  (M, N-4, ATL_rone, A+incA, lda, X+4, 1, ATL_rone, Y, 1);
333
 
         #endif
334
 
      }
335
 
      else gemvMlt8(M, N, A, lda, X, beta, Y);
336
 
   }
337
 
   else if (M) gemvNle4(M, N, A, lda, X, beta, Y);
338
 
}
339
 
 
340
 
static void gemvMlt8(const int M, const int N, const TYPE *A, const int lda,
341
 
                     const TYPE *X, const SCALAR beta, TYPE *Y)
342
 
{
343
 
   int i;
344
 
   register TYPE y0;
345
 
   for (i=M; i; i--)
346
 
   {
347
 
      #ifdef BETA0
348
 
         y0 = Mjoin(PATL,dot)(N, A, lda, X, 1);
349
 
      #else
350
 
         Yget(y0, *Y, beta);
351
 
         y0 += Mjoin(PATL,dot)(N, A, lda, X, 1);
352
 
      #endif
353
 
      *Y++ = y0;
354
 
      A++;
355
 
   }
356
 
}
357
 
static void gemvNle4(const int M, const int N, const TYPE *A, const int lda,
358
 
                     const TYPE *X, const SCALAR beta, TYPE *Y)
359
 
{
360
 
   int i;
361
 
   const TYPE *A0 = A, *A1 = A+lda, *A2 = A1+lda, *A3 = A2+lda;
362
 
   register TYPE x0, x1, x2, x3;
363
 
   #ifdef BETAX
364
 
      const register TYPE bet=beta;
365
 
   #endif
366
 
 
367
 
   switch(N)
368
 
   {
369
 
   case 1:
370
 
      #if defined(BETA0)
371
 
         Mjoin(PATL,cpsc)(M, *X, A, 1, Y, 1);
372
 
      #elif defined(BETAX)
373
 
         Mjoin(PATL,axpby)(M, *X, A, 1, beta, Y, 1);
374
 
      #else
375
 
         Mjoin(PATL,axpy)(M, *X, A, 1, Y, 1);
376
 
      #endif
377
 
      break;
378
 
   case 2:
379
 
      x0 = *X; x1 = X[1];
380
 
      for (i=0; i != M; i++)
381
 
      #ifdef BETA0
382
 
         Y[i] = A0[i] * x0 + A1[i] * x1;
383
 
      #elif defined(BETAX)
384
 
         Y[i] = Y[i]*bet + A0[i] * x0 + A1[i] * x1;
385
 
      #else
386
 
         Y[i] += A0[i] * x0 + A1[i] * x1;
387
 
      #endif
388
 
      break;
389
 
   case 3:
390
 
      x0 = *X; x1 = X[1]; x2 = X[2];
391
 
      for (i=0; i != M; i++)
392
 
      #ifdef BETA0
393
 
         Y[i] = A0[i] * x0 + A1[i] * x1 + A2[i] * x2;
394
 
      #elif defined(BETAX)
395
 
         Y[i] = Y[i]*bet + A0[i] * x0 + A1[i] * x1 + A2[i] * x2;
396
 
      #else
397
 
         Y[i] += A0[i] * x0 + A1[i] * x1 + A2[i] * x2;
398
 
      #endif
399
 
      break;
400
 
   case 4:
401
 
      if (M >= 32) gemv32x4(M, 4, A, lda, X, beta, Y);
402
 
      else
403
 
      {
404
 
         x0 = *X; x1 = X[1]; x2 = X[2]; x3 = X[3];
405
 
         for (i=0; i != M; i++)
406
 
         #ifdef BETA0
407
 
            Y[i] = A0[i] * x0 + A1[i] * x1 + A2[i] * x2 + A3[i] * x3;
408
 
         #elif defined(BETAX)
409
 
            Y[i] = Y[i]*bet + A0[i] * x0 + A1[i] * x1 + A2[i] * x2 + A3[i] * x3;
410
 
         #else
411
 
            Y[i] += A0[i] * x0 + A1[i] * x1 + A2[i] * x2 + A3[i] * x3;
412
 
         #endif
413
 
      }
414
 
      break;
415
 
   default:
416
 
      ATL_assert(!N);
417
 
   }
418
 
}
419
 
 
420
 
static void gemv16x2(const int M, const int N, const TYPE *A, const int lda,
421
 
                     const TYPE *X, const SCALAR beta, TYPE *Y)
422
 
/*
423
 
 * 16x2 with feeble prefetch
424
 
 */
425
 
{
426
 
   const int M16 = (M>>4)<<4, N2 = (N>>1)<<1, nr = N-N2;
427
 
   const int incA = lda << 1, incAm = 16 - N2*lda;
428
 
   const TYPE *stX = X + N2 - 2, *x;
429
 
   const TYPE *A0 = A, *A1 = A + lda;
430
 
   TYPE *stY = Y + M16;
431
 
   register TYPE x0, x1;
432
 
   register TYPE y0, y1, y2, y3, y4, y5, y6, y7;
433
 
   register TYPE y8, y9, y10, y11, y12, y13, y14, y15;
434
 
   register TYPE p0, p1;
435
 
 
436
 
   if (N > 4)
437
 
   {
438
 
      if (M16)
439
 
      {
440
 
         do
441
 
         {
442
 
            #ifdef BETA0
443
 
               y0 = y1 = y2 = y3 = y4 = y5 = y6 = y7 =
444
 
               y8 = y9 = y10 = y11 = y12 = y13 = y14 = y15 = ATL_rzero;
445
 
            #elif defined BETAX
446
 
               x0 = beta;
447
 
               y0 = *Y; y1 = Y[1]; y2 = Y[2];  y3 = Y[3];
448
 
               y4 = Y[4]; y5 = Y[5]; y6 = Y[6];  y7 = Y[7];
449
 
               y8 = Y[8]; y9 = Y[9]; y10 = Y[10];  y11 = Y[11];
450
 
               y12 = Y[12]; y13 = Y[13]; y14 = Y[14];  y15 = Y[15];
451
 
               y0 *= x0; y1 *= x0; y2 *= x0; y3 *= x0;
452
 
               y4 *= x0; y5 *= x0; y6 *= x0; y7 *= x0;
453
 
               y8 *= x0; y9 *= x0; y10 *= x0; y11 *= x0;
454
 
               y12 *= x0; y13 *= x0; y14 *= x0; y15 *= x0;
455
 
            #else
456
 
               y0 = *Y; y1 = Y[1]; y2 = Y[2];  y3 = Y[3];
457
 
               y4 = Y[4]; y5 = Y[5]; y6 = Y[6];  y7 = Y[7];
458
 
               y8 = Y[8]; y9 = Y[9]; y10 = Y[10];  y11 = Y[11];
459
 
               y12 = Y[12]; y13 = Y[13]; y14 = Y[14];  y15 = Y[15];
460
 
            #endif
461
 
            p0 = *A0;
462
 
            p1 = A1[1];
463
 
            x = X;
464
 
            do
465
 
            {
466
 
               x0 = *x; x1 = x[1]; x += 2;
467
 
               y0 += x0 * p0; p0 = A0[incA];
468
 
               y1 += x1 * p1; p1 = A1[incA+1];
469
 
               y8 += x0 * A0[8];
470
 
               y9 += x1 * A1[9];
471
 
               y4 += x0 * A0[4];
472
 
               y5 += x1 * A1[5];
473
 
               y12 += x0 * A0[12];
474
 
               y13 += x1 * A1[13];
475
 
               y2 += x0 * A0[2];
476
 
               y3 += x1 * A1[3];
477
 
               y10 += x0 * A0[10];
478
 
               y11 += x1 * A1[11];
479
 
               y6 += x0 * A0[6];
480
 
               y7 += x1 * A1[7];
481
 
               y14 += x0 * A0[14];
482
 
               y15 += x1 * A1[15];
483
 
 
484
 
               y0  += x1 * *A1;
485
 
               y1  += x0 * A0[1];
486
 
               y8  += x1 * A1[8];
487
 
               y9  += x0 * A0[9];
488
 
               y4  += x1 * A1[4];
489
 
               y5  += x0 * A0[5];
490
 
               y12 += x1 * A1[12];
491
 
               y13 += x0 * A0[13];
492
 
               y2  += x1 * A1[2];
493
 
               y3  += x0 * A0[3];
494
 
               y10 += x1 * A1[10];
495
 
               y11 += x0 * A0[11];
496
 
               y6  += x1 * A1[6];
497
 
               y7  += x0 * A0[7];
498
 
               y14 += x1 * A1[14]; A1 += incA;
499
 
               y15 += x0 * A0[15]; A0 += incA;
500
 
            }
501
 
            while (x != stX);
502
 
            if (!nr) /* 2 cols left */
503
 
            {
504
 
               x0 = *x; x1 = x[1]; x += 2;
505
 
               y0 += x0 * p0;
506
 
               y1 += x1 * p1;
507
 
               y8 += x0 * A0[8];
508
 
               y9 += x1 * A1[9];
509
 
               y4 += x0 * A0[4];
510
 
               y5 += x1 * A1[5];
511
 
               y12 += x0 * A0[12];
512
 
               y13 += x1 * A1[13];
513
 
               y2 += x0 * A0[2];
514
 
               y3 += x1 * A1[3];
515
 
               y10 += x0 * A0[10];
516
 
               y11 += x1 * A1[11];
517
 
               y6 += x0 * A0[6];
518
 
               y7 += x1 * A1[7];
519
 
               y14 += x0 * A0[14];
520
 
               y15 += x1 * A1[15];
521
 
 
522
 
               y0  += x1 * *A1;
523
 
               y1  += x0 * A0[1];
524
 
               y8  += x1 * A1[8];
525
 
               y9  += x0 * A0[9];
526
 
               y4  += x1 * A1[4];
527
 
               y5  += x0 * A0[5];
528
 
               y12 += x1 * A1[12];
529
 
               y13 += x0 * A0[13];
530
 
               y2  += x1 * A1[2];
531
 
               y3  += x0 * A0[3];
532
 
               y10 += x1 * A1[10];
533
 
               y11 += x0 * A0[11];
534
 
               y6  += x1 * A1[6];
535
 
               y7  += x0 * A0[7];
536
 
               y14 += x1 * A1[14]; A1 += incA;
537
 
               y15 += x0 * A0[15]; A0 += incA;
538
 
            }
539
 
            else     /* 3 cols left */
540
 
            {
541
 
               x0 = *x; x1 = x[1]; x += 2;
542
 
               y0 += x0 * p0; p0 = A0[incA];
543
 
               y1 += x1 * p1;
544
 
               y8 += x0 * A0[8];
545
 
               y9 += x1 * A1[9];
546
 
               y4 += x0 * A0[4];
547
 
               y5 += x1 * A1[5];
548
 
               y12 += x0 * A0[12];
549
 
               y13 += x1 * A1[13];
550
 
               y2 += x0 * A0[2];
551
 
               y3 += x1 * A1[3];
552
 
               y10 += x0 * A0[10];
553
 
               y11 += x1 * A1[11];
554
 
               y6 += x0 * A0[6];
555
 
               y7 += x1 * A1[7];
556
 
               y14 += x0 * A0[14];
557
 
               y15 += x1 * A1[15];
558
 
 
559
 
               y0  += x1 * *A1;
560
 
               y1  += x0 * A0[1];
561
 
               y8  += x1 * A1[8];
562
 
               y9  += x0 * A0[9];
563
 
               y4  += x1 * A1[4];
564
 
               y5  += x0 * A0[5];
565
 
               y12 += x1 * A1[12];
566
 
               y13 += x0 * A0[13];
567
 
               y2  += x1 * A1[2];
568
 
               y3  += x0 * A0[3];
569
 
               y10 += x1 * A1[10];
570
 
               y11 += x0 * A0[11];
571
 
               y6  += x1 * A1[6];
572
 
               y7  += x0 * A0[7];
573
 
               y14 += x1 * A1[14]; A1 += incA;
574
 
               y15 += x0 * A0[15]; A0 += incA;
575
 
 
576
 
               x0 = *x;
577
 
               y0  += x0 * *A0;
578
 
               y1  += x0 * A0[1];
579
 
               y2  += x0 * A0[2];
580
 
               y3  += x0 * A0[3];
581
 
               y4  += x0 * A0[4];
582
 
               y5  += x0 * A0[5];
583
 
               y6  += x0 * A0[6];
584
 
               y7  += x0 * A0[7];
585
 
               y8  += x0 * A0[8];
586
 
               y9  += x0 * A0[9];
587
 
               y10 += x0 * A0[10];
588
 
               y11 += x0 * A0[11];
589
 
               y12 += x0 * A0[12];
590
 
               y13 += x0 * A0[13];
591
 
               y14 += x0 * A0[14];
592
 
               y15 += x0 * A0[15];
593
 
            }
594
 
            *Y = y0;
595
 
            A0 += incAm;
596
 
            Y[1] = y1;
597
 
            A1 += incAm;
598
 
            Y[2] = y2;
599
 
            Y[3] = y3;
600
 
            Y[4] = y4;
601
 
            Y[5] = y5;
602
 
            Y[6] = y6;
603
 
            Y[7] = y7;
604
 
            Y[8] = y8;
605
 
            Y[9] = y9;
606
 
            Y[10] = y10;
607
 
            Y[11] = y11;
608
 
            Y[12] = y12;
609
 
            Y[13] = y13;
610
 
            Y[14] = y14;
611
 
            Y[15] = y15;
612
 
            Y += 16;
613
 
         }
614
 
         while (Y != stY);
615
 
      }
616
 
      if (M-M16) gemvMlt8(M-M16, N, A0, lda, X, beta, Y);
617
 
   }
618
 
   else if (M) gemvNle4(M, N, A, lda, X, beta, Y);
619
 
}
620
 
void Mjoin(Mjoin(Mjoin(Mjoin(Mjoin(PATL,gemvN),NM),_x1),BNM),_y1)
621
 
   (const int M, const int N, const SCALAR alpha, const TYPE *A, const int lda,
622
 
    const TYPE *X, const int incX, const SCALAR beta, TYPE *Y, const int incY)
623
 
{
624
 
   gemv16x2(M, N, A, lda, X, beta, Y);
625
 
}