]> git.ipfire.org Git - thirdparty/gcc.git/blobdiff - libgfortran/m4/matmul.m4
Update copyright years.
[thirdparty/gcc.git] / libgfortran / m4 / matmul.m4
index 6a58cd91adbe93dd0a8bf58b4d3f91642e7408e1..7fc1f5fa75fb1ce4f4607e1d2e2e2310dfe37f28 100644 (file)
@@ -1,5 +1,5 @@
 `/* Implementation of the MATMUL intrinsic
-   Copyright 2002, 2005, 2006, 2007, 2009, 2012 Free Software Foundation, Inc.
+   Copyright (C) 2002-2024 Free Software Foundation, Inc.
    Contributed by Paul Brook <paul@nowt.org>
 
 This file is part of the GNU Fortran runtime library (libgfortran).
@@ -24,7 +24,6 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 <http://www.gnu.org/licenses/>.  */
 
 #include "libgfortran.h"
-#include <stdlib.h>
 #include <string.h>
 #include <assert.h>'
 
@@ -33,7 +32,7 @@ include(iparm.m4)dnl
 `#if defined (HAVE_'rtype_name`)
 
 /* Prototype for the BLAS ?gemm subroutine, a pointer to which can be
-   passed to us by the front-end, in which case we''`ll call it for large
+   passed to us by the front-end, in which case we call it for large
    matrices.  */
 
 typedef void (*blas_call)(const char *, const char *, const int *, const int *,
@@ -76,303 +75,143 @@ extern void matmul_'rtype_code` ('rtype` * const restrict retarray,
        int blas_limit, blas_call gemm);
 export_proto(matmul_'rtype_code`);
 
