`/* Implementation of the MATMUL intrinsic
- Copyright (C) 2002-2016 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 runtime library (libgfortran).
<http://www.gnu.org/licenses/>. */
#include "libgfortran.h"
-#include <stdlib.h>
#include <string.h>
#include <assert.h>'
`#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 *,
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
- = xmallocarray (size0 ((array_t *) retarray), sizeof ('rtype_name`));
- 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;
+#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
+'