`/* 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 *,
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
+'