← 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, 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:
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
}