-void
-matmul_'rtype_code` ('rtype` * const restrict retarray, 
-       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
-       int blas_limit, blas_call gemm)
-{
-  const 'rtype_name` * restrict abase;
-  const 'rtype_name` * restrict bbase;
-  'rtype_name` * restrict dest;
-
-  index_type rxstride, rystride, axstride, aystride, bxstride, bystride;
-  index_type x, y, n, count, xcount, ycount;
-
-  assert (GFC_DESCRIPTOR_RANK (a) == 2
-          || GFC_DESCRIPTOR_RANK (b) == 2);
-
-/* C[xcount,ycount] = A[xcount, count] * B[count,ycount]
-
-   Either A or B (but not both) can be rank 1:
-
-   o One-dimensional argument A is implicitly treated as a row matrix
-     dimensioned [1,count], so xcount=1.
-
-   o One-dimensional argument B is implicitly treated as a column matrix
-     dimensioned [count, 1], so ycount=1.
-  */
-
-  if (retarray->base_addr == NULL)
-    {
-      if (GFC_DESCRIPTOR_RANK (a) == 1)
-        {
-         GFC_DIMENSION_SET(retarray->dim[0], 0,
-                           GFC_DESCRIPTOR_EXTENT(b,1) - 1, 1);
-        }
-      else if (GFC_DESCRIPTOR_RANK (b) == 1)
-        {
-         GFC_DIMENSION_SET(retarray->dim[0], 0,
-                           GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
-        }
-      else
-        {
-         GFC_DIMENSION_SET(retarray->dim[0], 0,
-                           GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
-
-          GFC_DIMENSION_SET(retarray->dim[1], 0,
-                           GFC_DESCRIPTOR_EXTENT(b,1) - 1,
-                           GFC_DESCRIPTOR_EXTENT(retarray,0));
-        }
-
-      retarray->base_addr
-       = xmalloc (sizeof ('rtype_name`) * size0 ((array_t *) retarray));
-      retarray->offset = 0;
-    }
-    else if (unlikely (compile_options.bounds_check))
-      {
-       index_type ret_extent, arg_extent;
-
-       if (GFC_DESCRIPTOR_RANK (a) == 1)
-         {
-           arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
-           ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
-           if (arg_extent != ret_extent)
-             runtime_error ("Incorrect extent in return array in"
-                            " MATMUL intrinsic: is %ld, should be %ld",
-                            (long int) ret_extent, (long int) arg_extent);
-         }
-       else if (GFC_DESCRIPTOR_RANK (b) == 1)
-         {
-           arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
-           ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
-           if (arg_extent != ret_extent)
-             runtime_error ("Incorrect extent in return array in"
-                            " MATMUL intrinsic: is %ld, should be %ld",
-                            (long int) ret_extent, (long int) arg_extent);         
-         }
-       else
-         {
-           arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
-           ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
-           if (arg_extent != ret_extent)
-             runtime_error ("Incorrect extent in return array in"
-                            " MATMUL intrinsic for dimension 1:"
-                            " is %ld, should be %ld",
-                            (long int) ret_extent, (long int) arg_extent);
-
-           arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
-           ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
-           if (arg_extent != ret_extent)
-             runtime_error ("Incorrect extent in return array in"
-                            " MATMUL intrinsic for dimension 2:"
-                            " is %ld, should be %ld",
-                            (long int) ret_extent, (long int) arg_extent);
-         }
-      }
-'
-sinclude(`matmul_asm_'rtype_code`.m4')dnl
-`
-  if (GFC_DESCRIPTOR_RANK (retarray) == 1)
-    {
-      /* One-dimensional result may be addressed in the code below
-        either as a row or a column matrix. We want both cases to
-        work. */
-      rxstride = rystride = GFC_DESCRIPTOR_STRIDE(retarray,0);
-    }
-  else
-    {
-      rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
-      rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
-    }
+/* Put exhaustive list of possible architectures here here, ORed together.  */
 
+#if defined(HAVE_AVX) || defined(HAVE_AVX2) || defined(HAVE_AVX512F)
 
-  if (GFC_DESCRIPTOR_RANK (a) == 1)
-    {
-      /* Treat it as a a row matrix A[1,count]. */
-      axstride = GFC_DESCRIPTOR_STRIDE(a,0);
-      aystride = 1;
-
-      xcount = 1;
-      count = GFC_DESCRIPTOR_EXTENT(a,0);
-    }
-  else
-    {
-      axstride = GFC_DESCRIPTOR_STRIDE(a,0);
-      aystride = GFC_DESCRIPTOR_STRIDE(a,1);
-
-      count = GFC_DESCRIPTOR_EXTENT(a,1);
-      xcount = GFC_DESCRIPTOR_EXTENT(a,0);
-    }
-
-  if (count != GFC_DESCRIPTOR_EXTENT(b,0))
-    {
-      if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
-       runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
-    }
+#ifdef HAVE_AVX
+'define(`matmul_name',`matmul_'rtype_code`_avx')dnl
+`static void
+'matmul_name` ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm) __attribute__((__target__("avx")));
+static' include(matmul_internal.m4)dnl
+`#endif /* HAVE_AVX */
+
+#ifdef HAVE_AVX2
+'define(`matmul_name',`matmul_'rtype_code`_avx2')dnl
+`static void
+'matmul_name` ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm) __attribute__((__target__("avx2,fma")));
+static' include(matmul_internal.m4)dnl
+`#endif /* HAVE_AVX2 */
+
+#ifdef HAVE_AVX512F
+'define(`matmul_name',`matmul_'rtype_code`_avx512f')dnl
+`static void
+'matmul_name` ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm) __attribute__((__target__("avx512f")));
+static' include(matmul_internal.m4)dnl
+`#endif  /* HAVE_AVX512F */
 
-  if (GFC_DESCRIPTOR_RANK (b) == 1)
-    {
-      /* Treat it as a column matrix B[count,1] */
-      bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
-
-      /* bystride should never be used for 1-dimensional b.
-        in case it is we want it to cause a segfault, rather than
-        an incorrect result. */
-      bystride = 0xDEADBEEF;
-      ycount = 1;
-    }
-  else
-    {
-      bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
-      bystride = GFC_DESCRIPTOR_STRIDE(b,1);
-      ycount = GFC_DESCRIPTOR_EXTENT(b,1);
-    }
+/* AMD-specifix funtions with AVX128 and FMA3/FMA4.  */
 
-  abase = a->base_addr;
-  bbase = b->base_addr;
-  dest = retarray->base_addr;
+#if defined(HAVE_AVX) && defined(HAVE_FMA3) && defined(HAVE_AVX128)
+'define(`matmul_name',`matmul_'rtype_code`_avx128_fma3')dnl
+`void
+'matmul_name` ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm) __attribute__((__target__("avx,fma")));
+internal_proto('matmul_name`);
+#endif
+
+#if defined(HAVE_AVX) && defined(HAVE_FMA4) && defined(HAVE_AVX128)
+'define(`matmul_name',`matmul_'rtype_code`_avx128_fma4')dnl
+`void
+'matmul_name` ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm) __attribute__((__target__("avx,fma4")));
+internal_proto('matmul_name`);
+#endif
 
+/* Function to fall back to if there is no special processor-specific version.  */
+'define(`matmul_name',`matmul_'rtype_code`_vanilla')dnl
+`static' include(matmul_internal.m4)dnl
 
-  /* Now that everything is set up, we''`re performing the multiplication
-     itself.  */
+`/* Compiling main function, with selection code for the processor.  */
 
-#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x)))
+/* Currently, this is i386 only.  Adjust for other architectures.  */
 
-  if (try_blas && rxstride == 1 && (axstride == 1 || aystride == 1)
-      && (bxstride == 1 || bystride == 1)
-      && (((float) xcount) * ((float) ycount) * ((float) count)
-          > POW3(blas_limit)))
-  {
-    const int m = xcount, n = ycount, k = count, ldc = rystride;
-    const 'rtype_name` one = 1, zero = 0;
-    const int lda = (axstride == 1) ? aystride : axstride,
-              ldb = (bxstride == 1) ? bystride : bxstride;
+void matmul_'rtype_code` ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm)
+{
+  static void (*matmul_p) ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm);
 
-    if (lda > 0 && ldb > 0 && ldc > 0 && m > 1 && n > 1 && k > 1)
-      {
-        assert (gemm != NULL);
-        gemm (axstride == 1 ? "N" : "T", bxstride == 1 ? "N" : "T", &m, &n, &k,
-              &one, abase, &lda, bbase, &ldb, &zero, dest, &ldc, 1, 1);
-        return;
-      }
-  }
+  void (*matmul_fn) ('rtype` * const restrict retarray, 
+       'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+       int blas_limit, blas_call gemm);
 
-  if (rxstride == 1 && axstride == 1 && bxstride == 1)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      const 'rtype_name` * restrict bbase_y;
-      'rtype_name` * restrict dest_y;
-      const 'rtype_name` * restrict abase_n;
-      'rtype_name` bbase_yn;
-
-      if (rystride == xcount)
-       memset (dest, 0, (sizeof ('rtype_name`) * xcount * ycount));
-      else
-       {
-         for (y = 0; y < ycount; y++)
-           for (x = 0; x < xcount; x++)
-             dest[x + y*rystride] = ('rtype_name`)0;
-       }
-
-      for (y = 0; y < ycount; y++)
+      matmul_fn = matmul_'rtype_code`_vanilla;
+      if (__builtin_cpu_is ("intel"))
        {
-         bbase_y = bbase + y*bystride;
-         dest_y = dest + y*rystride;
-         for (n = 0; n < count; n++)
+          /* Run down the available processors in order of preference.  */
+#ifdef HAVE_AVX512F
+         if (__builtin_cpu_supports ("avx512f"))
            {
-             abase_n = abase + n*aystride;
-             bbase_yn = bbase_y[n];
-             for (x = 0; x < xcount; x++)
-               {
-                 dest_y[x] += abase_n[x] * bbase_yn;
-               }
+             matmul_fn = matmul_'rtype_code`_avx512f;
+             goto store;
            }
-       }
-    }
-  else if (rxstride == 1 && aystride == 1 && bxstride == 1)
-    {
-      if (GFC_DESCRIPTOR_RANK (a) != 1)
-       {
-         const 'rtype_name` *restrict abase_x;
-         const 'rtype_name` *restrict bbase_y;
-         'rtype_name` *restrict dest_y;
-         'rtype_name` s;
 
-         for (y = 0; y < ycount; y++)
-           {
-             bbase_y = &bbase[y*bystride];
-             dest_y = &dest[y*rystride];
-             for (x = 0; x < xcount; x++)
-               {
-                 abase_x = &abase[x*axstride];
-                 s = ('rtype_name`) 0;
-                 for (n = 0; n < count; n++)
-                   s += abase_x[n] * bbase_y[n];
-                 dest_y[x] = s;
-               }
-           }
-       }
-      else
-       {
-         const 'rtype_name` *restrict bbase_y;
-         'rtype_name` s;
+#endif  /* HAVE_AVX512F */
 
-         for (y = 0; y < ycount; y++)
+#ifdef HAVE_AVX2
+         if (__builtin_cpu_supports ("avx2")
+             && __builtin_cpu_supports ("fma"))
            {
-             bbase_y = &bbase[y*bystride];
-             s = ('rtype_name`) 0;
-             for (n = 0; n < count; n++)
-               s += abase[n*axstride] * bbase_y[n];
-             dest[y*rystride] = s;
+             matmul_fn = matmul_'rtype_code`_avx2;
+             goto store;
            }
-       }
-    }
-  else if (axstride < aystride)
-    {
-      for (y = 0; y < ycount; y++)
-       for (x = 0; x < xcount; x++)
-         dest[x*rxstride + y*rystride] = ('rtype_name`)0;
-
-      for (y = 0; y < ycount; y++)
-       for (n = 0; n < count; n++)
-         for (x = 0; x < xcount; x++)
-           /* dest[x,y] += a[x,n] * b[n,y] */
-           dest[x*rxstride + y*rystride] += abase[x*axstride + n*aystride] * bbase[n*bxstride + y*bystride];
-    }
-  else if (GFC_DESCRIPTOR_RANK (a) == 1)
-    {
-      const 'rtype_name` *restrict bbase_y;
-      'rtype_name` s;
 
-      for (y = 0; y < ycount; y++)
-       {
-         bbase_y = &bbase[y*bystride];
-         s = ('rtype_name`) 0;
-         for (n = 0; n < count; n++)
-           s += abase[n*axstride] * bbase_y[n*bxstride];
-         dest[y*rxstride] = s;
-       }
-    }
-  else
-    {
-      const 'rtype_name` *restrict abase_x;
-      const 'rtype_name` *restrict bbase_y;
-      'rtype_name` *restrict dest_y;
-      'rtype_name` s;
+#endif
 
-      for (y = 0; y < ycount; y++)
-       {
-         bbase_y = &bbase[y*bystride];
-         dest_y = &dest[y*rystride];
-         for (x = 0; x < xcount; x++)
-           {
-             abase_x = &abase[x*axstride];
-             s = ('rtype_name`) 0;
-             for (n = 0; n < count; n++)
-               s += abase_x[n*aystride] * bbase_y[n*bxstride];
-             dest_y[x*rxstride] = s;
+#ifdef HAVE_AVX
+         if (__builtin_cpu_supports ("avx"))
+           {
+              matmul_fn = matmul_'rtype_code`_avx;
+             goto store;
            }
-       }
-    }
+#endif  /* HAVE_AVX */
+        }
+    else if (__builtin_cpu_is ("amd"))
+      {
+#if defined(HAVE_AVX) && defined(HAVE_FMA3) && defined(HAVE_AVX128)
+       if (__builtin_cpu_supports ("avx")
+           && __builtin_cpu_supports ("fma"))
+         {
+            matmul_fn = matmul_'rtype_code`_avx128_fma3;
+           goto store;
+         }
+#endif
+#if defined(HAVE_AVX) && defined(HAVE_FMA4) && defined(HAVE_AVX128)
+       if (__builtin_cpu_supports ("avx")
+           && __builtin_cpu_supports ("fma4"))
+         {
+            matmul_fn = matmul_'rtype_code`_avx128_fma4;
+           goto store;
+         }
+#endif
+
+      }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
+   }
+
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
-#endif'
+#else  /* Just the vanilla function.  */
+
+'define(`matmul_name',`matmul_'rtype_code)dnl
+define(`target_attribute',`')dnl
+include(matmul_internal.m4)dnl
+`#endif
+#endif
+'