]> git.ipfire.org Git - thirdparty/gcc.git/blobdiff - libgfortran/m4/matmul.m4
Update copyright years.
[thirdparty/gcc.git] / libgfortran / m4 / matmul.m4
index 6235f63e2b039c79711e1a86fd08fa620663177d..83f4ae63339b007a33fb34c08f28b0a33c6ccf1b 100644 (file)
@@ -1,44 +1,38 @@
 `/* Implementation of the MATMUL intrinsic
-   Copyright 2002, 2005, 2006 Free Software Foundation, Inc.
+   Copyright (C) 2002-2020 Free Software Foundation, Inc.
    Contributed by Paul Brook <paul@nowt.org>
 
-This file is part of the GNU Fortran 95 runtime library (libgfortran).
+This file is part of the GNU Fortran runtime library (libgfortran).
 
 Libgfortran is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public
 License as published by the Free Software Foundation; either
-version 2 of the License, or (at your option) any later version.
-
-In addition to the permissions in the GNU General Public License, the
-Free Software Foundation gives you unlimited permission to link the
-compiled version of this file into combinations with other programs,
-and to distribute those combinations without any restriction coming
-from the use of this file.  (The General Public License restrictions
-do apply in other respects; for example, they cover modification of
-the file, and distribution when not linked into a combine
-executable.)
+version 3 of the License, or (at your option) any later version.
 
 Libgfortran is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 
-You should have received a copy of the GNU General Public
-License along with libgfortran; see the file COPYING.  If not,
-write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
-Boston, MA 02110-1301, USA.  */
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
 
-#include "config.h"
-#include <stdlib.h>
+#include "libgfortran.h"
 #include <string.h>
-#include <assert.h>
-#include "libgfortran.h"'
+#include <assert.h>'
+
 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 *,
@@ -81,262 +75,144 @@ 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);
+/* Put exhaustive list of possible architectures here here, ORed together.  */
 
-/* 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->data == NULL)
-    {
-      if (GFC_DESCRIPTOR_RANK (a) == 1)
-        {
-          retarray->dim[0].lbound = 0;
-          retarray->dim[0].ubound = b->dim[1].ubound - b->dim[1].lbound;
-          retarray->dim[0].stride = 1;
-        }
-      else if (GFC_DESCRIPTOR_RANK (b) == 1)
-        {
-          retarray->dim[0].lbound = 0;
-          retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
-          retarray->dim[0].stride = 1;
-        }
-      else
-        {
-          retarray->dim[0].lbound = 0;
-          retarray->dim[0].ubound = a->dim[0].ubound - a->dim[0].lbound;
-          retarray->dim[0].stride = 1;
-
-          retarray->dim[1].lbound = 0;
-          retarray->dim[1].ubound = b->dim[1].ubound - b->dim[1].lbound;
-          retarray->dim[1].stride = retarray->dim[0].ubound+1;
-        }
+#if defined(HAVE_AVX) || defined(HAVE_AVX2) || defined(HAVE_AVX512F)
 
-      retarray->data
-       = internal_malloc_size (sizeof ('rtype_name`) * size0 ((array_t *) retarray));
-      retarray->offset = 0;
-    }
-'
-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 = retarray->dim[0].stride;
-    }
-  else
-    {
-      rxstride = retarray->dim[0].stride;
-      rystride = retarray->dim[1].stride;
-    }
-
-
-  if (GFC_DESCRIPTOR_RANK (a) == 1)
-    {
-      /* Treat it as a a row matrix A[1,count]. */
-      axstride = a->dim[0].stride;
-      aystride = 1;
-
-      xcount = 1;
-      count = a->dim[0].ubound + 1 - a->dim[0].lbound;
-    }
-  else
-    {
-      axstride = a->dim[0].stride;
-      aystride = a->dim[1].stride;
-
-      count = a->dim[1].ubound + 1 - a->dim[1].lbound;
-      xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
-    }
+#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 (count != b->dim[0].ubound + 1 - b->dim[0].lbound)
-    runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
+/* AMD-specifix funtions with AVX128 and FMA3/FMA4.  */
 
-  if (GFC_DESCRIPTOR_RANK (b) == 1)
-    {
-      /* Treat it as a column matrix B[count,1] */
-      bxstride = b->dim[0].stride;
-
-      /* 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 = b->dim[0].stride;
-      bystride = b->dim[1].stride;
-      ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
-    }
-
-  abase = a->data;
-  bbase = b->data;
-  dest = retarray->data;
+#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;
+#include <config/i386/cpuinfo.h>
+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 (__cpu_model.__cpu_vendor == VENDOR_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 (__cpu_model.__cpu_features[0] & (1 << FEATURE_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 ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
+            && (__cpu_model.__cpu_features[0] & (1 << FEATURE_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 (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
+           {
+              matmul_fn = matmul_'rtype_code`_avx;
+             goto store;
            }
-       }
-    }
+#endif  /* HAVE_AVX */
+        }
+    else if (__cpu_model.__cpu_vendor == VENDOR_AMD)
+      {
+#if defined(HAVE_AVX) && defined(HAVE_FMA3) && defined(HAVE_AVX128)
+        if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
+           && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
+         {
+            matmul_fn = matmul_'rtype_code`_avx128_fma3;
+           goto store;
+         }
+#endif
+#if defined(HAVE_AVX) && defined(HAVE_FMA4) && defined(HAVE_AVX128)
+        if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
+            && (__cpu_model.__cpu_features[0] & (1 << FEATURE_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
+'