← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/blas/level2/ATL_hemv.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-06-11 15:58:16 UTC
  • mfrom: (1.1.3 upstream)
  • mto: (2.2.21 experimental)
  • mto: This revision was merged to the branch mainline in revision 26.
  • Revision ID: package-import@ubuntu.com-20130611155816-b72z8f621tuhbzn0
Tags: upstream-3.10.1
Import upstream version 3.10.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/blas/level2/ATL_hemv.c
1
 
/* ---------------------------------------------------------------------
2
 
 *
3
 
 * -- Automatically Tuned Linear Algebra Software (ATLAS)
4
 
 *    (C) Copyright 2000 All Rights Reserved
5
 
 *
6
 
 * -- ATLAS routine -- Version 3.2 -- December 25, 2000
7
 
 *
8
 
 * Author         : Antoine P. Petitet
9
 
 * Contributor(s) : R. Clint Whaley
10
 
 * Originally developed at the University of Tennessee,
11
 
 * Innovative Computing Laboratory, Knoxville TN, 37996-1301, USA.
12
 
 *
13
 
 * ---------------------------------------------------------------------
14
 
 *
15
 
 * -- Copyright notice and Licensing terms:
16
 
 *
17
 
 *  Redistribution  and  use in  source and binary forms, with or without
18
 
 *  modification, are  permitted provided  that the following  conditions
19
 
 *  are met:
20
 
 *
21
 
 * 1. Redistributions  of  source  code  must retain the above copyright
22
 
 *    notice, this list of conditions and the following disclaimer.
23
 
 * 2. Redistributions in binary form must reproduce  the above copyright
24
 
 *    notice,  this list of conditions, and the  following disclaimer in
25
 
 *    the documentation and/or other materials provided with the distri-
26
 
 *    bution.
27
 
 * 3. The name of the University,  the ATLAS group,  or the names of its
28
 
 *    contributors  may not be used to endorse or promote products deri-
29
 
 *    ved from this software without specific written permission.
30
 
 *
31
 
 * -- Disclaimer:
32
 
 *
33
 
 * THIS  SOFTWARE  IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34
 
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,  INCLUDING,  BUT NOT
35
 
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36
 
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
37
 
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,  INDIRECT, INCIDENTAL, SPE-
38
 
 * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
39
 
 * TO,  PROCUREMENT  OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
40
 
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
41
 
 * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  (IN-
42
 
 * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
43
 
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
44
 
 *
45
 
 * ---------------------------------------------------------------------
46
 
 */
47
1
/*
48
 
 * Include files
 
2
 *             Automatically Tuned Linear Algebra Software v3.10.1
 
3
 * Copyright (C) 2010, 2012 R. Clint Whaley
 
4
 *
 
5
 * Code contributers : R. Clint Whaley, Antoine P. Petitet
 
6
 *
 
7
 * Redistribution and use in source and binary forms, with or without
 
8
 * modification, are permitted provided that the following conditions
 
9
 * are met:
 
10
 *   1. Redistributions of source code must retain the above copyright
 
11
 *      notice, this list of conditions and the following disclaimer.
 
12
 *   2. Redistributions in binary form must reproduce the above copyright
 
13
 *      notice, this list of conditions, and the following disclaimer in the
 
14
 *      documentation and/or other materials provided with the distribution.
 
15
 *   3. The name of the ATLAS group or the names of its contributers may
 
16
 *      not be used to endorse or promote products derived from this
 
17
 *      software without specific written permission.
 
18
 *
 
19
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 
20
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 
21
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 
22
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 
23
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 
24
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 
25
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 
26
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 
27
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 
28
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 
29
 * POSSIBILITY OF SUCH DAMAGE.
 
30
 *
49
31
 */
50
32
#include "atlas_misc.h"
51
33
#include "atlas_level1.h"
52
34
#include "atlas_kernel2.h"
53
35
#include "atlas_lvl2.h"
54
 
#include "atlas_mv.h"
55
 
#include "atlas_r1.h"
56
36
 
57
37
#include "atlas_reflvl2.h"          /* temporary for building purposes */
58
38
#include "atlas_reflevel2.h"        /* used for gbmv, gpmv and gpr.    */
59
 
 
60
 
void Mjoin( PATL, hemv )
61
 
