kernel/mips/dtrsm_kernel_LT_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_lt_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 res_c0, res_c1, res_c2, res_c3, res_c4, res_c5, res_c6, res_c7;
  37     v2f64 res_c8, res_c9, res_c10, res_c11, res_c12, res_c13, res_c14, res_c15;
  38     v2f64 src_a0, src_a1, src_a2, src_a3, src_a4, src_a5, src_a6, src_a7;
  39     v2f64 src_a9, src_a10, src_a11, src_a12, src_a13, src_a14, src_a15, src_a18;
  40     v2f64 src_a19, src_a20, src_a21, src_a22, src_a23, src_a27, src_a28;
  41     v2f64 src_a29, src_a30, src_a31, src_a36, src_a37, src_a38, src_a39;
  42     v2f64 src_a45, src_a46, src_a47, src_a54, src_a55, src_a63;
  43     FLOAT *c_nxt1line = c + ldc;
  44     FLOAT *c_nxt2line = c + 2 * ldc;
  45     FLOAT *c_nxt3line = c + 3 * ldc;
  46
  47     a += bk * 8;
  48     PREF_OFFSET(a,   0);
  49     PREF_OFFSET(a,  32);
  50     PREF_OFFSET(a,  72);
  51     PREF_OFFSET(a, 104);
  52     PREF_OFFSET(a, 144);
  53     PREF_OFFSET(a, 176);
  54     PREF_OFFSET(a, 216);
  55     PREF_OFFSET(a, 248);
  56     PREF_OFFSET(a, 288);
  57     PREF_OFFSET(a, 360);
  58     PREF_OFFSET(a, 504);
  59     PREF_OFFSET(a, 432);
  60     a -= bk * 8;
  61
  62     LD_DP4(c, 2, src_c0, src_c1, src_c2, src_c3);
  63     LD_DP4(c_nxt1line, 2, src_c4, src_c5, src_c6, src_c7);
  64     LD_DP4(c_nxt2line, 2, src_c8, src_c9, src_c10, src_c11);
  65     LD_DP4(c_nxt3line, 2, src_c12, src_c13, src_c14, src_c15);
  66
  67     if (bk)
  68     {
  69         BLASLONG i, pref_offset;
  70         FLOAT *pa0_pref;
  71         v2f64 src_b, src_b0, src_b1;
  72
  73         pref_offset = (uintptr_t)a & (L1_DATA_LINESIZE - 1);
  74
  75         if (pref_offset)
  76         {
  77             pref_offset = L1_DATA_LINESIZE - pref_offset;
  78             pref_offset = pref_offset / sizeof(FLOAT);
  79         }
  80
  81         pa0_pref = a + pref_offset;
  82
  83         for (i = (bk >> 1); i--;)
  84         {
  85             PREF_OFFSET(pa0_pref, 128);
  86             PREF_OFFSET(pa0_pref, 160);
  87             PREF_OFFSET(pa0_pref, 192);
  88             PREF_OFFSET(pa0_pref, 224);
  89
  90             LD_DP4_INC(a, 2, src_a0, src_a1, src_a2, src_a3);
  91             LD_DP2_INC(b, 2, src_b0, src_b1);
  92
  93             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
  94             src_c0 -= src_a0 * src_b;
  95             src_c1 -= src_a1 * src_b;
  96             src_c2 -= src_a2 * src_b;
  97             src_c3 -= src_a3 * src_b;
  98
  99             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 100             src_c4 -= src_a0 * src_b;
 101             src_c5 -= src_a1 * src_b;
 102             src_c6 -= src_a2 * src_b;
 103             src_c7 -= src_a3 * src_b;
 104
 105             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
 106             src_c8  -= src_a0 * src_b;
 107             src_c9  -= src_a1 * src_b;
 108             src_c10 -= src_a2 * src_b;
 109             src_c11 -= src_a3 * src_b;
 110
 111             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
 112             src_c12 -= src_a0 * src_b;
 113             src_c13 -= src_a1 * src_b;
 114             src_c14 -= src_a2 * src_b;
 115             src_c15 -= src_a3 * src_b;
 116
 117             LD_DP4_INC(a, 2, src_a0, src_a1, src_a2, src_a3);
 118             LD_DP2_INC(b, 2, src_b0, src_b1);
 119
 120             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 121             src_c0 -= src_a0 * src_b;
 122             src_c1 -= src_a1 * src_b;
 123             src_c2 -= src_a2 * src_b;
 124             src_c3 -= src_a3 * src_b;
 125
 126             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 127             src_c4 -= src_a0 * src_b;
 128             src_c5 -= src_a1 * src_b;
 129             src_c6 -= src_a2 * src_b;
 130             src_c7 -= src_a3 * src_b;
 131
 132             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
 133             src_c8  -= src_a0 * src_b;
 134             src_c9  -= src_a1 * src_b;
 135             src_c10 -= src_a2 * src_b;
 136             src_c11 -= src_a3 * src_b;
 137
 138             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
 139             src_c12 -= src_a0 * src_b;
 140             src_c13 -= src_a1 * src_b;
 141             src_c14 -= src_a2 * src_b;
 142             src_c15 -= src_a3 * src_b;
 143
 144             pa0_pref += 16;
 145         }
 146
 147         if (bk & 1)
 148         {
 149             LD_DP4_INC(a, 2, src_a0, src_a1, src_a2, src_a3);
 150             LD_DP2_INC(b, 2, src_b0, src_b1);
 151
 152             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 153             src_c0 -= src_a0 * src_b;
 154             src_c1 -= src_a1 * src_b;
 155             src_c2 -= src_a2 * src_b;
 156             src_c3 -= src_a3 * src_b;
 157
 158             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 159             src_c4 -= src_a0 * src_b;
 160             src_c5 -= src_a1 * src_b;
 161             src_c6 -= src_a2 * src_b;
 162             src_c7 -= src_a3 * src_b;
 163
 164             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
 165             src_c8  -= src_a0 * src_b;
 166             src_c9  -= src_a1 * src_b;
 167             src_c10 -= src_a2 * src_b;
 168             src_c11 -= src_a3 * src_b;
 169
 170             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
 171             src_c12 -= src_a0 * src_b;
 172             src_c13 -= src_a1 * src_b;
 173             src_c14 -= src_a2 * src_b;
 174             src_c15 -= src_a3 * src_b;
 175         }
 176     }
 177
 178     ILVRL_D2_DP(src_c4, src_c0, res_c0, res_c1);
 179     ILVRL_D2_DP(src_c5, src_c1, res_c2, res_c3);
 180     ILVRL_D2_DP(src_c6, src_c2, res_c4, res_c5);
 181     ILVRL_D2_DP(src_c7, src_c3, res_c6, res_c7);
 182     ILVRL_D2_DP(src_c12, src_c8, res_c8, res_c9);
 183     ILVRL_D2_DP(src_c13, src_c9, res_c10, res_c11);
 184     ILVRL_D2_DP(src_c14, src_c10, res_c12, res_c13);
 185     ILVRL_D2_DP(src_c15, src_c11, res_c14, res_c15);
 186
 187     src_a0 = LD_DP(a + 0);
 188     src_a1 = (v2f64) __msa_splati_d((v2i64) src_a0, 1);
 189     src_a0 = (v2f64) __msa_splati_d((v2i64) src_a0, 0);
 190     src_a2 = LD_DP(a + 2);
 191     src_a3 = (v2f64) __msa_splati_d((v2i64) src_a2, 1);
 192     src_a2 = (v2f64) __msa_splati_d((v2i64) src_a2, 0);
 193     src_a4 = LD_DP(a + 4);
 194     src_a5 = (v2f64) __msa_splati_d((v2i64) src_a4, 1);
 195     src_a4 = (v2f64) __msa_splati_d((v2i64) src_a4, 0);
 196     src_a6 = LD_DP(a + 6);
 197     src_a7 = (v2f64) __msa_splati_d((v2i64) src_a6, 1);
 198     src_a6 = (v2f64) __msa_splati_d((v2i64) src_a6, 0);
 199
 200     res_c0 *= src_a0;
 201     res_c1 -= res_c0 * src_a1;
 202     res_c2 -= res_c0 * src_a2;
 203     res_c3 -= res_c0 * src_a3;
 204     res_c4 -= res_c0 * src_a4;
 205     res_c5 -= res_c0 * src_a5;
 206     res_c6 -= res_c0 * src_a6;
 207     res_c7 -= res_c0 * src_a7;
 208
 209     res_c8 *= src_a0;
 210     res_c9 -= res_c8 * src_a1;
 211     res_c10 -= res_c8 * src_a2;
 212     res_c11 -= res_c8 * src_a3;
 213     res_c12 -= res_c8 * src_a4;
 214     res_c13 -= res_c8 * src_a5;
 215     res_c14 -= res_c8 * src_a6;
 216     res_c15 -= res_c8 * src_a7;
 217
 218     src_a9 = __msa_cast_to_vector_double(*(a + 9));
 219     src_a9 = (v2f64) __msa_splati_d((v2i64) src_a9, 0);
 220     src_a10 = LD_DP(a + 10);
 221     src_a11 = (v2f64) __msa_splati_d((v2i64) src_a10, 1);
 222     src_a10 = (v2f64) __msa_splati_d((v2i64) src_a10, 0);
 223     src_a12 = LD_DP(a + 12);
 224     src_a13 = (v2f64) __msa_splati_d((v2i64) src_a12, 1);
 225     src_a12 = (v2f64) __msa_splati_d((v2i64) src_a12, 0);
 226     src_a14 = LD_DP(a + 14);
 227     src_a15 = (v2f64) __msa_splati_d((v2i64) src_a14, 1);
 228     src_a14 = (v2f64) __msa_splati_d((v2i64) src_a14, 0);
 229
 230     res_c1 *= src_a9;
 231     res_c2 -= res_c1 * src_a10;
 232     res_c3 -= res_c1 * src_a11;
 233     res_c4 -= res_c1 * src_a12;
 234     res_c5 -= res_c1 * src_a13;
 235     res_c6 -= res_c1 * src_a14;
 236     res_c7 -= res_c1 * src_a15;
 237
 238     res_c9 *= src_a9;
 239     res_c10 -= res_c9 * src_a10;
 240     res_c11 -= res_c9 * src_a11;
 241     res_c12 -= res_c9 * src_a12;
 242     res_c13 -= res_c9 * src_a13;
 243     res_c14 -= res_c9 * src_a14;
 244     res_c15 -= res_c9 * src_a15;
 245
 246     ST_DP(res_c0, b + 0);
 247     ST_DP(res_c8, b + 2);
 248     ST_DP(res_c1, b + 4);
 249     ST_DP(res_c9, b + 6);
 250
 251     ILVRL_D2_DP(res_c1, res_c0, src_c0, src_c4);
 252     ILVRL_D2_DP(res_c9, res_c8, src_c8, src_c12);
 253
 254     ST_DP(src_c0, c);
 255     ST_DP(src_c4, c_nxt1line);
 256     ST_DP(src_c8, c_nxt2line);
 257     ST_DP(src_c12, c_nxt3line);
 258
 259     src_a18 = LD_DP(a + 18);
 260     src_a19 = (v2f64) __msa_splati_d((v2i64) src_a18, 1);
 261     src_a18 = (v2f64) __msa_splati_d((v2i64) src_a18, 0);
 262     src_a20 = LD_DP(a + 20);
 263     src_a21 = (v2f64) __msa_splati_d((v2i64) src_a20, 1);
 264     src_a20 = (v2f64) __msa_splati_d((v2i64) src_a20, 0);
 265     src_a22 = LD_DP(a + 22);
 266     src_a23 = (v2f64) __msa_splati_d((v2i64) src_a22, 1);
 267     src_a22 = (v2f64) __msa_splati_d((v2i64) src_a22, 0);
 268
 269     res_c2 *= src_a18;
 270     res_c3 -= res_c2 * src_a19;
 271     res_c4 -= res_c2 * src_a20;
 272     res_c5 -= res_c2 * src_a21;
 273     res_c6 -= res_c2 * src_a22;
 274     res_c7 -= res_c2 * src_a23;
 275
 276     res_c10 *= src_a18;
 277     res_c11 -= res_c10 * src_a19;
 278     res_c12 -= res_c10 * src_a20;
 279     res_c13 -= res_c10 * src_a21;
 280     res_c14 -= res_c10 * src_a22;
 281     res_c15 -= res_c10 * src_a23;
 282
 283     src_a27 = __msa_cast_to_vector_double(*(a + 27));
 284     src_a27 = (v2f64) __msa_splati_d((v2i64) src_a27, 0);
 285     src_a28 = LD_DP(a + 28);
 286     src_a29 = (v2f64) __msa_splati_d((v2i64) src_a28, 1);
 287     src_a28 = (v2f64) __msa_splati_d((v2i64) src_a28, 0);
 288     src_a30 = LD_DP(a + 30);
 289     src_a31 = (v2f64) __msa_splati_d((v2i64) src_a30, 1);
 290     src_a30 = (v2f64) __msa_splati_d((v2i64) src_a30, 0);
 291
 292     res_c3 *= src_a27;
 293     res_c4 -= res_c3 * src_a28;
 294     res_c5 -= res_c3 * src_a29;
 295     res_c6 -= res_c3 * src_a30;
 296     res_c7 -= res_c3 * src_a31;
 297
 298     res_c11 *= src_a27;
 299     res_c12 -= res_c11 * src_a28;
 300     res_c13 -= res_c11 * src_a29;
 301     res_c14 -= res_c11 * src_a30;
 302     res_c15 -= res_c11 * src_a31;
 303
 304     ST_DP(res_c2, b + 8);
 305     ST_DP(res_c10, b + 10);
 306     ST_DP(res_c3, b + 12);
 307     ST_DP(res_c11, b + 14);
 308
 309     ILVRL_D2_DP(res_c3, res_c2, src_c1, src_c5);
 310     ILVRL_D2_DP(res_c11, res_c10, src_c9, src_c13);
 311
 312     src_a36 = LD_DP(a + 36);
 313     src_a37 = (v2f64) __msa_splati_d((v2i64) src_a36, 1);
 314     src_a36 = (v2f64) __msa_splati_d((v2i64) src_a36, 0);
 315     src_a38 = LD_DP(a + 38);
 316     src_a39 = (v2f64) __msa_splati_d((v2i64) src_a38, 1);
 317     src_a38 = (v2f64) __msa_splati_d((v2i64) src_a38, 0);
 318
 319     res_c4 *= src_a36;
 320     res_c5 -= res_c4 * src_a37;
 321     res_c6 -= res_c4 * src_a38;
 322     res_c7 -= res_c4 * src_a39;
 323
 324     res_c12 *= src_a36;
 325     res_c13 -= res_c12 * src_a37;
 326     res_c14 -= res_c12 * src_a38;
 327     res_c15 -= res_c12 * src_a39;
 328
 329     src_a45 = __msa_cast_to_vector_double(*(a + 45));
 330     src_a45 = (v2f64) __msa_splati_d((v2i64) src_a45, 0);
 331     src_a46 = LD_DP(a + 46);
 332     src_a47 = (v2f64) __msa_splati_d((v2i64) src_a46, 1);
 333     src_a46 = (v2f64) __msa_splati_d((v2i64) src_a46, 0);
 334
 335     res_c5 *= src_a45;
 336     res_c6 -= res_c5 * src_a46;
 337     res_c7 -= res_c5 * src_a47;
 338
 339     res_c13 *= src_a45;
 340     res_c14 -= res_c13 * src_a46;
 341     res_c15 -= res_c13 * src_a47;
 342
 343     ST_DP(src_c1, c + 2);
 344     ST_DP(src_c5, c_nxt1line + 2);
 345     ST_DP(src_c9, c_nxt2line + 2);
 346     ST_DP(src_c13, c_nxt3line + 2);
 347
 348     ST_DP(res_c4, b + 16);
 349     ST_DP(res_c12, b + 18);
 350     ST_DP(res_c5, b + 20);
 351     ST_DP(res_c13, b + 22);
 352
 353     ILVRL_D2_DP(res_c5, res_c4, src_c2, src_c6);
 354     ILVRL_D2_DP(res_c13, res_c12, src_c10, src_c14);
 355
 356     src_a63 = __msa_cast_to_vector_double(*(a + 63));
 357     src_a63 = (v2f64) __msa_splati_d((v2i64) src_a63, 0);
 358     src_a54 = LD_DP(a + 54);
 359     src_a55 = (v2f64) __msa_splati_d((v2i64) src_a54, 1);
 360     src_a54 = (v2f64) __msa_splati_d((v2i64) src_a54, 0);
 361
 362     res_c6 *= src_a54;
 363     res_c7 -= res_c6 * src_a55;
 364
 365     res_c14 *= src_a54;
 366     res_c15 -= res_c14 * src_a55;
 367
 368     res_c7 *= src_a63;
 369     res_c15 *= src_a63;
 370
 371     ST_DP(src_c2, c + 4);
 372     ST_DP(src_c6, c_nxt1line + 4);
 373     ST_DP(src_c10, c_nxt2line + 4);
 374     ST_DP(src_c14, c_nxt3line + 4);
 375
 376     ST_DP(res_c6, b + 24);
 377     ST_DP(res_c14, b + 26);
 378     ST_DP(res_c7, b + 28);
 379     ST_DP(res_c15, b + 30);
 380
 381     ILVRL_D2_DP(res_c7, res_c6, src_c3, src_c7);
 382     ILVRL_D2_DP(res_c15, res_c14, src_c11, src_c15);
 383
 384     ST_DP(src_c3, c + 6);
 385     ST_DP(src_c7, c_nxt1line + 6);
 386     ST_DP(src_c11, c_nxt2line + 6);
 387     ST_DP(src_c15, c_nxt3line + 6);
 388 }
 389
 390 static void dsolve_8x2_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 391 {
 392     v2f64 src_c0, src_c1, src_c2, src_c3, src_c4, src_c5, src_c6, src_c7;
 393     v2f64 res_c0, res_c1, res_c2, res_c3, res_c4, res_c5, res_c6, res_c7;
 394     v2f64 src_a0, src_a1, src_a2, src_a3, src_a4, src_a5, src_a6, src_a7;
 395     v2f64 src_a9, src_a10, src_a11, src_a12, src_a13, src_a14, src_a15, src_a18;
 396     v2f64 src_a19, src_a20, src_a21, src_a22, src_a23, src_a27, src_a28;
 397     v2f64 src_a29, src_a30, src_a31, src_a36, src_a37, src_a38, src_a39;
 398     v2f64 src_a45, src_a46, src_a47, src_a54, src_a55, src_a63;
 399
 400     LD_DP4(c, 2, src_c0, src_c1, src_c2, src_c3);
 401     LD_DP4(c + ldc, 2, src_c4, src_c5, src_c6, src_c7);
 402
 403     if (bk)
 404     {
 405         BLASLONG i;
 406         v2f64 src_b, src_b0, src_b1;
 407
 408         LD_DP4(a, 2, src_a0, src_a1, src_a2, src_a3);
 409         src_b0 = LD_DP(b);
 410
 411         a += 8;
 412         b += 2;
 413
 414         for (i = (bk - 1); i--;)
 415         {
 416             LD_DP4(a, 2, src_a4, src_a5, src_a6, src_a7);
 417             src_b1 = LD_DP(b);
 418
 419             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 420             src_c0 -= src_a0 * src_b;
 421             src_c1 -= src_a1 * src_b;
 422             src_c2 -= src_a2 * src_b;
 423             src_c3 -= src_a3 * src_b;
 424
 425             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 426             src_c4 -= src_a0 * src_b;
 427             src_c5 -= src_a1 * src_b;
 428             src_c6 -= src_a2 * src_b;
 429             src_c7 -= src_a3 * src_b;
 430
 431             src_a0 = src_a4;
 432             src_a1 = src_a5;
 433             src_a2 = src_a6;
 434             src_a3 = src_a7;
 435             src_b0 = src_b1;
 436
 437             a += 8;
 438             b += 2;
 439         }
 440
 441         src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 442         src_c0 -= src_a0 * src_b;
 443         src_c1 -= src_a1 * src_b;
 444         src_c2 -= src_a2 * src_b;
 445         src_c3 -= src_a3 * src_b;
 446
 447         src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 448         src_c4 -= src_a0 * src_b;
 449         src_c5 -= src_a1 * src_b;
 450         src_c6 -= src_a2 * src_b;
 451         src_c7 -= src_a3 * src_b;
 452     }
 453
 454     ILVRL_D2_DP(src_c4, src_c0, res_c0, res_c1);
 455     ILVRL_D2_DP(src_c5, src_c1, res_c2, res_c3);
 456     ILVRL_D2_DP(src_c6, src_c2, res_c4, res_c5);
 457     ILVRL_D2_DP(src_c7, src_c3, res_c6, res_c7);
 458
 459     src_a0 = LD_DP(a + 0);
 460     src_a1 = (v2f64) __msa_splati_d((v2i64) src_a0, 1);
 461     src_a0 = (v2f64) __msa_splati_d((v2i64) src_a0, 0);
 462     src_a2 = LD_DP(a + 2);
 463     src_a3 = (v2f64) __msa_splati_d((v2i64) src_a2, 1);
 464     src_a2 = (v2f64) __msa_splati_d((v2i64) src_a2, 0);
 465     src_a4 = LD_DP(a + 4);
 466     src_a5 = (v2f64) __msa_splati_d((v2i64) src_a4, 1);
 467     src_a4 = (v2f64) __msa_splati_d((v2i64) src_a4, 0);
 468     src_a6 = LD_DP(a + 6);
 469     src_a7 = (v2f64) __msa_splati_d((v2i64) src_a6, 1);
 470     src_a6 = (v2f64) __msa_splati_d((v2i64) src_a6, 0);
 471
 472     res_c0 *= src_a0;
 473     res_c1 -= res_c0 * src_a1;
 474     res_c2 -= res_c0 * src_a2;
 475     res_c3 -= res_c0 * src_a3;
 476     res_c4 -= res_c0 * src_a4;
 477     res_c5 -= res_c0 * src_a5;
 478     res_c6 -= res_c0 * src_a6;
 479     res_c7 -= res_c0 * src_a7;
 480
 481     src_a9 = __msa_cast_to_vector_double(*(a + 9));
 482     src_a9 = (v2f64) __msa_splati_d((v2i64) src_a9, 0);
 483     src_a10 = LD_DP(a + 10);
 484     src_a11 = (v2f64) __msa_splati_d((v2i64) src_a10, 1);
 485     src_a10 = (v2f64) __msa_splati_d((v2i64) src_a10, 0);
 486     src_a12 = LD_DP(a + 12);
 487     src_a13 = (v2f64) __msa_splati_d((v2i64) src_a12, 1);
 488     src_a12 = (v2f64) __msa_splati_d((v2i64) src_a12, 0);
 489     src_a14 = LD_DP(a + 14);
 490     src_a15 = (v2f64) __msa_splati_d((v2i64) src_a14, 1);
 491     src_a14 = (v2f64) __msa_splati_d((v2i64) src_a14, 0);
 492
 493     res_c1 *= src_a9;
 494     res_c2 -= res_c1 * src_a10;
 495     res_c3 -= res_c1 * src_a11;
 496     res_c4 -= res_c1 * src_a12;
 497     res_c5 -= res_c1 * src_a13;
 498     res_c6 -= res_c1 * src_a14;
 499     res_c7 -= res_c1 * src_a15;
 500
 501     src_a18 = LD_DP(a + 18);
 502     src_a19 = (v2f64) __msa_splati_d((v2i64) src_a18, 1);
 503     src_a18 = (v2f64) __msa_splati_d((v2i64) src_a18, 0);
 504     src_a20 = LD_DP(a + 20);
 505     src_a21 = (v2f64) __msa_splati_d((v2i64) src_a20, 1);
 506     src_a20 = (v2f64) __msa_splati_d((v2i64) src_a20, 0);
 507     src_a22 = LD_DP(a + 22);
 508     src_a23 = (v2f64) __msa_splati_d((v2i64) src_a22, 1);
 509     src_a22 = (v2f64) __msa_splati_d((v2i64) src_a22, 0);
 510
 511     res_c2 *= src_a18;
 512     res_c3 -= res_c2 * src_a19;
 513     res_c4 -= res_c2 * src_a20;
 514     res_c5 -= res_c2 * src_a21;
 515     res_c6 -= res_c2 * src_a22;
 516     res_c7 -= res_c2 * src_a23;
 517
 518     src_a27 = __msa_cast_to_vector_double(*(a + 27));
 519     src_a27 = (v2f64) __msa_splati_d((v2i64) src_a27, 0);
 520     src_a28 = LD_DP(a + 28);
 521     src_a29 = (v2f64) __msa_splati_d((v2i64) src_a28, 1);
 522     src_a28 = (v2f64) __msa_splati_d((v2i64) src_a28, 0);
 523     src_a30 = LD_DP(a + 30);
 524     src_a31 = (v2f64) __msa_splati_d((v2i64) src_a30, 1);
 525     src_a30 = (v2f64) __msa_splati_d((v2i64) src_a30, 0);
 526
 527     res_c3 *= src_a27;
 528     res_c4 -= res_c3 * src_a28;
 529     res_c5 -= res_c3 * src_a29;
 530     res_c6 -= res_c3 * src_a30;
 531     res_c7 -= res_c3 * src_a31;
 532
 533     ST_DP(res_c0, b + 0);
 534     ST_DP(res_c1, b + 2);
 535     ST_DP(res_c2, b + 4);
 536     ST_DP(res_c3, b + 6);
 537
 538     ILVRL_D2_DP(res_c1, res_c0, src_c0, src_c4);
 539     ILVRL_D2_DP(res_c3, res_c2, src_c1, src_c5);
 540
 541     ST_DP2(src_c0, src_c1, c, 2);
 542     ST_DP2(src_c4, src_c5, c + ldc, 2);
 543
 544     src_a36 = LD_DP(a + 36);
 545     src_a37 = (v2f64) __msa_splati_d((v2i64) src_a36, 1);
 546     src_a36 = (v2f64) __msa_splati_d((v2i64) src_a36, 0);
 547     src_a38 = LD_DP(a + 38);
 548     src_a39 = (v2f64) __msa_splati_d((v2i64) src_a38, 1);
 549     src_a38 = (v2f64) __msa_splati_d((v2i64) src_a38, 0);
 550
 551     res_c4 *= src_a36;
 552     res_c5 -= res_c4 * src_a37;
 553     res_c6 -= res_c4 * src_a38;
 554     res_c7 -= res_c4 * src_a39;
 555
 556     src_a45 = __msa_cast_to_vector_double(*(a + 45));
 557     src_a45 = (v2f64) __msa_splati_d((v2i64) src_a45, 0);
 558     src_a46 = LD_DP(a + 46);
 559     src_a47 = (v2f64) __msa_splati_d((v2i64) src_a46, 1);
 560     src_a46 = (v2f64) __msa_splati_d((v2i64) src_a46, 0);
 561
 562     res_c5 *= src_a45;
 563     res_c6 -= res_c5 * src_a46;
 564     res_c7 -= res_c5 * src_a47;
 565
 566     src_a63 = __msa_cast_to_vector_double(*(a + 63));
 567     src_a63 = (v2f64) __msa_splati_d((v2i64) src_a63, 0);
 568     src_a54 = LD_DP(a + 54);
 569     src_a55 = (v2f64) __msa_splati_d((v2i64) src_a54, 1);
 570     src_a54 = (v2f64) __msa_splati_d((v2i64) src_a54, 0);
 571
 572     res_c6 *= src_a54;
 573     res_c7 -= res_c6 * src_a55;
 574
 575     res_c7 *= src_a63;
 576
 577     ST_DP(res_c4, b + 8);
 578     ST_DP(res_c5, b + 10);
 579     ST_DP(res_c6, b + 12);
 580     ST_DP(res_c7, b + 14);
 581
 582     ILVRL_D2_DP(res_c5, res_c4, src_c2, src_c6);
 583     ILVRL_D2_DP(res_c7, res_c6, src_c3, src_c7);
 584
 585     ST_DP2(src_c2, src_c3, c + 4, 2);
 586     ST_DP2(src_c6, src_c7, c + 4 + ldc, 2);
 587 }
 588
 589 static void dsolve_8x1_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
 590 {
 591     FLOAT a0, a1, a2, a3, a4, a5, a6, a7, a9, a10, a11, a12, a13, a14, a15, a18;
 592     FLOAT a19, a20, a21, a22, a23, a27, a28, a29, a30, a31, a36, a37, a38, a39;
 593     FLOAT a45, a46, a47, a54, a55, a63, c0, c1, c2, c3, c4, c5, c6, c7;
 594
 595     c0 = *(c + 0);
 596     c1 = *(c + 1);
 597     c2 = *(c + 2);
 598     c3 = *(c + 3);
 599     c4 = *(c + 4);
 600     c5 = *(c + 5);
 601     c6 = *(c + 6);
 602     c7 = *(c + 7);
 603
 604     if (bk)
 605     {
 606         BLASLONG i;
 607
 608         for (i = bk; i--; )
 609         {
 610             c0 -= a[0] * b[0];
 611             c1 -= a[1] * b[0];
 612             c2 -= a[2] * b[0];
 613             c3 -= a[3] * b[0];
 614             c4 -= a[4] * b[0];
 615             c5 -= a[5] * b[0];
 616             c6 -= a[6] * b[0];
 617             c7 -= a[7] * b[0];
 618
 619             a += 8;
 620             b += 1;
 621         }
 622     }
 623
 624     a0 = *(a + 0);
 625     a1 = *(a + 1);
 626     a2 = *(a + 2);
 627     a3 = *(a + 3);
 628     a4 = *(a + 4);
 629     a5 = *(a + 5);
 630     a6 = *(a + 6);
 631     a7 = *(a + 7);
 632     a9 = *(a + 9);
 633     a10 = *(a + 10);
 634     a11 = *(a + 11);
 635     a12 = *(a + 12);
 636     a13 = *(a + 13);
 637     a14 = *(a + 14);
 638     a15 = *(a + 15);
 639     a18 = *(a + 18);
 640     a19 = *(a + 19);
 641     a20 = *(a + 20);
 642     a21 = *(a + 21);
 643     a22 = *(a + 22);
 644     a23 = *(a + 23);
 645     a27 = *(a + 27);
 646     a28 = *(a + 28);
 647     a29 = *(a + 29);
 648     a30 = *(a + 30);
 649     a31 = *(a + 31);
 650     a36 = *(a + 36);
 651     a37 = *(a + 37);
 652     a38 = *(a + 38);
 653     a39 = *(a + 39);
 654     a45 = *(a + 45);
 655     a46 = *(a + 46);
 656     a47 = *(a + 47);
 657     a54 = *(a + 54);
 658     a55 = *(a + 55);
 659     a63 = *(a + 63);
 660
 661     c0 *= a0;
 662
 663     c1 -= c0 * a1;
 664     c1 *= a9;
 665
 666     c2 -= c0 * a2;
 667     c2 -= c1 * a10;
 668     c2 *= a18;
 669
 670     c3 -= c0 * a3;
 671     c3 -= c1 * a11;
 672     c3 -= c2 * a19;
 673     c3 *= a27;
 674
 675     c4 -= c0 * a4;
 676     c4 -= c1 * a12;
 677     c4 -= c2 * a20;
 678     c4 -= c3 * a28;
 679     c4 *= a36;
 680
 681     c5 -= c0 * a5;
 682     c5 -= c1 * a13;
 683     c5 -= c2 * a21;
 684     c5 -= c3 * a29;
 685     c5 -= c4 * a37;
 686     c5 *= a45;
 687
 688     c6 -= c0 * a6;
 689     c6 -= c1 * a14;
 690     c6 -= c2 * a22;
 691     c6 -= c3 * a30;
 692     c6 -= c4 * a38;
 693     c6 -= c5 * a46;
 694     c6 *= a54;
 695
 696     c7 -= c0 * a7;
 697     c7 -= c1 * a15;
 698     c7 -= c2 * a23;
 699     c7 -= c3 * a31;
 700     c7 -= c4 * a39;
 701     c7 -= c5 * a47;
 702     c7 -= c6 * a55;
 703     c7 *= a63;
 704
 705     *(c + 0) = c0;
 706     *(c + 1) = c1;
 707     *(c + 2) = c2;
 708     *(c + 3) = c3;
 709     *(c + 4) = c4;
 710     *(c + 5) = c5;
 711     *(c + 6) = c6;
 712     *(c + 7) = c7;
 713
 714     *(b + 0) = c0;
 715     *(b + 1) = c1;
 716     *(b + 2) = c2;
 717     *(b + 3) = c3;
 718     *(b + 4) = c4;
 719     *(b + 5) = c5;
 720     *(b + 6) = c6;
 721     *(b + 7) = c7;
 722 }
 723
 724 static void dsolve_4x4_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 725 {
 726     v2f64 src_c0, src_c1, src_c2, src_c3, src_c4, src_c5, src_c6, src_c7;
 727     v2f64 res_c0, res_c1, res_c2, res_c3, res_c4, res_c5, res_c6, res_c7;
 728     v2f64 src_a0, src_a1, src_a2, src_a3, src_a5, src_a6, src_a7;
 729     v2f64 src_a10, src_a11, src_a15;
 730
 731     LD_DP2(c, 2, src_c0, src_c1);
 732     LD_DP2(c + ldc, 2, src_c2, src_c3);
 733     LD_DP2(c + 2 * ldc, 2, src_c4, src_c5);
 734     LD_DP2(c + 3 * ldc, 2, src_c6, src_c7);
 735
 736     if (bk)
 737     {
 738         BLASLONG i;
 739         v2f64 src_a0, src_a1, src_b, src_b0, src_b1;
 740
 741         for (i = bk; i--;)
 742         {
 743             LD_DP2(a, 2, src_a0, src_a1);
 744             LD_DP2(b, 2, src_b0, src_b1);
 745
 746             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 747             src_c0 -= src_a0 * src_b;
 748             src_c1 -= src_a1 * src_b;
 749
 750             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 751             src_c2 -= src_a0 * src_b;
 752             src_c3 -= src_a1 * src_b;
 753
 754             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b1, (v2i64) src_b1);
 755             src_c4 -= src_a0 * src_b;
 756             src_c5 -= src_a1 * src_b;
 757
 758             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b1, (v2i64) src_b1);
 759             src_c6 -= src_a0 * src_b;
 760             src_c7 -= src_a1 * src_b;
 761
 762             a += 4;
 763             b += 4;
 764         }
 765     }
 766
 767     ILVRL_D2_DP(src_c2, src_c0, res_c0, res_c1);
 768     ILVRL_D2_DP(src_c3, src_c1, res_c2, res_c3);
 769     ILVRL_D2_DP(src_c6, src_c4, res_c4, res_c5);
 770     ILVRL_D2_DP(src_c7, src_c5, res_c6, res_c7);
 771
 772     src_a0 = LD_DP(a + 0);
 773     src_a1 = (v2f64) __msa_splati_d((v2i64) src_a0, 1);
 774     src_a0 = (v2f64) __msa_splati_d((v2i64) src_a0, 0);
 775     src_a2 = LD_DP(a + 2);
 776     src_a3 = (v2f64) __msa_splati_d((v2i64) src_a2, 1);
 777     src_a2 = (v2f64) __msa_splati_d((v2i64) src_a2, 0);
 778
 779     res_c0 *= src_a0;
 780     res_c1 -= res_c0 * src_a1;
 781     res_c2 -= res_c0 * src_a2;
 782     res_c3 -= res_c0 * src_a3;
 783
 784     res_c4 *= src_a0;
 785     res_c5 -= res_c4 * src_a1;
 786     res_c6 -= res_c4 * src_a2;
 787     res_c7 -= res_c4 * src_a3;
 788
 789     src_a5 = __msa_cast_to_vector_double(*(a + 5));
 790     src_a5 = (v2f64) __msa_splati_d((v2i64) src_a5, 0);
 791     src_a6 = LD_DP(a + 6);
 792     src_a7 = (v2f64) __msa_splati_d((v2i64) src_a6, 1);
 793     src_a6 = (v2f64) __msa_splati_d((v2i64) src_a6, 0);
 794
 795     res_c1 *= src_a5;
 796     res_c2 -= res_c1 * src_a6;
 797     res_c3 -= res_c1 * src_a7;
 798
 799     res_c5 *= src_a5;
 800     res_c6 -= res_c5 * src_a6;
 801     res_c7 -= res_c5 * src_a7;
 802
 803     src_a10 = LD_DP(a + 10);
 804     src_a11 = (v2f64) __msa_splati_d((v2i64) src_a10, 1);
 805     src_a10 = (v2f64) __msa_splati_d((v2i64) src_a10, 0);
 806     src_a15 = __msa_cast_to_vector_double(*(a + 15));
 807     src_a15 = (v2f64) __msa_splati_d((v2i64) src_a15, 0);
 808
 809     res_c2 *= src_a10;
 810     res_c3 -= res_c2 * src_a11;
 811     res_c3 *= src_a15;
 812
 813     res_c6 *= src_a10;
 814     res_c7 -= res_c6 * src_a11;
 815     res_c7 *= src_a15;
 816
 817     ST_DP(res_c0, b + 0);
 818     ST_DP(res_c4, b + 2);
 819     ST_DP(res_c1, b + 4);
 820     ST_DP(res_c5, b + 6);
 821     ST_DP(res_c2, b + 8);
 822     ST_DP(res_c6, b + 10);
 823     ST_DP(res_c3, b + 12);
 824     ST_DP(res_c7, b + 14);
 825
 826     ILVRL_D2_DP(res_c1, res_c0, src_c0, src_c2);
 827     ILVRL_D2_DP(res_c3, res_c2, src_c1, src_c3);
 828     ILVRL_D2_DP(res_c5, res_c4, src_c4, src_c6);
 829     ILVRL_D2_DP(res_c7, res_c6, src_c5, src_c7);
 830
 831     ST_DP2(src_c0, src_c1, c, 2);
 832     ST_DP2(src_c2, src_c3, c + ldc, 2);
 833     ST_DP2(src_c4, src_c5, c + 2 * ldc, 2);
 834     ST_DP2(src_c6, src_c7, c + 3 * ldc, 2);
 835 }
 836
 837 static void dsolve_4x2_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 838 {
 839     v2f64 src_c0, src_c1, src_c2, src_c3, res_c0, res_c1, res_c2, res_c3;
 840     v2f64 src_a0, src_a1, src_a2, src_a3, src_a5, src_a6, src_a7;
 841     v2f64 src_a10, src_a11, src_a15;
 842
 843     LD_DP2(c, 2, src_c0, src_c1);
 844     LD_DP2(c + ldc, 2, src_c2, src_c3);
 845
 846     if (bk)
 847     {
 848         BLASLONG i;
 849         v2f64 src_a0, src_a1, src_b, src_b0;
 850
 851         for (i = bk; i--;)
 852         {
 853             LD_DP2(a, 2, src_a0, src_a1);
 854             src_b0 = LD_DP(b);
 855
 856             src_b = (v2f64) __msa_ilvr_d((v2i64) src_b0, (v2i64) src_b0);
 857             src_c0 -= src_a0 * src_b;
 858             src_c1 -= src_a1 * src_b;
 859
 860             src_b = (v2f64) __msa_ilvl_d((v2i64) src_b0, (v2i64) src_b0);
 861             src_c2 -= src_a0 * src_b;
 862             src_c3 -= src_a1 * src_b;
 863
 864             a += 4;
 865             b += 2;
 866         }
 867     }
 868
 869     ILVRL_D2_DP(src_c2, src_c0, res_c0, res_c1);
 870     ILVRL_D2_DP(src_c3, src_c1, res_c2, res_c3);
 871
 872     src_a0 = LD_DP(a + 0);
 873     src_a1 = (v2f64) __msa_splati_d((v2i64) src_a0, 1);
 874     src_a0 = (v2f64) __msa_splati_d((v2i64) src_a0, 0);
 875     src_a2 = LD_DP(a + 2);
 876     src_a3 = (v2f64) __msa_splati_d((v2i64) src_a2, 1);
 877     src_a2 = (v2f64) __msa_splati_d((v2i64) src_a2, 0);
 878
 879     res_c0 *= src_a0;
 880     res_c1 -= res_c0 * src_a1;
 881     res_c2 -= res_c0 * src_a2;
 882     res_c3 -= res_c0 * src_a3;
 883
 884     src_a5 = __msa_cast_to_vector_double(*(a + 5));
 885     src_a5 = (v2f64) __msa_splati_d((v2i64) src_a5, 0);
 886     src_a6 = LD_DP(a + 6);
 887     src_a7 = (v2f64) __msa_splati_d((v2i64) src_a6, 1);
 888     src_a6 = (v2f64) __msa_splati_d((v2i64) src_a6, 0);
 889
 890     res_c1 *= src_a5;
 891     res_c2 -= res_c1 * src_a6;
 892     res_c3 -= res_c1 * src_a7;
 893
 894     src_a10 = LD_DP(a + 10);
 895     src_a11 = (v2f64) __msa_splati_d((v2i64) src_a10, 1);
 896     src_a10 = (v2f64) __msa_splati_d((v2i64) src_a10, 0);
 897     src_a15 = __msa_cast_to_vector_double(*(a + 15));
 898     src_a15 = (v2f64) __msa_splati_d((v2i64) src_a15, 0);
 899
 900     res_c2 *= src_a10;
 901     res_c3 -= res_c2 * src_a11;
 902     res_c3 *= src_a15;
 903
 904     ST_DP(res_c0, b + 0);
 905     ST_DP(res_c1, b + 2);
 906     ST_DP(res_c2, b + 4);
 907     ST_DP(res_c3, b + 6);
 908
 909     ILVRL_D2_DP(res_c1, res_c0, src_c0, src_c2);
 910     ILVRL_D2_DP(res_c3, res_c2, src_c1, src_c3);
 911
 912     ST_DP2(src_c0, src_c1, c, 2);
 913     ST_DP2(src_c2, src_c3, c + ldc, 2);
 914 }
 915
 916 static void dsolve_4x1_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
 917 {
 918     FLOAT a0, a1, a2, a3, a5, a6, a7, a10, a11, a15, c0, c1, c2, c3;
 919
 920     c0 = *(c + 0);
 921     c1 = *(c + 1);
 922     c2 = *(c + 2);
 923     c3 = *(c + 3);
 924
 925     if (bk)
 926     {
 927         BLASLONG i;
 928
 929         for (i = bk; i--;)
 930         {
 931             c0 -= a[0] * b[0];
 932             c1 -= a[1] * b[0];
 933             c2 -= a[2] * b[0];
 934             c3 -= a[3] * b[0];
 935
 936             a += 4;
 937             b += 1;
 938         }
 939     }
 940
 941     a0 = *(a + 0);
 942     a1 = *(a + 1);
 943     a2 = *(a + 2);
 944     a3 = *(a + 3);
 945     a5 = *(a + 5);
 946     a6 = *(a + 6);
 947     a7 = *(a + 7);
 948     a10 = *(a + 10);
 949     a11 = *(a + 11);
 950     a15 = *(a + 15);
 951
 952     c0 *= a0;
 953
 954     c1 -= c0 * a1;
 955     c1 *= a5;
 956
 957     c2 -= c0 * a2;
 958     c2 -= c1 * a6;
 959     c2 *= a10;
 960
 961     c3 -= c0 * a3;
 962     c3 -= c1 * a7;
 963     c3 -= c2 * a11;
 964     c3 *= a15;
 965
 966     *(b + 0) = c0;
 967     *(b + 1) = c1;
 968     *(b + 2) = c2;
 969     *(b + 3) = c3;
 970
 971     *(c + 0) = c0;
 972     *(c + 1) = c1;
 973     *(c + 2) = c2;
 974     *(c + 3) = c3;
 975 }
 976
 977 static void dsolve_2x4_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
 978 {
 979     FLOAT a0, a1, a3, c0, c1, c0_nxt1, c1_nxt1;
 980     FLOAT c0_nxt2, c1_nxt2, c0_nxt3, c1_nxt3;
 981
 982     c0 = *(c + 0);
 983     c1 = *(c + 1);
 984     c0_nxt1 = *(c + ldc);
 985     c1_nxt1 = *(c + 1 + ldc);
 986     c0_nxt2 = *(c + 2 * ldc);
 987     c1_nxt2 = *(c + 1 + 2 * ldc);
 988     c0_nxt3 = *(c + 3 * ldc);
 989     c1_nxt3 = *(c + 1 + 3 * ldc);
 990
 991     if (bk)
 992     {
 993         BLASLONG i;
 994
 995         for (i = bk; i--;)
 996         {
 997             c0 -= a[0] * b[0];
 998             c1 -= a[1] * b[0];
 999             c0_nxt1 -= a[0] * b[1];
1000             c1_nxt1 -= a[1] * b[1];
1001             c0_nxt2 -= a[0] * b[2];
1002             c1_nxt2 -= a[1] * b[2];
1003             c0_nxt3 -= a[0] * b[3];
1004             c1_nxt3 -= a[1] * b[3];
1005
1006             a += 2;
1007             b += 4;
1008         }
1009     }
1010
1011     a0 = *a;
1012     a1 = *(a + 1);
1013     a3 = *(a + 3);
1014
1015     c0 *= a0;
1016     c1 -= c0 * a1;
1017     c1 *= a3;
1018
1019     c0_nxt1 *= a0;
1020     c1_nxt1 -= c0_nxt1 * a1;
1021     c1_nxt1 *= a3;
1022
1023     c0_nxt2 *= a0;
1024     c1_nxt2 -= c0_nxt2 * a1;
1025     c1_nxt2 *= a3;
1026
1027     c0_nxt3 *= a0;
1028     c1_nxt3 -= c0_nxt3 * a1;
1029     c1_nxt3 *= a3;
1030
1031     *(b + 0) = c0;
1032     *(b + 1) = c0_nxt1;
1033     *(b + 2) = c0_nxt2;
1034     *(b + 3) = c0_nxt3;
1035     *(b + 4) = c1;
1036     *(b + 5) = c1_nxt1;
1037     *(b + 6) = c1_nxt2;
1038     *(b + 7) = c1_nxt3;
1039
1040     *(c + 0) = c0;
1041     *(c + 1) = c1;
1042     *(c + 0 + ldc) = c0_nxt1;
1043     *(c + 1 + ldc) = c1_nxt1;
1044     *(c + 0 + 2 * ldc) = c0_nxt2;
1045     *(c + 1 + 2 * ldc) = c1_nxt2;
1046     *(c + 0 + 3 * ldc) = c0_nxt3;
1047     *(c + 1 + 3 * ldc) = c1_nxt3;
1048 }
1049
1050 static void dsolve_2x2_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
1051 {
1052     FLOAT a0, a1, a3, c0, c1, c0_nxt, c1_nxt;
1053
1054     c0 = *(c + 0);
1055     c1 = *(c + 1);
1056
1057     c0_nxt = *(c + ldc);
1058     c1_nxt = *(c + 1 + ldc);
1059
1060     if (bk)
1061     {
1062         BLASLONG i;
1063
1064         for (i = bk; i--;)
1065         {
1066             c0 -= a[0] * b[0];
1067             c1 -= a[1] * b[0];
1068
1069             c0_nxt -= a[0] * b[1];
1070             c1_nxt -= a[1] * b[1];
1071
1072             a += 2;
1073             b += 2;
1074         }
1075     }
1076
1077     a0 = *a;
1078     a1 = *(a + 1);
1079     a3 = *(a + 3);
1080
1081     c0 *= a0;
1082     c1 -= c0 * a1;
1083     c1 *= a3;
1084
1085     c0_nxt *= a0;
1086     c1_nxt -= c0_nxt * a1;
1087     c1_nxt *= a3;
1088
1089     *(b + 0) = c0;
1090     *(b + 1) = c0_nxt;
1091     *(b + 2) = c1;
1092     *(b + 3) = c1_nxt;
1093
1094     *(c + 0) = c0;
1095     *(c + 1) = c1;
1096
1097     *(c + 0 + ldc) = c0_nxt;
1098     *(c + 1 + ldc) = c1_nxt;
1099 }
1100
1101 static void dsolve_2x1_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
1102 {
1103     FLOAT a0, a1, a3, c0, c1;
1104
1105     c0 = *(c + 0);
1106     c1 = *(c + 1);
1107
1108     if (bk)
1109     {
1110         BLASLONG i;
1111
1112         for (i = bk; i--;)
1113         {
1114             c0 -= a[0] * b[0];
1115             c1 -= a[1] * b[0];
1116
1117             a += 2;
1118             b += 1;
1119         }
1120     }
1121
1122     a0 = *(a + 0);
1123     a1 = *(a + 1);
1124     a3 = *(a + 3);
1125
1126     c0 *= a0;
1127     c1 -= c0 * a1;
1128     c1 *= a3;
1129
1130     *(b + 0) = c0;
1131     *(b + 1) = c1;
1132
1133     *(c + 0) = c0;
1134     *(c + 1) = c1;
1135 }
1136
1137 static void dsolve_1x4_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
1138 {
1139     FLOAT c0, c1, c2, c3;
1140
1141     c0 = *(c + 0);
1142     c1 = *(c + 1 * ldc);
1143     c2 = *(c + 2 * ldc);
1144     c3 = *(c + 3 * ldc);
1145
1146     if (bk)
1147     {
1148         BLASLONG i;
1149
1150         for (i = bk; i--;)
1151         {
1152             c0 -= a[0] * b[0];
1153             c1 -= a[0] * b[1];
1154             c2 -= a[0] * b[2];
1155             c3 -= a[0] * b[3];
1156
1157             a += 1;
1158             b += 4;
1159         }
1160     }
1161
1162     c0 *= *a;
1163     c1 *= *a;
1164     c2 *= *a;
1165     c3 *= *a;
1166
1167     *(c + 0 * ldc) = c0;
1168     *(c + 1 * ldc) = c1;
1169     *(c + 2 * ldc) = c2;
1170     *(c + 3 * ldc) = c3;
1171
1172     *(b + 0) = c0;
1173     *(b + 1) = c1;
1174     *(b + 2) = c2;
1175     *(b + 3) = c3;
1176 }
1177
1178 static void dsolve_1x2_lt_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG bk)
1179 {
1180     FLOAT c0, c1;
1181
1182     c0 = *c;
1183     c1 = *(c + ldc);
1184
1185     if (bk)
1186     {
1187         BLASLONG i;
1188
1189         for (i = bk; i--;)
1190         {
1191             c0 -= *a * b[0];
1192             c1 -= *a * b[1];
1193
1194             a += 1;
1195             b += 2;
1196         }
1197     }
1198
1199     c0 *= *a;
1200     c1 *= *a;
1201
1202     *(b + 0) = c0;
1203     *(b + 1) = c1;
1204
1205     *(c + 0) = c0;
1206     *(c + ldc) = c1;
1207 }
1208
1209 static void dgmm_dsolve_1x1_msa(FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG bk)
1210 {
1211     if (bk)
1212     {
1213         BLASLONG i;
1214
1215         for (i = bk; i--;)
1216         {
1217             *c -= *a * *b;
1218
1219             a += 1;
1220             b += 1;
1221         }
1222     }
1223
1224     *c *= *a;
1225     *b = *c;
1226 }
1227
1228 int CNAME(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT dummy1, FLOAT *a, FLOAT *b,
1229           FLOAT *c, BLASLONG ldc, BLASLONG offset)
1230 {
1231     BLASLONG i, j, kk;
1232     FLOAT *aa, *cc;
1233
1234     for (j = (n >> 2); j--;)
1235     {
1236         kk = offset;
1237         aa = a;
1238         cc = c;
1239
1240         for (i = (m >> 3); i--;)
1241         {
1242             dsolve_8x4_lt_msa(aa, b, cc, ldc, kk);
1243
1244             aa += 8 * k;
1245             cc += 8;
1246             kk += 8;
1247         }
1248
1249         if (m & 7)
1250         {
1251             if (m & 4)
1252             {
1253                 dsolve_4x4_lt_msa(aa, b, cc, ldc, kk);
1254
1255                 aa += 4 * k;
1256                 cc += 4;
1257                 kk += 4;
1258             }
1259
1260             if (m & 2)
1261             {
1262                 dsolve_2x4_lt_msa(aa, b, cc, ldc, kk);
1263
1264                 aa += 2 * k;
1265                 cc += 2;
1266                 kk += 2;
1267             }
1268
1269             if (m & 1)
1270             {
1271                 dsolve_1x4_lt_msa(aa, b, cc, ldc, kk);
1272
1273                 aa += k;
1274                 cc += 1;
1275                 kk += 1;
1276             }
1277         }
1278
1279         b += 4 * k;
1280         c += 4 * ldc;
1281     }
1282
1283     if (n & 3)
1284     {
1285         if (n & 2)
1286         {
1287             kk = offset;
1288             aa = a;
1289             cc = c;
1290
1291             for (i = (m >> 3); i--;)
1292             {
1293                 dsolve_8x2_lt_msa(aa, b, cc, ldc, kk);
1294
1295                 aa += 8 * k;
1296                 cc += 8;
1297                 kk += 8;
1298             }
1299
1300             if (m & 7)
1301             {
1302                 if (m & 4)
1303                 {
1304                     dsolve_4x2_lt_msa(aa, b, cc, ldc, kk);
1305
1306                     aa += 4 * k;
1307                     cc += 4;
1308                     kk += 4;
1309                 }
1310
1311                 if (m & 2)
1312                 {
1313                     dsolve_2x2_lt_msa(aa, b, cc, ldc, kk);
1314
1315                     aa += 2 * k;
1316                     cc += 2;
1317                     kk += 2;
1318                 }
1319
1320                 if (m & 1)
1321                 {
1322                     dsolve_1x2_lt_msa(aa, b, cc, ldc, kk);
1323
1324                     aa += k;
1325                     cc += 1;
1326                     kk += 1;
1327                 }
1328             }
1329
1330             b += 2 * k;
1331             c += 2 * ldc;
1332         }
1333
1334         if (n & 1)
1335         {
1336             kk = offset;
1337             aa = a;
1338             cc = c;
1339
1340             for (i = (m >> 3); i--;)
1341             {
1342                 dsolve_8x1_lt_msa(aa, b, cc, kk);
1343
1344                 aa += 8 * k;
1345                 cc += 8;
1346                 kk += 8;
1347             }
1348
1349             if (m & 7)
1350             {
1351                 if (m & 4)
1352                 {
1353                     dsolve_4x1_lt_msa(aa, b, cc, kk);
1354
1355                     aa += 4 * k;
1356                     cc += 4;
1357                     kk += 4;
1358                 }
1359
1360                 if (m & 2)
1361                 {
1362                     dsolve_2x1_lt_msa(aa, b, cc, kk);
1363
1364                     aa += 2 * k;
1365                     cc += 2;
1366                     kk += 2;
1367                 }
1368
1369                 if (m & 1)
1370                 {
1371                     dgmm_dsolve_1x1_msa(aa, b, cc, kk);
1372
1373                     aa += k;
1374                     cc += 1;
1375                     kk += 1;
1376                 }
1377             }
1378
1379             b += k;
1380             c += ldc;
1381         }
1382     }
1383
1384     return 0;
1385 }