kernel/zarch/dgemv_n_4.c

   1 /***************************************************************************
   2 Copyright (c) 2017, The OpenBLAS Project
   3 All rights reserved.
   4 Redistribution and use in source and binary forms, with or without
   5 modification, are permitted provided that the following conditions are
   6 met:
   7 1. Redistributions of source code must retain the above copyright
   8 notice, this list of conditions and the following disclaimer.
   9 2. Redistributions in binary form must reproduce the above copyright
  10 notice, this list of conditions and the following disclaimer in
  11 the documentation and/or other materials provided with the
  12 distribution.
  13 3. Neither the name of the OpenBLAS project nor the names of
  14 its contributors may be used to endorse or promote products
  15 derived from this software without specific prior written permission.
  16 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  17 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  18 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  19 ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
  20 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  21 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  22 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  23 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  24 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
  25 USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26 *****************************************************************************/
  27
  28
  29 #include "common.h"
  30
  31 #define NBMAX 2048
  32
  33 #define HAVE_KERNEL_4x4_VEC 1
  34 #define HAVE_KERNEL_4x2_VEC 1
  35 #define HAVE_KERNEL_4x1_VEC 1
  36
  37 #if defined(HAVE_KERNEL_4x4_VEC) || defined(HAVE_KERNEL_4x2_VEC) || defined(HAVE_KERNEL_4x1_VEC)
  38  #include <vecintrin.h>
  39 #endif
  40
  41 #ifdef HAVE_KERNEL_4x4
  42
  43 #elif HAVE_KERNEL_4x4_VEC
  44
  45 static void dgemv_kernel_4x4(BLASLONG n, FLOAT **ap, FLOAT *xo, FLOAT *y, FLOAT *alpha)
  46 {
  47     BLASLONG i;
  48     FLOAT x0,x1,x2,x3;
  49     x0 = xo[0] * *alpha;
  50     x1 = xo[1] * *alpha;
  51     x2 = xo[2] * *alpha;
  52     x3 = xo[3] * *alpha;
  53     __vector double   v_x0 = {x0,x0};
  54     __vector double   v_x1 = {x1,x1};
  55     __vector double   v_x2 = {x2,x2};
  56     __vector double   v_x3 = {x3,x3};
  57     __vector double* v_y =(__vector double*)y;
  58     __vector double* va0 = (__vector double*)ap[0];
  59     __vector double* va1 = (__vector double*)ap[1];
  60     __vector double* va2 = (__vector double*)ap[2];
  61     __vector double* va3 = (__vector double*)ap[3];
  62
  63     for ( i=0; i< n/2; i+=2 )
  64     {
  65         v_y[i]   += v_x0 * va0[i]   +  v_x1 * va1[i]   + v_x2 * va2[i]   + v_x3 * va3[i] ;
  66         v_y[i+1] += v_x0 * va0[i+1] +  v_x1 * va1[i+1] + v_x2 * va2[i+1] + v_x3 * va3[i+1] ;
  67     }
  68 }
  69
  70 #else
  71
  72 static void dgemv_kernel_4x4(BLASLONG n, FLOAT **ap, FLOAT *xo, FLOAT *y, FLOAT *alpha)
  73 {
  74     BLASLONG i;
  75     FLOAT *a0,*a1,*a2,*a3;
  76     FLOAT x[4]  __attribute__ ((aligned (16)));
  77     a0 = ap[0];
  78     a1 = ap[1];
  79     a2 = ap[2];
  80     a3 = ap[3];
  81
  82     for ( i=0; i<4; i++)
  83         x[i] = xo[i] * *alpha;
  84
  85     for ( i=0; i< n; i+=4 )
  86     {
  87         y[i] += a0[i]*x[0] + a1[i]*x[1] + a2[i]*x[2] + a3[i]*x[3];
  88         y[i+1] += a0[i+1]*x[0] + a1[i+1]*x[1] + a2[i+1]*x[2] + a3[i+1]*x[3];
  89         y[i+2] += a0[i+2]*x[0] + a1[i+2]*x[1] + a2[i+2]*x[2] + a3[i+2]*x[3];
  90         y[i+3] += a0[i+3]*x[0] + a1[i+3]*x[1] + a2[i+3]*x[2] + a3[i+3]*x[3];
  91     }
  92 }
  93
  94
  95 #endif
  96
  97 #ifdef HAVE_KERNEL_4x2
  98
  99 #elif HAVE_KERNEL_4x2_VEC
 100
 101 static void dgemv_kernel_4x2(BLASLONG n, FLOAT **ap, FLOAT *xo, FLOAT *y, FLOAT *alpha)
 102 {
 103     BLASLONG i;
 104     FLOAT x0,x1;
 105     x0 = xo[0] * *alpha;
 106     x1 = xo[1] * *alpha;
 107     __vector double   v_x0 = {x0,x0};
 108     __vector double   v_x1 = {x1,x1};
 109     __vector double* v_y =(__vector double*)y;
 110     __vector double* va0 = (__vector double*)ap[0];
 111     __vector double* va1 = (__vector double*)ap[1];
 112
 113     for ( i=0; i< n/2; i+=2 )
 114     {
 115         v_y[i]   += v_x0 * va0[i] +  v_x1 * va1[i]   ;
 116         v_y[i+1] += v_x0 * va0[i+1] +  v_x1 * va1[i+1]  ;
 117     }
 118 }
 119 #else
 120
 121 static void dgemv_kernel_4x2(BLASLONG n, FLOAT **ap, FLOAT *xo, FLOAT *y, FLOAT *alpha)
 122 {
 123     BLASLONG i;
 124     FLOAT *a0,*a1;
 125     FLOAT x[4]  __attribute__ ((aligned (16)));
 126     a0 = ap[0];
 127     a1 = ap[1];
 128
 129     for ( i=0; i<2; i++)
 130         x[i] = xo[i] * *alpha;
 131
 132     for ( i=0; i< n; i+=4 )
 133     {
 134         y[i] += a0[i]*x[0] + a1[i]*x[1];
 135         y[i+1] += a0[i+1]*x[0] + a1[i+1]*x[1];
 136         y[i+2] += a0[i+2]*x[0] + a1[i+2]*x[1];
 137         y[i+3] += a0[i+3]*x[0] + a1[i+3]*x[1];
 138     }
 139 }
 140
 141
 142 #endif
 143
 144 #ifdef HAVE_KERNEL_4x1
 145
 146 #elif HAVE_KERNEL_4x1_VEC
 147 static void dgemv_kernel_4x1(BLASLONG n, FLOAT *ap, FLOAT *xo, FLOAT *y, FLOAT *alpha)
 148 {
 149
 150     BLASLONG i;
 151     FLOAT x0;
 152     x0 = xo[0] * *alpha;
 153     __vector double   v_x0 = {x0,x0};
 154     __vector double* v_y =(__vector double*)y;
 155     __vector double* va0 = (__vector double*)ap;
 156
 157     for ( i=0; i< n/2; i+=2 )
 158     {
 159         v_y[i] += v_x0 * va0[i]    ;
 160         v_y[i+1] += v_x0 * va0[i+1]  ;
 161     }
 162
 163
 164 }
 165
 166 #else
 167 static void dgemv_kernel_4x1(BLASLONG n, FLOAT *ap, FLOAT *xo, FLOAT *y, FLOAT *alpha)
 168 {
 169     BLASLONG i;
 170     FLOAT *a0;
 171     FLOAT x[4]  __attribute__ ((aligned (16)));
 172     a0 = ap;
 173
 174     for ( i=0; i<1; i++)
 175         x[i] = xo[i] * *alpha;
 176
 177     for ( i=0; i< n; i+=4 )
 178     {
 179         y[i] += a0[i]*x[0];
 180         y[i+1] += a0[i+1]*x[0];
 181         y[i+2] += a0[i+2]*x[0];
 182         y[i+3] += a0[i+3]*x[0];
 183     }
 184 }
 185
 186
 187 #endif
 188
 189
 190
 191 static void add_y(BLASLONG n, FLOAT *src, FLOAT *dest, BLASLONG inc_dest)
 192 {
 193     BLASLONG i;
 194
 195     for ( i=0; i<n; i++ ){
 196             *dest += *src;
 197             src++;
 198             dest += inc_dest;
 199     }
 200     return;
 201
 202 }
 203
 204 int CNAME(BLASLONG m, BLASLONG n, BLASLONG dummy1, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y, FLOAT *buffer)
 205 {
 206     BLASLONG i;
 207     BLASLONG j;
 208     FLOAT *a_ptr;
 209     FLOAT *x_ptr;
 210     FLOAT *y_ptr;
 211     FLOAT *ap[4];
 212     BLASLONG n1;
 213     BLASLONG m1;
 214     BLASLONG m2;
 215     BLASLONG m3;
 216     BLASLONG n2;
 217     BLASLONG lda4 =  lda << 2;
 218     FLOAT xbuffer[8],*ybuffer;
 219
 220     if ( m < 1 ) return(0);
 221     if ( n < 1 ) return(0);
 222
 223     ybuffer = buffer;
 224
 225     n1 = n >> 2 ;
 226     n2 = n &  3 ;
 227
 228     m3 = m & 3  ;
 229     m1 = m & -4 ;
 230     m2 = (m & (NBMAX-1)) - m3 ;
 231
 232     y_ptr = y;
 233
 234     BLASLONG NB = NBMAX;
 235
 236     while ( NB == NBMAX )
 237     {
 238
 239         m1 -= NB;
 240         if ( m1 < 0)
 241         {
 242             if ( m2 == 0 ) break;
 243             NB = m2;
 244         }
 245
 246         a_ptr = a;
 247         x_ptr = x;
 248
 249         ap[0] = a_ptr;
 250         ap[1] = a_ptr + lda;
 251         ap[2] = ap[1] + lda;
 252         ap[3] = ap[2] + lda;
 253
 254         if ( inc_y != 1 )
 255             memset(ybuffer,0,NB*8);
 256         else
 257             ybuffer = y_ptr;
 258
 259         if ( inc_x == 1 )
 260         {
 261
 262
 263             for( i = 0; i < n1 ; i++)
 264             {
 265                 dgemv_kernel_4x4(NB,ap,x_ptr,ybuffer,&alpha);
 266                 ap[0] += lda4;
 267                 ap[1] += lda4;
 268                 ap[2] += lda4;
 269                 ap[3] += lda4;
 270                 a_ptr += lda4;
 271                 x_ptr += 4;
 272             }
 273
 274             if ( n2 & 2 )
 275             {
 276                 dgemv_kernel_4x2(NB,ap,x_ptr,ybuffer,&alpha);
 277                 a_ptr += lda*2;
 278                 x_ptr += 2;
 279             }
 280
 281
 282             if ( n2 & 1 )
 283             {
 284                 dgemv_kernel_4x1(NB,a_ptr,x_ptr,ybuffer,&alpha);
 285                 a_ptr += lda;
 286                 x_ptr += 1;
 287
 288             }
 289
 290
 291         }
 292         else
 293         {
 294
 295             for( i = 0; i < n1 ; i++)
 296             {
 297                 xbuffer[0] = x_ptr[0];
 298                 x_ptr += inc_x;
 299                 xbuffer[1] =  x_ptr[0];
 300                 x_ptr += inc_x;
 301                 xbuffer[2] =  x_ptr[0];
 302                 x_ptr += inc_x;
 303                 xbuffer[3] = x_ptr[0];
 304                 x_ptr += inc_x;
 305                 dgemv_kernel_4x4(NB,ap,xbuffer,ybuffer,&alpha);
 306                 ap[0] += lda4;
 307                 ap[1] += lda4;
 308                 ap[2] += lda4;
 309                 ap[3] += lda4;
 310                 a_ptr += lda4;
 311             }
 312
 313             for( i = 0; i < n2 ; i++)
 314             {
 315                 xbuffer[0] = x_ptr[0];
 316                 x_ptr += inc_x;
 317                 dgemv_kernel_4x1(NB,a_ptr,xbuffer,ybuffer,&alpha);
 318                 a_ptr += lda;
 319
 320             }
 321
 322         }
 323
 324         a     += NB;
 325         if ( inc_y != 1 )
 326         {
 327             add_y(NB,ybuffer,y_ptr,inc_y);
 328             y_ptr += NB * inc_y;
 329         }
 330         else
 331             y_ptr += NB ;
 332
 333     }
 334
 335     if ( m3 == 0 ) return(0);
 336
 337     if ( m3 == 3 )
 338     {
 339         a_ptr = a;
 340         x_ptr = x;
 341         FLOAT temp0 = 0.0;
 342         FLOAT temp1 = 0.0;
 343         FLOAT temp2 = 0.0;
 344         if ( lda == 3 && inc_x ==1 )
 345         {
 346
 347             for( i = 0; i < ( n & -4 ); i+=4 )
 348             {
 349
 350                 temp0 += a_ptr[0] * x_ptr[0] + a_ptr[3] * x_ptr[1];
 351                 temp1 += a_ptr[1] * x_ptr[0] + a_ptr[4] * x_ptr[1];
 352                 temp2 += a_ptr[2] * x_ptr[0] + a_ptr[5] * x_ptr[1];
 353
 354                 temp0 += a_ptr[6] * x_ptr[2] + a_ptr[9]  * x_ptr[3];
 355                 temp1 += a_ptr[7] * x_ptr[2] + a_ptr[10] * x_ptr[3];
 356                 temp2 += a_ptr[8] * x_ptr[2] + a_ptr[11] * x_ptr[3];
 357
 358                 a_ptr += 12;
 359                 x_ptr += 4;
 360             }
 361
 362             for( ; i < n; i++ )
 363             {
 364                 temp0 += a_ptr[0] * x_ptr[0];
 365                 temp1 += a_ptr[1] * x_ptr[0];
 366                 temp2 += a_ptr[2] * x_ptr[0];
 367                 a_ptr += 3;
 368                 x_ptr ++;
 369             }
 370
 371         }
 372         else
 373         {
 374
 375             for( i = 0; i < n; i++ )
 376             {
 377                 temp0 += a_ptr[0] * x_ptr[0];
 378                 temp1 += a_ptr[1] * x_ptr[0];
 379                 temp2 += a_ptr[2] * x_ptr[0];
 380                 a_ptr += lda;
 381                 x_ptr += inc_x;
 382
 383
 384             }
 385
 386         }
 387         y_ptr[0] += alpha * temp0;
 388         y_ptr += inc_y;
 389         y_ptr[0] += alpha * temp1;
 390         y_ptr += inc_y;
 391         y_ptr[0] += alpha * temp2;
 392         return(0);
 393     }
 394
 395
 396     if ( m3 == 2 )
 397     {
 398         a_ptr = a;
 399         x_ptr = x;
 400         FLOAT temp0 = 0.0;
 401         FLOAT temp1 = 0.0;
 402         if ( lda == 2 && inc_x ==1 )
 403         {
 404
 405             for( i = 0; i < (n & -4) ; i+=4 )
 406             {
 407                 temp0 += a_ptr[0] * x_ptr[0] + a_ptr[2] * x_ptr[1];
 408                 temp1 += a_ptr[1] * x_ptr[0] + a_ptr[3] * x_ptr[1];
 409                 temp0 += a_ptr[4] * x_ptr[2] + a_ptr[6] * x_ptr[3];
 410                 temp1 += a_ptr[5] * x_ptr[2] + a_ptr[7] * x_ptr[3];
 411                 a_ptr += 8;
 412                 x_ptr += 4;
 413
 414             }
 415
 416
 417             for( ; i < n; i++ )
 418             {
 419                 temp0 += a_ptr[0]   * x_ptr[0];
 420                 temp1 += a_ptr[1]   * x_ptr[0];
 421                 a_ptr += 2;
 422                 x_ptr ++;
 423             }
 424
 425         }
 426         else
 427         {
 428
 429             for( i = 0; i < n; i++ )
 430             {
 431                 temp0 += a_ptr[0] * x_ptr[0];
 432                 temp1 += a_ptr[1] * x_ptr[0];
 433                 a_ptr += lda;
 434                 x_ptr += inc_x;
 435
 436
 437             }
 438
 439         }
 440         y_ptr[0] += alpha * temp0;
 441         y_ptr += inc_y;
 442         y_ptr[0] += alpha * temp1;
 443         return(0);
 444     }
 445
 446     if ( m3 == 1 )
 447     {
 448         a_ptr = a;
 449         x_ptr = x;
 450         FLOAT temp = 0.0;
 451         if ( lda == 1 && inc_x ==1 )
 452         {
 453
 454             for( i = 0; i < (n & -4); i+=4 )
 455             {
 456                 temp += a_ptr[i] * x_ptr[i] + a_ptr[i+1] * x_ptr[i+1] + a_ptr[i+2] * x_ptr[i+2] + a_ptr[i+3] * x_ptr[i+3];
 457
 458             }
 459
 460             for( ; i < n; i++ )
 461             {
 462                 temp += a_ptr[i] * x_ptr[i];
 463             }
 464
 465         }
 466         else
 467         {
 468
 469             for( i = 0; i < n; i++ )
 470             {
 471                 temp += a_ptr[0] * x_ptr[0];
 472                 a_ptr += lda;
 473                 x_ptr += inc_x;
 474             }
 475
 476         }
 477         y_ptr[0] += alpha * temp;
 478         return(0);
 479     }
 480
 481
 482     return(0);
 483 }
 484
 485