(
62
 
   const enum ATLAS_UPLO      UPLO,
63
 
   const int                  N,
64
 
   const SCALAR               ALPHA,
65
 
   const TYPE                 * A,
66
 
   const int                  LDA,
67
 
   const TYPE                 * X,
68
 
   const int                  INCX,
69
 
   const SCALAR               BETA,
70
 
   TYPE                       * Y,
71
 
   const int                  INCY
72
 
)
73
 
{
 
39
#include Mstr(Mjoin(Mjoin(atlas_,PRE),sysinfo.h))
 
40
#include "atlas_cacheedge.h"
 
41
 
 
42
#ifdef CacheEdge
 
43
   #if !defined(CacheEdge) || CacheEdge == 0
 
44
      #define MY_CE (4*ATL_MulBySize(ATL_L1elts))
 
45
   #else
 
46
      #define MY_CE CacheEdge
 
47
   #endif
 
48
#else
 
49
   #define MY_CE (4*ATL_MulBySize(ATL_L1elts))
 
50
#endif
 
51
 
 
52
/*
 
53
 * Guess as to whether GEMV-based SYRK will beat ref.  Will on any arch
 
54
 * where the compiler is not too whoopy (ARM) or where vectorization provides
 
55
 * faster memory load as well as computation (x86 with vectorization).
 
56
 * Single precision complex has enough flops/load to be almost always
 
57
 * compute-bound enough that the auto-GEMV-tuning outweighs the double load.
 
58
 */
 
59
#if defined(ATL_GAS_ARM)
 
60
   #define USE_GEMV_BASED 1
 
61
#elif defined(SREAL)
 
62
   #ifdef ATL_SSE1
 
63
      #define USE_GEMV_BASED 1
 
64
   #endif
 
65
#elif defined(DREAL)
 
66
   #ifdef ATL_SSE2
 
67
      #define USE_GEMV_BASED 1
 
68
   #endif
 
69
#elif defined(SCPLX)
 
70
   #define USE_GEMV_BASED 1
 
71
#elif defined(DCPLX)
 
72
   #ifdef ATL_SSE2
 
73
      #define USE_GEMV_BASED 1
 
74
   #endif
 
75
#endif
 
76
 
 
77
#ifdef USE_GEMV_BASED
 
78
typedef void (*ATL_symvK_t)
 
79
   (const enum ATLAS_UPLO, const int, const SCALAR, const TYPE*, const int,
 
80
    const TYPE*, const int, const SCALAR, TYPE*, const int);
 
81
 
 
82
static void ATL_symvL
 
83
(
 
84
   ATL_symvK_t symvK,
 
85
   const int NB,
 
86
   ATL_CINT N,
 
87
   const TYPE *A,
 
88
   ATL_CINT lda,
 
89
   const TYPE *x,
 
90
   TYPE  *y,
 
91
   const TYPE *xt,
 
92
   TYPE  *yt
 
93
)
 
94
{
 
95
   const TYPE one[2] = {ATL_rone, ATL_rzero};
 
96
   ATL_INT Mmb2;
 
97
   ATL_INT Mmb, mr, MB, j;
 
98
   const size_t incA = (NB SHIFT)*lda;
 
99
   const size_t opsize = ((size_t)(N+8)*(N+4))*(sizeof(TYPE)>>1)SHIFT;
 
100
   void (*gemvT)(ATL_CINT, ATL_CINT, const SCALAR, const TYPE*, ATL_CINT,
 
101
                 const TYPE*, ATL_CINT, const SCALAR, TYPE*, ATL_CINT);
 
102
   void (*gemvN)(ATL_CINT, ATL_CINT, const SCALAR, const TYPE*, ATL_CINT,
 
103
                 const TYPE*, ATL_CINT, const SCALAR, TYPE*, ATL_CINT);
 
104
 
 
105
   if (opsize > MY_CE)
 
106
   {
 
107
      gemvT = Mjoin(PATL,gemvT);
 
108
      gemvN = Mjoin(PATL,gemvN_L2);
 
109
   }
 
110
   else if (opsize <= ATL_MulBySize(ATL_L1elts))
 
111
   {
 
112
      gemvT = Mjoin(PATL,gemvT_L1);
 
113
      gemvN = Mjoin(PATL,gemvN_L1);
 
114
   }
 
115
   else
 
116
   {
 
117
      gemvT = Mjoin(PATL,gemvT_L2);
 
118
      gemvN = Mjoin(PATL,gemvN_L2);
 
119
   }
 
120
/*
 
121
 * Choose MB such that A is retained in L2 cache for second GEMV call
 
122
 * If partial block is tiny, absorbe it into last block since cache is not
 
123
 * precise anyway.
 
124
 */
 
125
   MB = ATL_DivBySize(MY_CE) / NB;
 
126
   MB = (MB > N || MB < 240) ? N : MB;
 
127
   Mmb2 = Mmb+Mmb;
 
128
   for (j=0; j < N; j += NB, A += incA)
 
129
   {
 
130
      const register size_t j2 = j+j;
 
131
      register int i, nb=N-j;
 
132
      nb = (nb >= NB) ? NB : nb;
 
133
      symvK(AtlasLower, nb, one, A+j2, lda, x+j2, 1, one, y+j2, 1);
 
134
      for (i=j+nb; i < N; i += MB)
 
135
      {
 
136
         const register size_t i2 = i+i;
 
137
         register int mb = N-i;
 
138
         mb = (mb >= MB) ? MB : mb;
 
139
         gemvT(mb, nb, one, A+i2, lda, xt+i2, 1, one, yt+j2, 1);
 
140
         gemvN(mb, nb, one, A+i2, lda, x+j2, 1, one, y+i2, 1);
 
141
      }
 
142
   }
 
143
}
 
144
static void ATL_symvU
 
145
(
 
146
   ATL_symvK_t symvK,
 
147
   const int NB,
 
148
   ATL_CINT N,
 
149
   const TYPE *A,
 
150
   ATL_CINT lda,
 
151
   const TYPE *x,
 
152
   TYPE  *y,
 
153
   const TYPE *xt,
 
154
   TYPE  *yt
 
155
)
 
156
{
 
157
   const TYPE one[2] = {ATL_rone, ATL_rzero};
 
158
   ATL_INT Mmb2;
 
159
   ATL_INT Mmb, mr, MB, j;
 
160
   const size_t incA = (NB SHIFT)*lda;
 
161
   const size_t opsize = ((size_t)(N+8)*(N+4))*(sizeof(TYPE)>>1)SHIFT;
 
162
   void (*gemvT)(ATL_CINT, ATL_CINT, const SCALAR, const TYPE*, ATL_CINT,
 
163
                 const TYPE*, ATL_CINT, const SCALAR, TYPE*, ATL_CINT);
 
164
   void (*gemvN)(ATL_CINT, ATL_CINT, const SCALAR, const TYPE*, ATL_CINT,
 
165
                 const TYPE*, ATL_CINT, const SCALAR, TYPE*, ATL_CINT);
 
166
 
 
167
   if (opsize > MY_CE)
 
168
   {
 
169
      gemvT = Mjoin(PATL,gemvT);
 
170
      gemvN = Mjoin(PATL,gemvN_L2);
 
171
   }
 
172
   else if (opsize <= ATL_MulBySize(ATL_L1elts))
 
173
   {
 
174
      gemvT = Mjoin(PATL,gemvT_L1);
 
175
      gemvN = Mjoin(PATL,gemvN_L1);
 
176
   }
 
177
   else
 
178
   {
 
179
      gemvT = Mjoin(PATL,gemvT_L2);
 
180
      gemvN = Mjoin(PATL,gemvN_L2);
 
181
   }
 
182
/*
 
183
 * Choose MB such that A is retained in L2 cache for second GEMV call
 
184
 * If partial block is tiny, absorbe it into last block since cache is not
 
185
 * precise anyway.
 
186
 */
 
187
   MB = ATL_DivBySize(MY_CE) / NB;
 
188
   MB = (MB > N || MB < 240) ? N : MB;
 
189
   Mmb2 = Mmb+Mmb;
 
190
   for (j=0; j < N; j += NB, A += incA)
 
191
   {
 
192
      const register size_t j2 = j+j;
 
193
      register int i, nb=N-j;
 
194
      nb = (nb >= NB) ? NB : nb;
 
195
      for (i=0; i < j; i += MB)
 
196
      {
 
197
         const register size_t i2 = i+i;
 
198
         register int mb = j-i;
 
199
         mb = (mb >= MB) ? MB : mb;
 
200
         gemvT(mb, nb, one, A+i2, lda, xt+i2, 1, one, yt+j2, 1);
 
201
         gemvN(mb, nb, one, A+i2, lda, x+j2, 1, one, y+i2, 1);
 
202
      }
 
203
      symvK(AtlasUpper, nb, one, A+j2, lda, x+j2, 1, one, y+j2, 1);
 
204
   }
 
205
}
 
206
#endif
 
207
 
 
208
void Mjoin(PATL,hemv)
 
209
(
 
210
   const enum ATLAS_UPLO Uplo,
 
211
   const int             N,
 
212
   const SCALAR          alpha,
 
213
   const TYPE            *A,
 
214
   const int             lda,
 
215
   const TYPE            *X,
 
216
   const int             incX,
 
217
   const SCALAR          beta,
 
218
   TYPE                  *Y,
 
219
   const int             incY
 
220
)
74
221
/*
75
222
 * Purpose
76
223
 * =======
88
235
 *
89
236
 * ---------------------------------------------------------------------
90
237
 */
91
 
/*
92
 
 * .. Local Variables ..
93
 
 */
94
 
   void                       (*gemv0)( const int, const int, const SCALAR,
95
 
                              const TYPE *, const int, const TYPE *, const int,
96
 
                              const SCALAR, TYPE *, const int );
97
 
   void                       (*gemv1)( const int, const int, const SCALAR,
98
 
                              const TYPE *, const int, const TYPE *, const int,
99
 
                              const SCALAR, TYPE *, const int );
100
 
   void                       (*gemvN)( const int, const int, const SCALAR,
101
 
                              const TYPE *, const int, const TYPE *, const int,
102
 
                              const SCALAR, TYPE *, const int );
103
 
#ifdef TREAL
104
 
   TYPE                       alphaY, beta0;
105
 
#define lda2                  LDA
106
 
#define one                   ATL_rone
107
 
#define zero                  ATL_rzero
108
 
   void                       * vx = NULL, * vy = NULL;
109
 
   TYPE                       * x, * y;
110
 
   TYPE                       * A0, * A1, * x0, * x1, * y00, * y0, * y1;
111
 
#else
112
 
   const int                  lda2 = ( LDA SHIFT );
113
 
   const TYPE                 * alphaY, * beta0;
114
 
   const TYPE                 one [2] = { ATL_rone,  ATL_rzero },
115
 
                              zero[2] = { ATL_rzero, ATL_rzero };
116
 
   void                       * vx = NULL, * vy = NULL;
117
 
   TYPE                       * x, * y;
118
 
   TYPE                       * A0, * A1, * x0, * x1, * y00, * y0, * y1;
 
238
{
 
239
   const int BETA0 = (*beta == ATL_rzero && beta[1] == ATL_rzero);
 
240
   const int BETA1 = (*beta == ATL_rone && beta[1] == ATL_rzero);
 
241
   const int ALPHA1 = (*alpha == ATL_rone && alpha[1] == ATL_rzero);
 
242
   const int ALPHA0 = (*alpha == ATL_rzero && alpha[1] == ATL_rzero);
 
243
   if (N <= 0 || (ALPHA0 && BETA1))
 
244
      return;
 
245
   if (ALPHA0)
 
246
   {
 
247
      if (BETA0)
 
248
         Mjoin(PATL,zero)(N, Y, incY);
 
249
      else
 
250
         Mjoin(PATL,scal)(N, beta, Y, incY);
 
251
      return;
 
252
   }
 
253
#ifdef USE_GEMV_BASED
 
254
   if (N >= 240)
 
255
   {
 
256
      void *vp=NULL;
 
257
      TYPE *x=(TYPE*)X, *y=Y, *xh, *yh;
 
258
      const size_t tX = (size_t)X, tY = (size_t)Y, N2 = N+N;
 
259
      const int COPYY = !(incY == 1 &&
 
260
                          (ATL_MulByCachelen(ATL_DivByCachelen(tY)) == tY));
 
261
      const int COPYX = !(incX == 1 && (COPYY || ALPHA1) &&
 
262
                          (ATL_MulByCachelen(ATL_DivByCachelen(tX)) == tX));
 
263
      const TYPE one[2] = {ATL_rone, ATL_rzero};
 
264
      const TYPE *calp=one, *cbet=one;
 
265
      TYPE *tp;
 
266
      tp = vp = malloc((COPYX+COPYY+2)*(ATL_Cachelen+ATL_MulBySize(N)));
 
267
      if (!vp)
 
268
      {
 
269
         Mjoin(PATL,refhemv)(Uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
 
270
         return;
 
271
      }
 
272
      yh = ATL_AlignPtr(tp);
 
273
      Mjoin(PATL,zero)(N, yh, 1);
 
274
      tp = yh + N2;
 
275
      xh = ATL_AlignPtr(tp);
 
276
      tp = xh + N2;
 
277
      if (COPYX)
 
278
      {
 
279
         x = ATL_AlignPtr(tp);
 
280
         if (COPYY || ALPHA1)
 
281
         {
 
282
            register ATL_INT i;
 
283
            register const size_t incX2 = incX+incX;
 
284
            const TYPE *xx=X;
 
285
            for (i=0; i < N2; i += 2, xx += incX2)
 
286
            {
 
287
               xh[i] = x[i] = *xx;
 
288
               xh[i+1] = -(x[i+1] = xx[1]);
 
289
            }
 
290
         }
 
291
         else if (alpha[1] == ATL_rzero)
 
292
         {
 
293
            register ATL_INT i;
 
294
            register const size_t incX2 = incX+incX;
 
295
            register const TYPE ra=(*alpha), ia=alpha[1];
 
296
            const TYPE *xx=X;
 
297
            for (i=0; i < N2; i += 2, xx += incX2)
 
298
            {
 
299
               register TYPE rx = *xx, ix = xx[1];
 
300
               x[i] = rx*ra;
 
301
               x[i+1] = ix*ra;
 
302
               xh[i] = rx;
 
303
               xh[i+1] = -ix;
 
304
            }
 
305
         }
 
306
         else
 
307
         {
 
308
            register ATL_INT i;
 
309
            register const size_t incX2 = incX+incX;
 
310
            register const TYPE ra=(*alpha), ia=alpha[1];
 
311
            const TYPE *xx=X;
 
312
            for (i=0; i < N2; i += 2, xx += incX2)
 
313
            {
 
314
               register TYPE rx = *xx, ix = xx[1];
 
315
               x[i] = rx*ra - ix*ia;
 
316
               x[i+1] = rx*ia + ix*ra;
 
317
               xh[i] = rx;
 
318
               xh[i+1] = -ix;
 
319
            }
 
320
         }
 
321
         tp = x + N2;
 
322
      }
 
323
      else
 
324
         Mjoin(PATL,copyConj)(N, X, incX, xh, 1);
 
325
      if (COPYY)
 
326
      {
 
327
         calp = alpha;
 
328
         cbet = beta;
 
329
         y = ATL_AlignPtr(tp);
 
330
         Mjoin(PATL,zero)(N, y, 1);
 
331
      }
 
332
      else if (BETA0)
 
333
         Mjoin(PATL,zero)(N, y, 1);
 
334
      else if (!BETA1)
 
335
         Mjoin(PATL,scal)(N, beta, y, 1);
 
336
      if (Uplo == AtlasLower)
 
337
         ATL_symvL(Mjoin(PATL,refhemv), 120, N, A, lda, x, y, xh, yh);
 
338
      else
 
339
         ATL_symvU(Mjoin(PATL,refhemv), 120, N, A, lda, x, y, xh, yh);
 
340
      if (COPYY)
 
341
      {
 
342
         Mjoin(PATL,axpbyConj)(N, alpha, yh, 1, calp, y, 1);
 
343
         Mjoin(PATL,axpby)(N, one, y, 1, cbet, Y, incY);
 
344
      }
 
345
      else
 
346
         Mjoin(PATL,axpyConj)(N, alpha, yh, 1, Y, incY);
 
347
      free(vp);
 
348
      return;
 
349
   }
119
350
#endif
120
 
   int                        incA, incA1, incXY, incXY1, j, jb, mb, mb1, n, nb;
121
 
/* ..
122
 
 * .. Executable Statements ..
123
 
 *
124
 
 */
125
 
   if( N == 0 ) return;
126
 
 
127
 
   if( SCALAR_IS_ZERO( ALPHA ) )
128
 
   {
129
 
      if( !( SCALAR_IS_ONE( BETA ) ) ) Mjoin( PATL, scal )( N, BETA, Y, INCY );
130
 
      return;
131
 
   }
132
 
 
133
 
   if( ( INCX != 1 ) || ( ( INCY == 1 ) && !( SCALAR_IS_ONE( ALPHA ) ) ) )
134
 
   {
135
 
      vx = (void *)malloc( ATL_Cachelen + ATL_MulBySize( N ) );
136
 
      ATL_assert( vx ); x = ATL_AlignPtr( vx );
137
 
      Mjoin( PATL, cpsc )( N, ALPHA, X, INCX, x, 1 );
138
 
      alphaY = one;
139
 
   }
140
 
   else { x = (TYPE *)(X); alphaY = ALPHA; }
141
 
 
142
 
   if( ( INCY != 1 ) || !( SCALAR_IS_ONE( alphaY ) ) )
143
 
   {
144
 
      vy = malloc( ATL_Cachelen + ATL_MulBySize( N ) );
145
 
      ATL_assert( vy ); y00 = y = ATL_AlignPtr( vy );
146
 
      beta0 = zero;
147
 
   }
148
 
   else { y00 = y = (TYPE *)(Y); beta0 = BETA; }
149
 
 
150
 
   if(      SCALAR_IS_ZERO( beta0 ) ) gemv0 = Mjoin( PATL, gemvC_a1_x1_b0_y1 );
151
 
   else if( SCALAR_IS_ONE ( beta0 ) ) gemv0 = Mjoin( PATL, gemvC_a1_x1_b1_y1 );
152
 
   else                               gemv0 = Mjoin( PATL, gemvC_a1_x1_bX_y1 );
153
 
   gemv1 = Mjoin( PATL, gemvC_a1_x1_b1_y1 );
154
 
   gemvN = Mjoin( PATL, gemvS_a1_x1_b1_y1 );
155
 
 
156
 
   ATL_GetPartSYMV( A, LDA, &mb, &nb );
157
 
 
158
 
   mb1   = N - ( ( N - 1 ) / mb ) * mb;
159
 
   incA1 = nb * lda2;  incXY1 = (nb SHIFT);
160
 
 
161
 
   if( UPLO == AtlasUpper )
162
 
   {
163
 
      incA  = ( incXY = (mb SHIFT) ) + mb * lda2;
164
 
      A0 = (TYPE *)(A) + mb * lda2; x0 = x + incXY; y0 = y + incXY;
165
 
 
166
 
      for( n  = N - mb; n > 0; n -= mb, A0 += incA, A += incA, x0 += incXY,
167
 
           x += incXY, y0 += incXY, y += incXY )
168
 
      {
169
 
         Mjoin( PATL, hemvU )( mb, A, LDA, x, beta0, y );
170
 
 
171
 
         for( j  =  0, A1 = A0,     x1  = x0,     y1  = y0;     j < n;
172
 
              j += nb, A1 += incA1, x1 += incXY1, y1 += incXY1 )
173
 
         {
174
 
            jb = n - j; jb = Mmin( jb, nb );
175
 
            gemv0( jb, mb, one, A1, LDA, x,  1, beta0, y1, 1 );
176
 
            gemvN( mb, jb, one, A1, LDA, x1, 1, one,   y,  1 );
177
 
         }
178
 
         beta0 = one; gemv0 = gemv1;
179
 
      }
180
 
      Mjoin( PATL, hemvU )( mb1, A, LDA, x, beta0, y );
181
 
   }
182
 
   else
183
 
   {
184
 
      incA = incXY = (mb SHIFT);
185
 
      A0 = (TYPE *)(A); x0 = x; y0 = y;
186
 
 
187
 
      for( n  = N - mb, A += ((N-mb) SHIFT), x += ((N-mb) SHIFT),
188
 
           y += ((N-mb) SHIFT); n > 0; n -= mb, A -= incA, x -= incXY,
189
 
           y -= incXY )
190
 
      {
191
 
         Mjoin( PATL, hemvL )( mb, A+n*lda2, LDA, x, beta0, y );
192
 
 
193
 
         for( j  =  0, A1  = (TYPE *)(A), x1  = x0,     y1  = y0;      j < n;
194
 
              j += nb, A1 += incA1,       x1 += incXY1, y1 += incXY1 )
195
 
         {
196
 
            jb = n - j; jb = Mmin( jb, nb );
197
 
            gemv0( jb, mb, one, A1, LDA, x,  1, beta0, y1, 1 );
198
 
            gemvN( mb, jb, one, A1, LDA, x1, 1, one,   y,  1 );
199
 
         }
200
 
         beta0 = one; gemv0 = gemv1;
201
 
      }
202
 
      Mjoin( PATL, hemvL )( mb1, A0, LDA, x0, beta0, y0 );
203
 
   }
204
 
 
205
 
   if( vx ) free( vx );
206
 
   if( vy )
207
 
   { Mjoin( PATL, axpby )( N, alphaY, y00, 1, BETA, Y, INCY ); free( vy ); }
208
 
/*
209
 
 * End of Mjoin( PATL, hemv )
210
 
 */
 
351
   Mjoin(PATL,refhemv)(Uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
211
352
}