tests/src/JIT/Performance/CodeQuality/BenchF/MatInv4/MatInv4.cs

   1 // Licensed to the .NET Foundation under one or more agreements.
   2 // The .NET Foundation licenses this file to you under the MIT license.
   3 // See the LICENSE file in the project root for more information.
   4
   5 using Microsoft.Xunit.Performance;
   6 using System;
   7 using System.Runtime.CompilerServices;
   8 using Xunit;
   9
  10 [assembly: OptimizeForBenchmarks]
  11
  12 namespace Benchstone.BenchF
  13 {
  14 public static class MatInv4
  15 {
  16 #if DEBUG
  17     public const int Iterations = 1;
  18 #else
  19     public const int Iterations = 60;
  20 #endif
  21
  22     private static float s_det;
  23
  24     private struct X
  25     {
  26         public float[] A;
  27         public X(int size)
  28         {
  29             A = new float[size];
  30         }
  31     }
  32
  33     [MethodImpl(MethodImplOptions.NoInlining)]
  34     private static bool Bench()
  35     {
  36         X a = new X(Iterations * Iterations);
  37         float[] b = new float[Iterations * Iterations];
  38         float[] c = new float[Iterations * Iterations];
  39         float[] d = new float[Iterations * Iterations];
  40         float[] l1 = new float[Iterations];
  41         float[] l2 = new float[Iterations];
  42
  43         int i, k, n, nsq;
  44
  45         n = Iterations;
  46         nsq = n * n;
  47         for (i = 0; i < n; ++i)
  48         {
  49             for (k = 0; k < n; ++k)
  50             {
  51                 if (i == k)
  52                 {
  53                     a.A[i * n + k] = 40.0F;
  54                 }
  55                 else
  56                 {
  57                     a.A[i * n + k] = 0.0F;
  58                 }
  59             }
  60         }
  61
  62         for (i = 0; i < n; ++i)
  63         {
  64             for (k = i; k < nsq; k += n)
  65             {
  66                 b[k] = a.A[k];
  67             }
  68         }
  69
  70         /*** second(&t1); ***/
  71
  72         MinV1(b, ref n, out s_det, l1, l2);
  73
  74         if (s_det == 0.0F)
  75         {
  76             goto L990;
  77         }
  78
  79         /*** second(&tx); ***/
  80
  81         MProd(b, a.A, c, ref n);
  82         for (k = 1; k <= nsq; ++k)
  83         {
  84             b[k - 1] = a.A[k - 1];
  85         }
  86
  87         /*** second(&tx); ***/
  88
  89         MinV2(b, ref n, out s_det, l1, l2);
  90
  91         if (s_det == 0.0F)
  92         {
  93             goto L990;
  94         }
  95
  96         /*** second(&ty); ***/
  97
  98         MProd(b, a.A, d, ref n);
  99         CompM(c, d, ref n);
 100
 101         /*** second(&t2); ***/
 102
 103         return true;
 104
 105     L990:
 106         {
 107         }
 108
 109         return true;
 110     }
 111
 112     private static void MinV1(float[] a, ref int n, out float d, float[] l, float[] m)
 113     {
 114         float biga, hold;
 115         int i, j, k, ij, ik, ji, jk, nk, ki, kj, kk, iz, jp, jq, jr;
 116
 117         d = 1.0F;
 118         ji = 0;
 119         hold = 0.0F;
 120         nk = -n;
 121         for (k = 1; k <= n; ++k)
 122         {
 123             nk = nk + n;
 124             l[k - 1] = k;
 125             m[k - 1] = k;
 126             kk = nk + k;
 127             biga = a[kk - 1];
 128             for (j = k; j <= n; ++j)
 129             {
 130                 // j <= n, so iz <= n^2 - n
 131                 iz = n * (j - 1);
 132                 for (i = k; i <= n; ++i)
 133                 {
 134                     // iz <= n^2 - n and i <= n, so ij <= n^2
 135                     ij = iz + i;
 136                     if (System.Math.Abs(biga) >= System.Math.Abs(a[ij - 1]))
 137                     {
 138                         continue;
 139                     }
 140                     // accessing up to n^2 - 1
 141                     biga = a[ij - 1];
 142                     l[k - 1] = i;
 143                     m[k - 1] = j;
 144                 }
 145             }
 146
 147             j = (int)l[k - 1];
 148
 149             if (j <= k)
 150             {
 151                 goto L35;
 152             }
 153
 154             // -n < ki <= 0
 155             ki = k - n;
 156             for (i = 1; i <= n; ++i)
 157             {
 158                 // i <= n, ki <= n + n + ... + n (n times) i.e. k <= n * n (when ki = 0 initially)
 159                 ki = ki + n;
 160                 // Accessing upto n^2 -1
 161                 hold = -a[ki - 1];
 162                 // ji <= n^2 - n + n (for ki = 0 initially when k = n and 0 < j <= n)
 163                 // Therefore ji <= n^2
 164                 ji = ki - k + j;
 165                 a[ki - 1] = a[ji - 1];
 166                 a[ji - 1] = hold;
 167             }
 168         L35:
 169             i = (int)m[k - 1];
 170             if (i <= k)
 171             {
 172                 goto L45;
 173             }
 174
 175             // 0 <= jp <= n^2 - n
 176             jp = n * (i - 1);
 177             for (j = 1; j <= n; ++j)
 178             {
 179                 // 0 < nk <= n * (n-1)
 180                 // jk <= n^2 - n + n
 181                 // jk <= n^2
 182                 jk = nk + j;
 183                 // jp <= n^2 - n
 184                 // ji <= n^2 - n + n or ji <= n^2 (since 0 < j <= n)
 185                 ji = jp + j;
 186                 hold = -a[jk - 1];
 187                 a[jk - 1] = a[ji - 1];
 188                 a[ji - 1] = hold;
 189             }
 190         L45:
 191             if (biga != 0.0F)
 192             {
 193                 goto L48;
 194             }
 195             d = 0.0F;
 196             return;
 197
 198         L48:
 199             for (i = 1; i <= n; ++i)
 200             {
 201                 if (i == k)
 202                 {
 203                     break;
 204                 }
 205                 // 0 < nk <= n * (n-1)
 206                 // 0 < ik <= n^2
 207                 ik = nk + i;
 208                 a[ik - 1] = a[ik - 1] / (-biga);
 209             }
 210
 211             for (i = 1; i <= n; ++i)
 212             {
 213                 if (i == k)
 214                 {
 215                     continue;
 216                 }
 217                 // 0 < nk <= n * (n-1)
 218                 // 0 < ik <= n^2
 219                 ik = nk + i;
 220                 hold = a[ik - 1];
 221                 // -n < ij <= 0
 222                 ij = i - n;
 223                 for (j = 1; j <= n; ++j)
 224                 {
 225                     // i <= n, ij <= n + n + ... + n (n times) or ij <= n * n
 226                     ij = ij + n;
 227                     if (j == k)
 228                     {
 229                         continue;
 230                     }
 231                     // if i=1, kj = (1 + (n-1) * n) - 1 + n ==> ij = n^2
 232                     // if i=n, kj = (n * n) - n + n ==> ij = n ^2
 233                     // So j <= n^2
 234                     kj = ij - i + k;
 235                     a[ij - 1] = hold * a[kj - 1] + a[ij - 1];
 236                 }
 237             }
 238             kj = k - n;
 239             for (j = 1; j <= n; ++j)
 240             {
 241                 // k <= n, kj <= n + n + ... + n (n times) or kj <= n * n
 242                 kj = kj + n;
 243                 if (j == k)
 244                 {
 245                     continue;
 246                 }
 247                 // Accessing upto n^2 - 1
 248                 a[kj - 1] = a[kj - 1] / biga;
 249             }
 250             d = d * biga;
 251             a[kk - 1] = 1.0F / biga;
 252         }
 253         k = n;
 254     L100:
 255         k = k - 1;
 256         if (k < 1)
 257         {
 258             return;
 259         }
 260         i = (int)l[k - 1];
 261         if (i <= k)
 262         {
 263             goto L120;
 264         }
 265
 266         // 0 <= jq <= n^2 - n
 267         // 0 <= jr <= n^2 - n
 268         jq = n * (k - 1);
 269         jr = n * (i - 1);
 270         for (j = 1; j <= n; ++j)
 271         {
 272             // jk <= n^2 - n + n
 273             // jk <= n^2
 274             jk = jq + j;
 275             hold = a[jk - 1];
 276             // ji <= n^2 - n + n
 277             // ji <= n^2
 278             ji = jr + j;
 279             a[jk - 1] = -a[ji - 1];
 280             a[ji - 1] = hold;
 281         }
 282     L120:
 283         j = (int)m[k - 1];
 284         if (j <= k)
 285         {
 286             goto L100;
 287         }
 288         // 0 <= jr <= n^2 - n
 289         ki = k - n;
 290         for (i = 1; i <= n; ++i)
 291         {
 292             // ki <= n + n + ... + n (n times) or ki <= n * n
 293             ki = ki + n;
 294             hold = a[ki - 1];
 295             // if i=1, ji = (1 + (n-1) * n) - 1 + n ==> ij = n^2
 296             // if i=n, ji = (n * n) - n + n ==> ij = n ^2
 297             // Therefore ji <= n^2
 298             ji = ki - k + j;
 299             a[ki - 1] = -a[ji - 1];
 300         }
 301         a[ji - 1] = hold;
 302         goto L100;
 303     }
 304
 305     private static void MinV2(float[] a, ref int n, out float d, float[] l, float[] m)
 306     {
 307         float biga, hold;
 308         int i, j, k;
 309
 310         d = 1.0F;
 311         for (k = 1; k <= n; ++k)
 312         {
 313             l[k - 1] = k;
 314             m[k - 1] = k;
 315             biga = a[(k - 1) * n + (k - 1)];
 316             for (j = k; j <= n; ++j)
 317             {
 318                 for (i = k; i <= n; ++i)
 319                 {
 320                     // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 321                     if (System.Math.Abs(biga) >= System.Math.Abs(a[(i - 1) * n + (j - 1)]))
 322                     {
 323                         continue;
 324                     }
 325                     biga = a[(i - 1) * n + (j - 1)];
 326                     l[k - 1] = i;
 327                     m[k - 1] = j;
 328                 }
 329             }
 330             j = (int)l[k - 1];
 331             if (l[k - 1] <= k)
 332             {
 333                 goto L200;
 334             }
 335             for (i = 1; i <= n; ++i)
 336             {
 337                 // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 338                 hold = -a[(k - 1) * n + (i - 1)];
 339                 a[(k - 1) * n + (i - 1)] = a[(j - 1) * n + (i - 1)];
 340                 a[(j - 1) * n + (i - 1)] = hold;
 341             }
 342         L200:
 343             i = (int)m[k - 1];
 344             if (m[k - 1] <= k)
 345             {
 346                 goto L250;
 347             }
 348             for (j = 1; j <= n; ++j)
 349             {
 350                 // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 351                 hold = -a[(j - 1) * n + (k - 1)];
 352                 a[(j - 1) * n + (k - 1)] = a[(j - 1) * n + (i - 1)];
 353                 a[(j - 1) * n + (i - 1)] = hold;
 354             }
 355         L250:
 356             if (biga != 0.0F)
 357             {
 358                 goto L300;
 359             }
 360             d = 0.0F;
 361             return;
 362
 363         L300:
 364             for (i = 1; i <= n; ++i)
 365             {
 366                 if (i != k)
 367                 {
 368                     // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 369                     a[(i - 1) * n + (k - 1)] = a[(i - 1) * n + (k - 1)] / (-biga);
 370                 }
 371             }
 372             for (i = 1; i <= n; ++i)
 373             {
 374                 if (i == k)
 375                 {
 376                     continue;
 377                 }
 378                 for (j = 1; j <= n; ++j)
 379                 {
 380                     if (j != k)
 381                     {
 382                         // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 383                         a[(i - 1) * n + (j - 1)] = a[(i - 1) * n + (k - 1)] * a[(k - 1) * n + (j - 1)] + a[(i - 1) * n + (j - 1)];
 384                     }
 385                 }
 386             }
 387             for (j = 1; j < n; ++j)
 388             {
 389                 if (j != k)
 390                 {
 391                     // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 392                     a[(k - 1) * n + (j - 1)] = a[(k - 1) * n + (j - 1)] / biga;
 393                 }
 394             }
 395             d = d * biga;
 396             a[(k - 1) * n + (k - 1)] = 1.0F / biga;
 397         }
 398         k = n;
 399     L400:
 400         k = k - 1;
 401         if (k < 1)
 402         {
 403             return;
 404         }
 405         i = (int)l[k - 1];
 406         if (i <= k)
 407         {
 408             goto L450;
 409         }
 410         for (j = 1; j <= n; ++j)
 411         {
 412             // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 413             hold = a[(j - 1) * n + (k - 1)];
 414             a[(j - 1) * n + (k - 1)] = -a[(j - 1) * n + (i - 1)];
 415             a[(j - 1) * n + (i - 1)] = hold;
 416         }
 417     L450:
 418         j = (int)m[k - 1];
 419         if (j <= k)
 420         {
 421             goto L400;
 422         }
 423         for (i = 1; i <= n; ++i)
 424         {
 425             // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 426             hold = a[(k - 1) * n + (i - 1)];
 427             a[(k - 1) * n + (i - 1)] = -a[(j - 1) * n + (i - 1)];
 428             a[(j - 1) * n + (i - 1)] = hold;
 429         }
 430         goto L400;
 431     }
 432
 433     private static void MProd(float[] a, float[] b, float[] c, ref int n)
 434     {
 435         int i, j, k;
 436
 437         for (i = 1; i <= n; ++i)
 438         {
 439             for (j = 1; j <= n; ++j)
 440             {
 441                 // Accessing upto n^2 - n + n - 1 ==> n^2 - 1
 442                 c[(i - 1) * n + (j - 1)] = 0.0F;
 443                 for (k = 1; k <= n; ++k)
 444                 {
 445                     c[(i - 1) * n + (j - 1)] = c[(i - 1) * n + (j - 1)] + a[(i - 1) * n + (k - 1)] * b[(k - 1) * n + (j - 1)];
 446                 }
 447             }
 448         }
 449         return;
 450     }
 451
 452     private static void CompM(float[] a, float[] b, ref int n)
 453     {
 454         int i, j;
 455         float x, sum = 0.0F;
 456
 457         //(starting compare.)
 458         for (i = 1; i <= n; ++i)
 459         {
 460             for (j = 1; j <= n; ++j)
 461             {
 462                 x = 0.0F;
 463                 if (i == j)
 464                 {
 465                     x = 1.0F;
 466                 }
 467                 sum = sum + System.Math.Abs(System.Math.Abs(a[(i - 1) * n + (j - 1)]) - x);
 468             }
 469         }
 470         return;
 471     }
 472
 473     [Benchmark]
 474     public static void Test()
 475     {
 476         foreach (var iteration in Benchmark.Iterations)
 477         {
 478             using (iteration.StartMeasurement())
 479             {
 480                 Bench();
 481             }
 482         }
 483     }
 484
 485     private static bool TestBase()
 486     {
 487         bool result = Bench();
 488         return result;
 489     }
 490
 491     public static int Main()
 492     {
 493         bool result = TestBase();
 494         return (result ? 100 : -1);
 495     }
 496 }
 497 }