kernel/mips/zdot_msa.c

   1 /*******************************************************************************
   2 Copyright (c) 2016, 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 #include "common.h"
  29 #include "macros_msa.h"
  30
  31 #if !defined(CONJ)
  32     #define OP1     -=
  33     #define OP2     +=
  34     #define OP3     -
  35     #define OP4     +
  36 #else
  37     #define OP1     +=
  38     #define OP2     -=
  39     #define OP3     +
  40     #define OP4     -
  41 #endif
  42
  43 OPENBLAS_COMPLEX_FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y)
  44 {
  45     BLASLONG i = 0;
  46     FLOAT dot[2];
  47     BLASLONG inc_x2, inc_y2;
  48     v2f64 vx0, vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8, vx9, vx10, vx11;
  49     v2f64 vy0, vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8, vy9, vy10, vy11;
  50     v2f64 vx0r, vx0i, vx1r, vx1i, vx2r, vx2i, vx3r, vx3i;
  51     v2f64 vy0r, vy0i, vy1r, vy1i, vy2r, vy2i, vy3r, vy3i;
  52     v2f64 dot0 = {0, 0};
  53     v2f64 dot1 = {0, 0};
  54     v2f64 dot2 = {0, 0};
  55     v2f64 dot3 = {0, 0};
  56     v2f64 dot4 = {0, 0};
  57     v2f64 dot5 = {0, 0};
  58     v2f64 dot6 = {0, 0};
  59     v2f64 dot7 = {0, 0};
  60     v2f64 zero = {0, 0};
  61     OPENBLAS_COMPLEX_FLOAT result;
  62
  63     dot[0] = 0.0;
  64     dot[1] = 0.0;
  65
  66     CREAL(result) = 0.0;
  67     CIMAG(result) = 0.0;
  68
  69     if (n < 1) return (result);
  70
  71     inc_x2 = 2 * inc_x;
  72     inc_y2 = 2 * inc_y;
  73
  74     if ((1 == inc_x) && (1 == inc_y))
  75     {
  76         if (n > 7)
  77         {
  78             FLOAT *x_pref, *y_pref;
  79             BLASLONG pref_offset;
  80
  81             pref_offset = (BLASLONG)x & (L1_DATA_LINESIZE - 1);
  82             if (pref_offset > 0)
  83             {
  84                 pref_offset = L1_DATA_LINESIZE - pref_offset;
  85                 pref_offset = pref_offset / sizeof(FLOAT);
  86             }
  87             x_pref = x + pref_offset + 32 + 8;
  88
  89             pref_offset = (BLASLONG)y & (L1_DATA_LINESIZE - 1);
  90             if (pref_offset > 0)
  91             {
  92                 pref_offset = L1_DATA_LINESIZE - pref_offset;
  93                 pref_offset = pref_offset / sizeof(FLOAT);
  94             }
  95             y_pref = y + pref_offset + 32 + 8;
  96
  97             LD_DP4_INC(x, 2, vx0, vx1, vx2, vx3);
  98             LD_DP4_INC(y, 2, vy0, vy1, vy2, vy3);
  99
 100             PCKEVOD_D2_DP(vx1, vx0, vx0r, vx0i);
 101             PCKEVOD_D2_DP(vy1, vy0, vy0r, vy0i);
 102
 103             for (i = (n >> 3) - 1; i--;)
 104             {
 105                 PREF_OFFSET(x_pref, 0);
 106                 PREF_OFFSET(x_pref, 32);
 107                 PREF_OFFSET(x_pref, 64);
 108                 PREF_OFFSET(x_pref, 96);
 109                 PREF_OFFSET(y_pref, 0);
 110                 PREF_OFFSET(y_pref, 32);
 111                 PREF_OFFSET(y_pref, 64);
 112                 PREF_OFFSET(y_pref, 96);
 113                 x_pref += 16;
 114                 y_pref += 16;
 115
 116                 vx4 = LD_DP(x); x += 2;
 117                 vx1r = (v2f64) __msa_pckev_d((v2i64) vx3, (v2i64) vx2);
 118                 dot0 += (vx0r * vy0r);
 119                 vx5 = LD_DP(x); x += 2;
 120                 vx1i = (v2f64) __msa_pckod_d((v2i64) vx3, (v2i64) vx2);
 121                 dot1 OP2 (vx0i * vy0r);
 122                 vy4 = LD_DP(y); y += 2;
 123                 vy1r = (v2f64) __msa_pckev_d((v2i64) vy3, (v2i64) vy2);
 124                 dot2 += (vx1r * vy1r);
 125                 vy5 = LD_DP(y); y += 2;
 126                 vy1i = (v2f64) __msa_pckod_d((v2i64) vy3, (v2i64) vy2);
 127                 dot3 OP2 (vx1i * vy1r);
 128                 vx6 = LD_DP(x); x += 2;
 129                 vx7 = LD_DP(x); x += 2;
 130                 vy6 = LD_DP(y); y += 2;
 131                 vy7 = LD_DP(y); y += 2;
 132                 vx8 = LD_DP(x); x += 2;
 133                 dot0 OP1 (vx0i * vy0i);
 134                 vx9 = LD_DP(x); x += 2;
 135                 vx2r = (v2f64) __msa_pckev_d((v2i64) vx5, (v2i64) vx4);
 136                 dot1 += (vx0r * vy0i);
 137                 vy8 = LD_DP(y); y += 2;
 138                 vx2i = (v2f64) __msa_pckod_d((v2i64) vx5, (v2i64) vx4);
 139                 dot2 OP1 (vx1i * vy1i);
 140                 vy9 = LD_DP(y); y += 2;
 141                 vy2r = (v2f64) __msa_pckev_d((v2i64) vy5, (v2i64) vy4);
 142                 dot3 += (vx1r * vy1i);
 143                 vx10 = LD_DP(x); x += 2;
 144                 vy2i = (v2f64) __msa_pckod_d((v2i64) vy5, (v2i64) vy4);
 145                 vx11 = LD_DP(x); x += 2;
 146                 vx3r = (v2f64) __msa_pckev_d((v2i64) vx7, (v2i64) vx6);
 147                 dot4 += (vx2r * vy2r);
 148                 vy10 = LD_DP(y); y += 2;
 149                 vx3i = (v2f64) __msa_pckod_d((v2i64) vx7, (v2i64) vx6);
 150                 dot5 OP2 (vx2i * vy2r);
 151                 vy11 = LD_DP(y); y += 2;
 152                 vy3r = (v2f64) __msa_pckev_d((v2i64) vy7, (v2i64) vy6);
 153                 vy3i = (v2f64) __msa_pckod_d((v2i64) vy7, (v2i64) vy6);
 154                 dot6 += (vx3r * vy3r);
 155                 vx0r = (v2f64) __msa_pckev_d((v2i64) vx9, (v2i64) vx8);
 156                 dot7 OP2 (vx3i * vy3r);
 157                 vx0i = (v2f64) __msa_pckod_d((v2i64) vx9, (v2i64) vx8);
 158                 vy0r = (v2f64) __msa_pckev_d((v2i64) vy9, (v2i64) vy8);
 159                 vx2 = vx10;
 160                 vy0i = (v2f64) __msa_pckod_d((v2i64) vy9, (v2i64) vy8);
 161                 vx3 = vx11;
 162                 dot4 OP1 (vx2i * vy2i);
 163                 vy2 = vy10;
 164                 dot5 += (vx2r * vy2i);
 165                 vy3 = vy11;
 166                 dot6 OP1 (vx3i * vy3i);
 167                 dot7 += (vx3r * vy3i);
 168             }
 169
 170             vx4 = LD_DP(x); x += 2;
 171             vx1r = (v2f64) __msa_pckev_d((v2i64) vx3, (v2i64) vx2);
 172             dot0 += (vx0r * vy0r);
 173             vx5 = LD_DP(x); x += 2;
 174             vx1i = (v2f64) __msa_pckod_d((v2i64) vx3, (v2i64) vx2);
 175             dot1 OP2 (vx0i * vy0r);
 176             vy4 = LD_DP(y); y += 2;
 177             vy1r = (v2f64) __msa_pckev_d((v2i64) vy3, (v2i64) vy2);
 178             dot2 += (vx1r * vy1r);
 179             vy5 = LD_DP(y); y += 2;
 180             vy1i = (v2f64) __msa_pckod_d((v2i64) vy3, (v2i64) vy2);
 181             dot3 OP2 (vx1i * vy1r);
 182             vx6 = LD_DP(x); x += 2;
 183             vx7 = LD_DP(x); x += 2;
 184             vy6 = LD_DP(y); y += 2;
 185             vy7 = LD_DP(y); y += 2;
 186             dot0 OP1 (vx0i * vy0i);
 187             vx2r = (v2f64) __msa_pckev_d((v2i64) vx5, (v2i64) vx4);
 188             dot1 += (vx0r * vy0i);
 189             vx2i = (v2f64) __msa_pckod_d((v2i64) vx5, (v2i64) vx4);
 190             dot2 OP1 (vx1i * vy1i);
 191             vy2r = (v2f64) __msa_pckev_d((v2i64) vy5, (v2i64) vy4);
 192             dot3 += (vx1r * vy1i);
 193             vy2i = (v2f64) __msa_pckod_d((v2i64) vy5, (v2i64) vy4);
 194             vx3r = (v2f64) __msa_pckev_d((v2i64) vx7, (v2i64) vx6);
 195             dot4 += (vx2r * vy2r);
 196             vx3i = (v2f64) __msa_pckod_d((v2i64) vx7, (v2i64) vx6);
 197             dot5 OP2 (vx2i * vy2r);
 198             vy3r = (v2f64) __msa_pckev_d((v2i64) vy7, (v2i64) vy6);
 199             vy3i = (v2f64) __msa_pckod_d((v2i64) vy7, (v2i64) vy6);
 200             dot6 += (vx3r * vy3r);
 201             dot7 OP2 (vx3i * vy3r);
 202             dot4 OP1 (vx2i * vy2i);
 203             dot5 += (vx2r * vy2i);
 204             dot6 OP1 (vx3i * vy3i);
 205             dot7 += (vx3r * vy3i);
 206         }
 207     }
 208     else if (n > 7)
 209     {
 210         LD_DP4_INC(x, inc_x2, vx0, vx1, vx2, vx3);
 211         LD_DP4_INC(y, inc_y2, vy0, vy1, vy2, vy3);
 212
 213         PCKEVOD_D2_DP(vx1, vx0, vx0r, vx0i);
 214         PCKEVOD_D2_DP(vy1, vy0, vy0r, vy0i);
 215
 216         for (i = (n >> 3) - 1; i--;)
 217         {
 218             vx4 = LD_DP(x); x += inc_x2;
 219             vx1r = (v2f64) __msa_pckev_d((v2i64) vx3, (v2i64) vx2);
 220             dot0 += (vx0r * vy0r);
 221             vx5 = LD_DP(x); x += inc_x2;
 222             vx1i = (v2f64) __msa_pckod_d((v2i64) vx3, (v2i64) vx2);
 223             dot1 OP2 (vx0i * vy0r);
 224             vy4 = LD_DP(y); y += inc_y2;
 225             vy1r = (v2f64) __msa_pckev_d((v2i64) vy3, (v2i64) vy2);
 226             dot2 += (vx1r * vy1r);
 227             vy5 = LD_DP(y); y += inc_y2;
 228             vy1i = (v2f64) __msa_pckod_d((v2i64) vy3, (v2i64) vy2);
 229             dot3 OP2 (vx1i * vy1r);
 230             vx6 = LD_DP(x); x += inc_x2;
 231             vx7 = LD_DP(x); x += inc_x2;
 232             vy6 = LD_DP(y); y += inc_y2;
 233             vy7 = LD_DP(y); y += inc_y2;
 234             vx8 = LD_DP(x); x += inc_x2;
 235             dot0 OP1 (vx0i * vy0i);
 236             vx9 = LD_DP(x); x += inc_x2;
 237             vx2r = (v2f64) __msa_pckev_d((v2i64) vx5, (v2i64) vx4);
 238             dot1 += (vx0r * vy0i);
 239             vy8 = LD_DP(y); y += inc_y2;
 240             vx2i = (v2f64) __msa_pckod_d((v2i64) vx5, (v2i64) vx4);
 241             dot2 OP1 (vx1i * vy1i);
 242             vy9 = LD_DP(y); y += inc_y2;
 243             vy2r = (v2f64) __msa_pckev_d((v2i64) vy5, (v2i64) vy4);
 244             dot3 += (vx1r * vy1i);
 245             vx10 = LD_DP(x); x += inc_x2;
 246             vy2i = (v2f64) __msa_pckod_d((v2i64) vy5, (v2i64) vy4);
 247             vx11 = LD_DP(x); x += inc_x2;
 248             vx3r = (v2f64) __msa_pckev_d((v2i64) vx7, (v2i64) vx6);
 249             dot4 += (vx2r * vy2r);
 250             vy10 = LD_DP(y); y += inc_y2;
 251             vx3i = (v2f64) __msa_pckod_d((v2i64) vx7, (v2i64) vx6);
 252             dot5 OP2 (vx2i * vy2r);
 253             vy11 = LD_DP(y); y += inc_y2;
 254             vy3r = (v2f64) __msa_pckev_d((v2i64) vy7, (v2i64) vy6);
 255             vy3i = (v2f64) __msa_pckod_d((v2i64) vy7, (v2i64) vy6);
 256             dot6 += (vx3r * vy3r);
 257             vx0r = (v2f64) __msa_pckev_d((v2i64) vx9, (v2i64) vx8);
 258             dot7 OP2 (vx3i * vy3r);
 259             vx0i = (v2f64) __msa_pckod_d((v2i64) vx9, (v2i64) vx8);
 260             vy0r = (v2f64) __msa_pckev_d((v2i64) vy9, (v2i64) vy8);
 261             vx2 = vx10;
 262             vy0i = (v2f64) __msa_pckod_d((v2i64) vy9, (v2i64) vy8);
 263             vx3 = vx11;
 264             dot4 OP1 (vx2i * vy2i);
 265             vy2 = vy10;
 266             dot5 += (vx2r * vy2i);
 267             vy3 = vy11;
 268             dot6 OP1 (vx3i * vy3i);
 269             dot7 += (vx3r * vy3i);
 270         }
 271
 272         vx4 = LD_DP(x); x += inc_x2;
 273         vx1r = (v2f64) __msa_pckev_d((v2i64) vx3, (v2i64) vx2);
 274         dot0 += (vx0r * vy0r);
 275         vx5 = LD_DP(x); x += inc_x2;
 276         vx1i = (v2f64) __msa_pckod_d((v2i64) vx3, (v2i64) vx2);
 277         dot1 OP2 (vx0i * vy0r);
 278         vy4 = LD_DP(y); y += inc_y2;
 279         vy1r = (v2f64) __msa_pckev_d((v2i64) vy3, (v2i64) vy2);
 280         dot2 += (vx1r * vy1r);
 281         vy5 = LD_DP(y); y += inc_y2;
 282         vy1i = (v2f64) __msa_pckod_d((v2i64) vy3, (v2i64) vy2);
 283         dot3 OP2 (vx1i * vy1r);
 284         vx6 = LD_DP(x); x += inc_x2;
 285         vx7 = LD_DP(x); x += inc_x2;
 286         vy6 = LD_DP(y); y += inc_y2;
 287         vy7 = LD_DP(y); y += inc_y2;
 288         dot0 OP1 (vx0i * vy0i);
 289         vx2r = (v2f64) __msa_pckev_d((v2i64) vx5, (v2i64) vx4);
 290         dot1 += (vx0r * vy0i);
 291         vx2i = (v2f64) __msa_pckod_d((v2i64) vx5, (v2i64) vx4);
 292         dot2 OP1 (vx1i * vy1i);
 293         vy2r = (v2f64) __msa_pckev_d((v2i64) vy5, (v2i64) vy4);
 294         dot3 += (vx1r * vy1i);
 295         vy2i = (v2f64) __msa_pckod_d((v2i64) vy5, (v2i64) vy4);
 296         vx3r = (v2f64) __msa_pckev_d((v2i64) vx7, (v2i64) vx6);
 297         dot4 += (vx2r * vy2r);
 298         vx3i = (v2f64) __msa_pckod_d((v2i64) vx7, (v2i64) vx6);
 299         dot5 OP2 (vx2i * vy2r);
 300         vy3r = (v2f64) __msa_pckev_d((v2i64) vy7, (v2i64) vy6);
 301         vy3i = (v2f64) __msa_pckod_d((v2i64) vy7, (v2i64) vy6);
 302         dot6 += (vx3r * vy3r);
 303         dot7 OP2 (vx3i * vy3r);
 304         dot4 OP1 (vx2i * vy2i);
 305         dot5 += (vx2r * vy2i);
 306         dot6 OP1 (vx3i * vy3i);
 307         dot7 += (vx3r * vy3i);
 308     }
 309
 310     if (n & 7)
 311     {
 312         if (n & 4)
 313         {
 314             LD_DP4_INC(x, inc_x2, vx0, vx1, vx2, vx3);
 315             LD_DP4_INC(y, inc_y2, vy0, vy1, vy2, vy3);
 316
 317             PCKEVOD_D2_DP(vx1, vx0, vx0r, vx0i);
 318             PCKEVOD_D2_DP(vx3, vx2, vx1r, vx1i);
 319
 320             PCKEVOD_D2_DP(vy1, vy0, vy0r, vy0i);
 321             PCKEVOD_D2_DP(vy3, vy2, vy1r, vy1i);
 322
 323             dot0 += (vx0r * vy0r);
 324             dot0 OP1 (vx0i * vy0i);
 325             dot1 OP2 (vx0i * vy0r);
 326             dot1 += (vx0r * vy0i);
 327
 328             dot2 += (vx1r * vy1r);
 329             dot2 OP1 (vx1i * vy1i);
 330             dot3 OP2 (vx1i * vy1r);
 331             dot3 += (vx1r * vy1i);
 332         }
 333
 334         if (n & 2)
 335         {
 336             LD_DP2_INC(x, inc_x2, vx0, vx1);
 337             LD_DP2_INC(y, inc_y2, vy0, vy1);
 338             PCKEVOD_D2_DP(vx1, vx0, vx0r, vx0i);
 339             PCKEVOD_D2_DP(vy1, vy0, vy0r, vy0i);
 340
 341             dot0 += (vx0r * vy0r);
 342             dot0 OP1 (vx0i * vy0i);
 343             dot1 OP2 (vx0i * vy0r);
 344             dot1 += (vx0r * vy0i);
 345         }
 346
 347         if (n & 1)
 348         {
 349             vx0 = LD_DP(x);
 350             vy0 = LD_DP(y);
 351             PCKEVOD_D2_DP(zero, vx0, vx0r, vx0i);
 352             PCKEVOD_D2_DP(zero, vy0, vy0r, vy0i);
 353
 354             dot0 += (vx0r * vy0r);
 355             dot0 OP1 (vx0i * vy0i);
 356             dot1 OP2 (vx0i * vy0r);
 357             dot1 += (vx0r * vy0i);
 358         }
 359     }
 360
 361     dot0 += dot2 + dot4 + dot6;
 362     dot1 += dot3 + dot5 + dot7;
 363
 364     dot[0] += (dot0[0] + dot0[1]);
 365     dot[1] += (dot1[0] + dot1[1]);
 366
 367     CREAL(result) = dot[0];
 368     CIMAG(result) = dot[1];
 369
 370     return (result);
 371 }