← 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_4x2_0.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2002-04-13 10:07:52 UTC
  • Revision ID: james.westby@ubuntu.com-20020413100752-va9zm0rd4gpurdkq
Tags: upstream-3.2.1ln
Import upstream version 3.2.1ln

Show diffs side-by-side

added added

removed removed

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