kernel/mips/dtrsm_kernel_RT_8x4_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 static __attribute__ ((noinline))
  32 void dsolve_8x4_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
  33 {
  34     v2f64 src_c0, src_c1, src_c2, src_c3, src_c4, src_c5, src_c6, src_c7;
  35     v2f64 src_c8, src_c9, src_c10, src_c11, src_c12, src_c13, src_c14, src_c15;
  36     v2f64 src_b0, src_b4, src_b5, src_b8, src_b9, src_b10, src_b12, src_b13;
  37     v2f64 src_b14, src_b15;
  38     FLOAT *c_nxt1line = c + ldc;
  39     FLOAT *c_nxt2line = c + 2 * ldc;
  40     FLOAT *c_nxt3line = c + 3 * ldc;
  41
  42     LD_DP4(c, 2, src_c0, src_c1, src_c2, src_c3);
  43     LD_DP4(c_nxt1line, 2, src_c4, src_c5, src_c6, src_c7);
  44     LD_DP4(c_nxt2line, 2, src_c8, src_c9, src_c10, src_c11);
  45     LD_DP4(c_nxt3line, 2, src_c12, src_c13, src_c14, src_c15);
  46
  47     if (bk > 0)
  48     {
  49         BLASLONG i, pref_offset;
  50         FLOAT *pba = a, *pbb = b, *pa0_pref;
  51         v2f64 src_b, src_b0, src_b1, src_a0, src_a1, src_a2, src_a3;
  52
  53         pref_offset = (uintptr_t)a & (L1_DATA_LINESIZE - 1);
  54
  55         if (pref_offset)
  56         {
  57             pref_offset = L1_DATA_LINESIZE - pref_offset;
  58             pref_offset = pref_offset / sizeof(FLOAT);
  59         }
  60
  61         pa0_pref = a + pref_offset;
  62
  63         for (i = (bk >> 1); i--;)
  64         {
  65             PREF_OFFSET(pa0_pref, 128);
  66             PREF_OFFSET(pa0_pref, 160);
  67             PREF_OFFSET(pa0_pref, 192);
  68             PREF_OFFSET(pa0_pref, 224);
  69
  70             LD_DP4_INC(pba, 2, src_a0, src_a1, src_a2, src_a3);
  71             LD_DP2_INC(pbb, 2, src_b0, src_b1);
  72
  73             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
  74             src_c0 -= src_a0 * src_b;
  75             src_c1 -= src_a1 * src_b;
  76             src_c2 -= src_a2 * src_b;
  77             src_c3 -= src_a3 * src_b;
  78
  79             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
  80             src_c4 -= src_a0 * src_b;
  81             src_c5 -= src_a1 * src_b;
  82             src_c6 -= src_a2 * src_b;
  83             src_c7 -= src_a3 * src_b;
  84
  85             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
  86             src_c8  -= src_a0 * src_b;
  87             src_c9  -= src_a1 * src_b;
  88             src_c10 -= src_a2 * src_b;
  89             src_c11 -= src_a3 * src_b;
  90
  91             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
  92             src_c12 -= src_a0 * src_b;
  93             src_c13 -= src_a1 * src_b;
  94             src_c14 -= src_a2 * src_b;
  95             src_c15 -= src_a3 * src_b;
  96
  97             LD_DP4_INC(pba, 2, src_a0, src_a1, src_a2, src_a3);
  98             LD_DP2_INC(pbb, 2, src_b0, src_b1);
  99
 100             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 101             src_c0 -= src_a0 * src_b;
 102             src_c1 -= src_a1 * src_b;
 103             src_c2 -= src_a2 * src_b;
 104             src_c3 -= src_a3 * src_b;
 105
 106             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 107             src_c4 -= src_a0 * src_b;
 108             src_c5 -= src_a1 * src_b;
 109             src_c6 -= src_a2 * src_b;
 110             src_c7 -= src_a3 * src_b;
 111
 112             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
 113             src_c8  -= src_a0 * src_b;
 114             src_c9  -= src_a1 * src_b;
 115             src_c10 -= src_a2 * src_b;
 116             src_c11 -= src_a3 * src_b;
 117
 118             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
 119             src_c12 -= src_a0 * src_b;
 120             src_c13 -= src_a1 * src_b;
 121             src_c14 -= src_a2 * src_b;
 122             src_c15 -= src_a3 * src_b;
 123
 124             pa0_pref += 16;
 125         }
 126
 127         if (bk & 1)
 128         {
 129             LD_DP4_INC(pba, 2, src_a0, src_a1, src_a2, src_a3);
 130             LD_DP2_INC(pbb, 2, src_b0, src_b1);
 131
 132             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 133             src_c0 -= src_a0 * src_b;
 134             src_c1 -= src_a1 * src_b;
 135             src_c2 -= src_a2 * src_b;
 136             src_c3 -= src_a3 * src_b;
 137
 138             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 139             src_c4 -= src_a0 * src_b;
 140             src_c5 -= src_a1 * src_b;
 141             src_c6 -= src_a2 * src_b;
 142             src_c7 -= src_a3 * src_b;
 143
 144             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
 145             src_c8  -= src_a0 * src_b;
 146             src_c9  -= src_a1 * src_b;
 147             src_c10 -= src_a2 * src_b;
 148             src_c11 -= src_a3 * src_b;
 149
 150             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
 151             src_c12 -= src_a0 * src_b;
 152             src_c13 -= src_a1 * src_b;
 153             src_c14 -= src_a2 * src_b;
 154             src_c15 -= src_a3 * src_b;
 155         }
 156     }
 157
 158     a -= 32;
 159     b -= 16;
 160
 161     src_b12 = LD_DP(b + 12);
 162     src_b13 = (v2f64) __msa_splati_d((v2i64) src_b12, 1);
 163     src_b12 = (v2f64) __msa_splati_d((v2i64) src_b12, 0);
 164     src_b14 = LD_DP(b + 14);
 165     src_b15 = (v2f64) __msa_splati_d((v2i64) src_b14, 1);
 166     src_b14 = (v2f64) __msa_splati_d((v2i64) src_b14, 0);
 167
 168     src_b8 = LD_DP(b + 8);
 169     src_b9 = (v2f64) __msa_splati_d((v2i64) src_b8, 1);
 170     src_b8 = (v2f64) __msa_splati_d((v2i64) src_b8, 0);
 171     src_b10 = __msa_cast_to_vector_double(*(b + 10));
 172     src_b10 = (v2f64) __msa_splati_d((v2i64) src_b10, 0);
 173
 174     src_b0 = __msa_cast_to_vector_double(*(b + 0));
 175     src_b0 = (v2f64) __msa_splati_d((v2i64) src_b0, 0);
 176     src_b4 = LD_DP(b + 4);
 177     src_b5 = (v2f64) __msa_splati_d((v2i64) src_b4, 1);
 178     src_b4 = (v2f64) __msa_splati_d((v2i64) src_b4, 0);
 179
 180     src_c12 *= src_b15;
 181     src_c13 *= src_b15;
 182     src_c14 *= src_b15;
 183     src_c15 *= src_b15;
 184
 185     src_c8 -= src_c12 * src_b14;
 186     src_c9 -= src_c13 * src_b14;
 187     src_c10 -= src_c14 * src_b14;
 188     src_c11 -= src_c15 * src_b14;
 189
 190     src_c8 *= src_b10;
 191     src_c9 *= src_b10;
 192     src_c10 *= src_b10;
 193     src_c11 *= src_b10;
 194
 195     src_c4 -= src_c12 * src_b13;
 196     src_c5 -= src_c13 * src_b13;
 197     src_c6 -= src_c14 * src_b13;
 198     src_c7 -= src_c15 * src_b13;
 199
 200     src_c4 -= src_c8 * src_b9;
 201     src_c5 -= src_c9 * src_b9;
 202     src_c6 -= src_c10 * src_b9;
 203     src_c7 -= src_c11 * src_b9;
 204
 205     src_c4 *= src_b5;
 206     src_c5 *= src_b5;
 207     src_c6 *= src_b5;
 208     src_c7 *= src_b5;
 209
 210     src_c0 -= src_c12 * src_b12;
 211     src_c1 -= src_c13 * src_b12;
 212     src_c2 -= src_c14 * src_b12;
 213     src_c3 -= src_c15 * src_b12;
 214
 215     src_c0 -= src_c8 * src_b8;
 216     src_c1 -= src_c9 * src_b8;
 217     src_c2 -= src_c10 * src_b8;
 218     src_c3 -= src_c11 * src_b8;
 219
 220     src_c0 -= src_c4 * src_b4;
 221     src_c1 -= src_c5 * src_b4;
 222     src_c2 -= src_c6 * src_b4;
 223     src_c3 -= src_c7 * src_b4;
 224
 225     src_c0 *= src_b0;
 226     src_c1 *= src_b0;
 227     src_c2 *= src_b0;
 228     src_c3 *= src_b0;
 229
 230     ST_DP4(src_c12, src_c13, src_c14, src_c15, c_nxt3line, 2);
 231     ST_DP4(src_c12, src_c13, src_c14, src_c15, a + 24, 2);
 232     ST_DP4(src_c8, src_c9, src_c10, src_c11, c_nxt2line, 2);
 233     ST_DP4(src_c8, src_c9, src_c10, src_c11, a + 16, 2);
 234     ST_DP4(src_c4, src_c5, src_c6, src_c7, c_nxt1line, 2);
 235     ST_DP4(src_c4, src_c5, src_c6, src_c7, a + 8, 2);
 236     ST_DP4(src_c0, src_c1, src_c2, src_c3, c, 2);
 237     ST_DP4(src_c0, src_c1, src_c2, src_c3, a, 2);
 238 }
 239
 240 static void dsolve_8x2_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 241 {
 242     v2f64 src_c0, src_c1, src_c2, src_c3, src_c4, src_c5, src_c6, src_c7;
 243     v2f64 src_b0, src_b2, src_b3;
 244
 245     LD_DP4(c, 2, src_c0, src_c1, src_c2, src_c3);
 246     LD_DP4(c + ldc, 2, src_c4, src_c5, src_c6, src_c7);
 247
 248     if (bk > 0)
 249     {
 250         BLASLONG i;
 251         FLOAT *pba = a, *pbb = b;
 252         v2f64 src_b, src_b1, src_a0, src_a1, src_a2, src_a3;
 253         v2f64 src_a4, src_a5, src_a6, src_a7;
 254
 255         LD_DP4(pba, 2, src_a0, src_a1, src_a2, src_a3);
 256         src_b0 = LD_DP(pbb);
 257
 258         for (i = bk - 1; i--;)
 259         {
 260             pba += 8;
 261             pbb += 2;
 262
 263             LD_DP4(pba, 2, src_a4, src_a5, src_a6, src_a7);
 264             src_b1 = LD_DP(pbb);
 265
 266             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 267             src_c0 -= src_a0 * src_b;
 268             src_c1 -= src_a1 * src_b;
 269             src_c2 -= src_a2 * src_b;
 270             src_c3 -= src_a3 * src_b;
 271
 272             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 273             src_c4 -= src_a0 * src_b;
 274             src_c5 -= src_a1 * src_b;
 275             src_c6 -= src_a2 * src_b;
 276             src_c7 -= src_a3 * src_b;
 277
 278             src_a0 = src_a4;
 279             src_a1 = src_a5;
 280             src_a2 = src_a6;
 281             src_a3 = src_a7;
 282             src_b0 = src_b1;
 283         }
 284
 285         src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 286         src_c0 -= src_a0 * src_b;
 287         src_c1 -= src_a1 * src_b;
 288         src_c2 -= src_a2 * src_b;
 289         src_c3 -= src_a3 * src_b;
 290
 291         src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 292         src_c4 -= src_a0 * src_b;
 293         src_c5 -= src_a1 * src_b;
 294         src_c6 -= src_a2 * src_b;
 295         src_c7 -= src_a3 * src_b;
 296     }
 297
 298     a -= 16;
 299     b -= 4;
 300
 301     src_b0 = __msa_cast_to_vector_double(*(b + 0));
 302     src_b0 = (v2f64) __msa_splati_d((v2i64) src_b0, 0);
 303     src_b2 = LD_DP(b + 2);
 304     src_b3 = (v2f64) __msa_splati_d((v2i64) src_b2, 1);
 305     src_b2 = (v2f64) __msa_splati_d((v2i64) src_b2, 0);
 306
 307     src_c4 *= src_b3;
 308     src_c5 *= src_b3;
 309     src_c6 *= src_b3;
 310     src_c7 *= src_b3;
 311
 312     src_c0 -= src_c4 * src_b2;
 313     src_c1 -= src_c5 * src_b2;
 314     src_c2 -= src_c6 * src_b2;
 315     src_c3 -= src_c7 * src_b2;
 316
 317     src_c0 *= src_b0;
 318     src_c1 *= src_b0;
 319     src_c2 *= src_b0;
 320     src_c3 *= src_b0;
 321
 322     ST_DP4(src_c0, src_c1, src_c2, src_c3, c, 2);
 323     ST_DP4(src_c4, src_c5, src_c6, src_c7, c + ldc, 2);
 324
 325     ST_DP4(src_c0, src_c1, src_c2, src_c3, a, 2);
 326     ST_DP4(src_c4, src_c5, src_c6, src_c7, a + 8, 2);
 327 }
 328
 329 static void dsolve_8x1_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
 330 {
 331     v2f64 src_c0, src_c1, src_c2, src_c3;
 332     v2f64 src_b0;
 333
 334     LD_DP4(c, 2, src_c0, src_c1, src_c2, src_c3);
 335
 336     if (bk > 0)
 337     {
 338         BLASLONG i;
 339         FLOAT *aa = a, *bb = b;
 340         v2f64 src_a0, src_a1, src_a2, src_a3, src_a4, src_a5, src_a6, src_a7;
 341         v2f64 src_b1;
 342
 343         LD_DP4(aa, 2, src_a0, src_a1, src_a2, src_a3);
 344         src_b0 = LD_DP(bb);
 345
 346         aa += 8;
 347         bb += 1;
 348
 349         for (i = (bk - 1); i--;)
 350         {
 351             LD_DP4(aa, 2, src_a4, src_a5, src_a6, src_a7);
 352             src_b1 = LD_DP(bb);
 353
 354             src_b0 = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 355             src_c0 -= src_a0 * src_b0;
 356             src_c1 -= src_a1 * src_b0;
 357             src_c2 -= src_a2 * src_b0;
 358             src_c3 -= src_a3 * src_b0;
 359
 360             src_a0 = src_a4;
 361             src_a1 = src_a5;
 362             src_a2 = src_a6;
 363             src_a3 = src_a7;
 364             src_b0 = src_b1;
 365
 366             aa += 8;
 367             bb += 1;
 368         }
 369
 370         src_b0 = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 371         src_c0 -= src_a0 * src_b0;
 372         src_c1 -= src_a1 * src_b0;
 373         src_c2 -= src_a2 * src_b0;
 374         src_c3 -= src_a3 * src_b0;
 375     }
 376
 377     a -= 8;
 378     b -= 1;
 379
 380     src_b0 = __msa_cast_to_vector_double(*b);
 381     src_b0 = (v2f64) __msa_splati_d((v2i64) src_b0, 0);
 382
 383     src_c0 *= src_b0;
 384     src_c1 *= src_b0;
 385     src_c2 *= src_b0;
 386     src_c3 *= src_b0;
 387
 388     ST_DP4(src_c0, src_c1, src_c2, src_c3, c, 2);
 389     ST_DP4(src_c0, src_c1, src_c2, src_c3, a, 2);
 390 }
 391
 392 static void dsolve_4x4_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 393 {
 394     v2f64 src_c0, src_c1, src_c2, src_c3, src_c4, src_c5, src_c6, src_c7;
 395     v2f64 src_b0, src_b4, src_b5, src_b8, src_b9, src_b10, src_b12, src_b13;
 396     v2f64 src_b14, src_b15;
 397
 398     LD_DP2(c, 2, src_c0, src_c1);
 399     LD_DP2(c + ldc, 2, src_c2, src_c3);
 400     LD_DP2(c + 2 * ldc, 2, src_c4, src_c5);
 401     LD_DP2(c + 3 * ldc, 2, src_c6, src_c7);
 402
 403     if (bk > 0)
 404     {
 405         BLASLONG i;
 406         FLOAT *aa = a, *bb = b;
 407         v2f64 src_a0, src_a1, src_b, src_b0, src_b1;
 408
 409         for (i = bk; i--;)
 410         {
 411             LD_DP2(aa, 2, src_a0, src_a1);
 412             LD_DP2(bb, 2, src_b0, src_b1);
 413
 414             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 415             src_c0 -= src_a0 * src_b;
 416             src_c1 -= src_a1 * src_b;
 417
 418             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 419             src_c2 -= src_a0 * src_b;
 420             src_c3 -= src_a1 * src_b;
 421
 422             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
 423             src_c4 -= src_a0 * src_b;
 424             src_c5 -= src_a1 * src_b;
 425
 426             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
 427             src_c6 -= src_a0 * src_b;
 428             src_c7 -= src_a1 * src_b;
 429
 430             aa += 4;
 431             bb += 4;
 432         }
 433     }
 434
 435     a -= 16;
 436     b -= 16;
 437
 438     src_b12 = LD_DP(b + 12);
 439     src_b13 = (v2f64) __msa_splati_d((v2i64) src_b12, 1);
 440     src_b12 = (v2f64) __msa_splati_d((v2i64) src_b12, 0);
 441     src_b14 = LD_DP(b + 14);
 442     src_b15 = (v2f64) __msa_splati_d((v2i64) src_b14, 1);
 443     src_b14 = (v2f64) __msa_splati_d((v2i64) src_b14, 0);
 444
 445     src_b8 = LD_DP(b + 8);
 446     src_b9 = (v2f64) __msa_splati_d((v2i64) src_b8, 1);
 447     src_b8 = (v2f64) __msa_splati_d((v2i64) src_b8, 0);
 448     src_b10 = __msa_cast_to_vector_double(*(b + 10));
 449     src_b10 = (v2f64) __msa_splati_d((v2i64) src_b10, 0);
 450
 451     src_b0 = __msa_cast_to_vector_double(*(b + 0));
 452     src_b0 = (v2f64) __msa_splati_d((v2i64) src_b0, 0);
 453     src_b4 = LD_DP(b + 4);
 454     src_b5 = (v2f64) __msa_splati_d((v2i64) src_b4, 1);
 455     src_b4 = (v2f64) __msa_splati_d((v2i64) src_b4, 0);
 456
 457     src_c6 *= src_b15;
 458     src_c7 *= src_b15;
 459
 460     src_c4 -= src_c6 * src_b14;
 461     src_c5 -= src_c7 * src_b14;
 462
 463     src_c4 *= src_b10;
 464     src_c5 *= src_b10;
 465
 466     src_c2 -= src_c6 * src_b13;
 467     src_c3 -= src_c7 * src_b13;
 468
 469     src_c2 -= src_c4 * src_b9;
 470     src_c3 -= src_c5 * src_b9;
 471
 472     src_c2 *= src_b5;
 473     src_c3 *= src_b5;
 474
 475     src_c0 -= src_c6 * src_b12;
 476     src_c1 -= src_c7 * src_b12;
 477
 478     src_c0 -= src_c4 * src_b8;
 479     src_c1 -= src_c5 * src_b8;
 480
 481     src_c0 -= src_c2 * src_b4;
 482     src_c1 -= src_c3 * src_b4;
 483
 484     src_c0 *= src_b0;
 485     src_c1 *= src_b0;
 486
 487     ST_DP2(src_c6, src_c7, c + 3 * ldc, 2);
 488     ST_DP2(src_c4, src_c5, c + 2 * ldc, 2);
 489     ST_DP2(src_c2, src_c3, c + ldc, 2);
 490     ST_DP2(src_c0, src_c1, c, 2);
 491
 492     ST_DP4(src_c4, src_c5, src_c6, src_c7, a + 8, 2);
 493     ST_DP4(src_c0, src_c1, src_c2, src_c3, a, 2);
 494 }
 495
 496 static void dsolve_4x2_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 497 {
 498     v2f64 src_c0, src_c1, src_c2, src_c3, src_b0, src_b2, src_b3;
 499
 500     LD_DP2(c, 2, src_c0, src_c1);
 501     LD_DP2(c + ldc, 2, src_c2, src_c3);
 502
 503     if (bk > 0)
 504     {
 505         BLASLONG i;
 506         FLOAT *aa = a, *bb = b;
 507         v2f64 src_a0, src_a1, src_b, src_b0;
 508
 509         for (i = bk; i--;)
 510         {
 511             LD_DP2(aa, 2, src_a0, src_a1);
 512             src_b0 = LD_DP(bb);
 513
 514             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 515             src_c0 -= src_a0 * src_b;
 516             src_c1 -= src_a1 * src_b;
 517
 518             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 519             src_c2 -= src_a0 * src_b;
 520             src_c3 -= src_a1 * src_b;
 521
 522             aa += 4;
 523             bb += 2;
 524         }
 525     }
 526
 527     a -= 8;
 528     b -= 4;
 529
 530     src_b0 = __msa_cast_to_vector_double(*(b + 0));
 531     src_b0 = (v2f64) __msa_splati_d((v2i64) src_b0, 0);
 532     src_b2 = LD_DP(b + 2);
 533     src_b3 = (v2f64) __msa_splati_d((v2i64) src_b2, 1);
 534     src_b2 = (v2f64) __msa_splati_d((v2i64) src_b2, 0);
 535
 536     src_c2 *= src_b3;
 537     src_c3 *= src_b3;
 538
 539     src_c0 -= src_c2 * src_b2;
 540     src_c1 -= src_c3 * src_b2;
 541
 542     src_c0 *= src_b0;
 543     src_c1 *= src_b0;
 544
 545     ST_DP2(src_c0, src_c1, c, 2);
 546     ST_DP2(src_c2, src_c3, c + ldc, 2);
 547
 548     ST_DP4(src_c0, src_c1, src_c2, src_c3, a, 2);
 549 }
 550
 551 static void dsolve_4x1_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
 552 {
 553     FLOAT b0, c0, c1, c2, c3;
 554
 555     c0 = *(c + 0);
 556     c1 = *(c + 1);
 557     c2 = *(c + 2);
 558     c3 = *(c + 3);
 559
 560     if (bk > 0)
 561     {
 562         BLASLONG i;
 563         FLOAT *aa = a, *bb = b;
 564
 565         for (i = bk; i--;)
 566         {
 567             c0 -= aa[0] * bb[0];
 568             c1 -= aa[1] * bb[0];
 569             c2 -= aa[2] * bb[0];
 570             c3 -= aa[3] * bb[0];
 571
 572             aa += 4;
 573             bb += 1;
 574         }
 575     }
 576
 577     a -= 4;
 578
 579     b0 = *(b - 1);
 580
 581     c0 *= b0;
 582     c1 *= b0;
 583     c2 *= b0;
 584     c3 *= b0;
 585
 586     *(a + 0) = c0;
 587     *(a + 1) = c1;
 588     *(a + 2) = c2;
 589     *(a + 3) = c3;
 590
 591     *(c + 0) = c0;
 592     *(c + 1) = c1;
 593     *(c + 2) = c2;
 594     *(c + 3) = c3;
 595 }
 596
 597 static void dsolve_2x4_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 598 {
 599     FLOAT b0, b4, b5, b8, b9, b10, b12, b13, b14, b15;
 600     FLOAT c0, c1, c0_nxt1, c1_nxt1, c0_nxt2, c1_nxt2, c0_nxt3, c1_nxt3;
 601
 602     c0 = *(c + 0);
 603     c1 = *(c + 1);
 604     c0_nxt1 = *(c + 0 + 1 * ldc);
 605     c1_nxt1 = *(c + 1 + 1 * ldc);
 606     c0_nxt2 = *(c + 0 + 2 * ldc);
 607     c1_nxt2 = *(c + 1 + 2 * ldc);
 608     c0_nxt3 = *(c + 0 + 3 * ldc);
 609     c1_nxt3 = *(c + 1 + 3 * ldc);
 610
 611     if (bk > 0)
 612     {
 613         BLASLONG i;
 614         FLOAT *aa = a, *bb = b;
 615
 616         for (i = bk; i--;)
 617         {
 618             c0 -= aa[0] * bb[0];
 619             c1 -= aa[1] * bb[0];
 620             c0_nxt1 -= aa[0] * bb[1];
 621             c1_nxt1 -= aa[1] * bb[1];
 622             c0_nxt2 -= aa[0] * bb[2];
 623             c1_nxt2 -= aa[1] * bb[2];
 624             c0_nxt3 -= aa[0] * bb[3];
 625             c1_nxt3 -= aa[1] * bb[3];
 626
 627             aa += 2;
 628             bb += 4;
 629         }
 630     }
 631
 632     a -= 8;
 633     b -= 16;
 634
 635     b0 = *b;
 636     b4 = *(b + 4);
 637     b5 = *(b + 5);
 638     b8 = *(b + 8);
 639     b9 = *(b + 9);
 640     b10 = *(b + 10);
 641     b12 = *(b + 12);
 642     b13 = *(b + 13);
 643     b14 = *(b + 14);
 644     b15 = *(b + 15);
 645
 646     c0_nxt3 *= b15;
 647     c1_nxt3 *= b15;
 648
 649     c0_nxt2 -= c0_nxt3 * b14;
 650     c1_nxt2 -= c1_nxt3 * b14;
 651     c0_nxt2 *= b10;
 652     c1_nxt2 *= b10;
 653
 654     c0_nxt1 -= c0_nxt3 * b13;
 655     c1_nxt1 -= c1_nxt3 * b13;
 656     c0_nxt1 -= c0_nxt2 * b9;
 657     c1_nxt1 -= c1_nxt2 * b9;
 658     c0_nxt1 *= b5;
 659     c1_nxt1 *= b5;
 660
 661     c0 -= c0_nxt3 * b12;
 662     c1 -= c1_nxt3 * b12;
 663     c0 -= c0_nxt2 * b8;
 664     c1 -= c1_nxt2 * b8;
 665     c0 -= c0_nxt1 * b4;
 666     c1 -= c1_nxt1 * b4;
 667     c0 *= b0;
 668     c1 *= b0;
 669
 670     *(a + 0) = c0;
 671     *(a + 1) = c1;
 672     *(a + 2) = c0_nxt1;
 673     *(a + 3) = c1_nxt1;
 674     *(a + 4) = c0_nxt2;
 675     *(a + 5) = c1_nxt2;
 676     *(a + 6) = c0_nxt3;
 677     *(a + 7) = c1_nxt3;
 678
 679     *(c + 0) = c0;
 680     *(c + 1) = c1;
 681     *(c + 0 + 1 * ldc) = c0_nxt1;
 682     *(c + 1 + 1 * ldc) = c1_nxt1;
 683     *(c + 0 + 2 * ldc) = c0_nxt2;
 684     *(c + 1 + 2 * ldc) = c1_nxt2;
 685     *(c + 0 + 3 * ldc) = c0_nxt3;
 686     *(c + 1 + 3 * ldc) = c1_nxt3;
 687 }
 688
 689 static void dsolve_2x2_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 690 {
 691     FLOAT b0, b2, b3, c0, c1, c0_nxt, c1_nxt;
 692
 693     c0 = *(c + 0);
 694     c1 = *(c + 1);
 695     c0_nxt = *(c + 0 + ldc);
 696     c1_nxt = *(c + 1 + ldc);
 697
 698     if (bk > 0)
 699     {
 700         BLASLONG i;
 701         FLOAT *aa = a, *bb = b;
 702
 703         for (i = bk; i--;)
 704         {
 705             c0 -= aa[0] * bb[0];
 706             c1 -= aa[1] * bb[0];
 707
 708             c0_nxt -= aa[0] * bb[1];
 709             c1_nxt -= aa[1] * bb[1];
 710
 711             aa += 2;
 712             bb += 2;
 713         }
 714     }
 715
 716     a -= 4;
 717     b -= 4;
 718
 719     b3 = *(b + 3);
 720     b2 = *(b + 2);
 721     b0 = *b;
 722
 723     c0_nxt *= b3;
 724     c1_nxt *= b3;
 725
 726     c0 -= c0_nxt * b2;
 727     c0 *= b0;
 728
 729     c1 -= c1_nxt * b2;
 730     c1 *= b0;
 731
 732     *(a + 0) = c0;
 733     *(a + 1) = c1;
 734     *(a + 2) = c0_nxt;
 735     *(a + 3) = c1_nxt;
 736
 737     *(c + 0) = c0;
 738     *(c + 1) = c1;
 739     *(c + 0 + ldc) = c0_nxt;
 740     *(c + 1 + ldc) = c1_nxt;
 741 }
 742
 743 static void dsolve_2x1_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
 744 {
 745     FLOAT b0, c0, c1;
 746
 747     c0 = *(c + 0);
 748     c1 = *(c + 1);
 749
 750     if (bk > 0)
 751     {
 752         BLASLONG i;
 753         FLOAT *aa = a, *bb = b;
 754
 755         for (i = bk; i--;)
 756         {
 757             c0 -= aa[0] * bb[0];
 758             c1 -= aa[1] * bb[0];
 759
 760             aa += 2;
 761             bb += 1;
 762         }
 763     }
 764
 765     b0 = *(b - 1);
 766
 767     c0 *= b0;
 768     c1 *= b0;
 769
 770     *(a - 2) = c0;
 771     *(a - 1) = c1;
 772
 773     *(c + 0) = c0;
 774     *(c + 1) = c1;
 775 }
 776
 777 static void dsolve_1x4_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 778 {
 779     FLOAT b0, b4, b5, b8, b9, b10, b12, b13, b14, b15, c0, c1, c2, c3;
 780
 781     c0 = *(c + 0);
 782     c1 = *(c + 1 * ldc);
 783     c2 = *(c + 2 * ldc);
 784     c3 = *(c + 3 * ldc);
 785
 786     if (bk > 0)
 787     {
 788         BLASLONG i;
 789         FLOAT *aa = a, *bb = b;
 790
 791         for (i = bk; i--;)
 792         {
 793             c0 -= aa[0] * bb[0];
 794             c1 -= aa[0] * bb[1];
 795             c2 -= aa[0] * bb[2];
 796             c3 -= aa[0] * bb[3];
 797
 798             aa += 1;
 799             bb += 4;
 800         }
 801     }
 802
 803     a -= 4;
 804     b -= 16;
 805
 806     b0 = *b;
 807     b4 = *(b + 4);
 808     b5 = *(b + 5);
 809     b8 = *(b + 8);
 810     b9 = *(b + 9);
 811     b10 = *(b + 10);
 812     b12 = *(b + 12);
 813     b13 = *(b + 13);
 814     b14 = *(b + 14);
 815     b15 = *(b + 15);
 816
 817     c3 *= b15;
 818
 819     c2 -= c3 * b14;
 820     c2 *= b10;
 821
 822     c1 -= c3 * b13;
 823     c1 -= c2 * b9;
 824     c1 *= b5;
 825
 826     c0 -= c3 * b12;
 827     c0 -= c2 * b8;
 828     c0 -= c1 * b4;
 829     c0 *= b0;
 830
 831     *(a + 0) = c0;
 832     *(a + 1) = c1;
 833     *(a + 2) = c2;
 834     *(a + 3) = c3;
 835
 836     *(c + 0 * ldc) = c0;
 837     *(c + 1 * ldc) = c1;
 838     *(c + 2 * ldc) = c2;
 839     *(c + 3 * ldc) = c3;
 840 }
 841
 842 static void dsolve_1x2_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 843 {
 844     FLOAT b0, b2, b3, c0, c1;
 845
 846     c0 = *(c + 0);
 847     c1 = *(c + ldc);
 848
 849     if (bk > 0)
 850     {
 851         BLASLONG i;
 852         FLOAT *aa = a, *bb = b;
 853
 854         for (i = bk; i--;)
 855         {
 856             c0 -= *aa * bb[0];
 857             c1 -= *aa * bb[1];
 858
 859             aa += 1;
 860             bb += 2;
 861         }
 862     }
 863
 864     a -= 2;
 865     b -= 4;
 866
 867     b3 = *(b + 3);
 868     b2 = *(b + 2);
 869     b0 = *b;
 870
 871     c1 *= b3;
 872
 873     c0 -= c1 * b2;
 874     c0 *= b0;
 875
 876     *(a + 0) = c0;
 877     *(a + 1) = c1;
 878
 879     *(c + 0) = c0;
 880     *(c + ldc) = c1;
 881 }
 882
 883 static void dsolve_1x1_rt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
 884 {
 885     if (bk > 0)
 886     {
 887         BLASLONG i;
 888
 889         for (i = 0; i < bk; i++)
 890         {
 891             *c -= a[i] * b[i];
 892         }
 893     }
 894
 895     *c *= *(b - 1);
 896     *(a - 1) = *c;
 897 }
 898
 899 int CNAME(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT dummy1, FLOAT *a, FLOAT *b,
 900           FLOAT *c, BLASLONG ldc, BLASLONG offset)
 901 {
 902     BLASLONG i, j, kk;
 903     FLOAT *aa, *cc, *bb;
 904
 905     kk = n - offset;
 906     c += n * ldc;
 907     b += n * k;
 908
 909     if (n & 3)
 910     {
 911         if (n & 1)
 912         {
 913             aa = a;
 914             c -= ldc;
 915             b -= k;
 916             bb = b + kk;
 917             cc = c;
 918
 919             for (i = (m >> 3); i--;)
 920             {
 921                 dsolve_8x1_rt_msa(aa + 8 * kk, bb, cc, (k - kk));
 922
 923                 aa += 8 * k;
 924                 cc += 8;
 925             }
 926
 927             if (m & 7)
 928             {
 929                 if (m & 4)
 930                 {
 931                     dsolve_4x1_rt_msa(aa + 4 * kk, bb, cc, (k - kk));
 932
 933                     aa += 4 * k;
 934                     cc += 4;
 935                 }
 936
 937                 if (m & 2)
 938                 {
 939                     dsolve_2x1_rt_msa(aa + 2 * kk, bb, cc, (k - kk));
 940
 941                     aa += 2 * k;
 942                     cc += 2;
 943                 }
 944
 945                 if (m & 1)
 946                 {
 947                     dsolve_1x1_rt_msa(aa + kk, bb, cc, (k - kk));
 948
 949                     aa += k;
 950                     cc += 1;
 951                 }
 952
 953             }
 954
 955             kk -= 1;
 956         }
 957
 958         if (n & 2)
 959         {
 960             aa = a;
 961             c -= 2 * ldc;
 962             b -= 2 * k;
 963             bb = b + 2 * kk;
 964             cc = c;
 965
 966             for (i = (m >> 3); i--;)
 967             {
 968                 dsolve_8x2_rt_msa(aa + 8 * kk, bb, cc, ldc, (k - kk));
 969
 970                 aa += 8 * k;
 971                 cc += 8;
 972             }
 973
 974             if (m & 7)
 975             {
 976                 if (m & 4)
 977                 {
 978                     dsolve_4x2_rt_msa(aa + 4 * kk, bb, cc, ldc, (k - kk));
 979
 980                     aa += 4 * k;
 981                     cc += 4;
 982                 }
 983
 984                 if (m & 2)
 985                 {
 986                     dsolve_2x2_rt_msa(aa + 2 * kk, bb, cc, ldc, (k - kk));
 987
 988                     aa += 2 * k;
 989                     cc += 2;
 990                 }
 991
 992                 if (m & 1)
 993                 {
 994                     dsolve_1x2_rt_msa(aa + kk, bb, cc, ldc, (k - kk));
 995
 996                     aa += k;
 997                     cc += 1;
 998                 }
 999             }
1000
1001             kk -= 2;
1002         }
1003     }
1004
1005     for (j = (n >> 2); j--;)
1006     {
1007         aa  = a;
1008         b -= 4 * k;
1009         bb = b + 4 * kk;
1010         c -= 4 * ldc;
1011         cc = c;
1012
1013         for (i = (m >> 3); i--;)
1014         {
1015             dsolve_8x4_rt_msa(aa + kk * 8, bb, cc, ldc, (k - kk));
1016
1017             aa += 8 * k;
1018             cc += 8;
1019         }
1020
1021         if (m & 7)
1022         {
1023             if (m & 4)
1024             {
1025                 dsolve_4x4_rt_msa(aa + kk * 4, bb, cc, ldc, (k - kk));
1026
1027                 aa += 4 * k;
1028                 cc += 4;
1029             }
1030
1031             if (m & 2)
1032             {
1033                 dsolve_2x4_rt_msa(aa + kk * 2, bb, cc, ldc, (k - kk));
1034
1035                 aa += 2 * k;
1036                 cc += 2;
1037             }
1038
1039             if (m & 1)
1040             {
1041                 dsolve_1x4_rt_msa(aa + kk, bb, cc, ldc, (k - kk));
1042
1043                 aa += k;
1044                 cc += 1;
1045             }
1046         }
1047
1048         kk -= 4;
1049     }
1050
1051     return 0;
1052 }