`/* Implementation of the MATMUL intrinsic
- Copyright 2002, 2005 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 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. */
-
-#include "config.h"
-#include <stdlib.h>
-#include <string.h>
-#include <assert.h>
-#include "libgfortran.h"'
-include(iparm.m4)dnl
-
-`#if defined (HAVE_'rtype_name`)'
+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.
-/* This is a C version of the following fortran pseudo-code. The key
- point is the loop order -- we access all arrays column-first, which
- improves the performance enough to boost galgel spec score by 50%.
+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/>. */
- DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
- C = 0
- DO J=1,N
- DO K=1,COUNT
- DO I=1,M
- C(I,J) = C(I,J)+A(I,K)*B(K,J)
-*/
-
-extern void matmul_`'rtype_code (rtype * retarray, rtype * a, rtype * b);
-export_proto(matmul_`'rtype_code);
+#include "libgfortran.h"
+#include <string.h>
+#include <assert.h>'
-void
-matmul_`'rtype_code (rtype * retarray, rtype * a, rtype * b)
-{
- rtype_name *abase;
- rtype_name *bbase;
- rtype_name *dest;
+include(iparm.m4)dnl
- index_type rxstride, rystride, axstride, aystride, bxstride, bystride;
- index_type x, y, n, count, xcount, ycount;
+`#if defined (HAVE_'rtype_name`)
- assert (GFC_DESCRIPTOR_RANK (a) == 2
- || GFC_DESCRIPTOR_RANK (b) == 2);
+/* Prototype for the BLAS ?gemm subroutine, a pointer to which can be
+ passed to us by the front-end, in which case we call it for large
+ matrices. */
-/* C[xcount,ycount] = A[xcount, count] * B[count,ycount]
+typedef void (*blas_call)(const char *, const char *, const int *, const int *,
+ const int *, const 'rtype_name` *, const 'rtype_name` *,
+ const int *, const 'rtype_name` *, const int *,
+ const 'rtype_name` *, 'rtype_name` *, const int *,
+ int, int);
- Either A or B (but not both) can be rank 1:
+/* The order of loops is different in the case of plain matrix
+ multiplication C=MATMUL(A,B), and in the frequent special case where
+ the argument A is the temporary result of a TRANSPOSE intrinsic:
+ C=MATMUL(TRANSPOSE(A),B). Transposed temporaries are detected by
+ looking at their strides.
- o One-dimensional argument A is implicitly treated as a row matrix
- dimensioned [1,count], so xcount=1.
+ The equivalent Fortran pseudo-code is:
- o One-dimensional argument B is implicitly treated as a column matrix
- dimensioned [count, 1], so ycount=1.
- */
+ DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
+ IF (.NOT.IS_TRANSPOSED(A)) THEN
+ C = 0
+ DO J=1,N
+ DO K=1,COUNT
+ DO I=1,M
+ C(I,J) = C(I,J)+A(I,K)*B(K,J)
+ ELSE
+ DO J=1,N
+ DO I=1,M
+ S = 0
+ DO K=1,COUNT
+ S = S+A(I,K)*B(K,J)
+ C(I,J) = S
+ ENDIF
+*/
- 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 try_blas is set to a nonzero value, then the matmul function will
+ see if there is a way to perform the matrix multiplication by a call
+ to the BLAS gemm function. */
+
+extern 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);
+export_proto(matmul_'rtype_code`);
+
+/* Put exhaustive list of possible architectures here here, ORed together. */
+
+#if defined(HAVE_AVX) || defined(HAVE_AVX2) || defined(HAVE_AVX512F)
+
+#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 */
+
+/* AMD-specifix funtions with AVX128 and FMA3/FMA4. */
+
+#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
- retarray->data
- = internal_malloc_size (sizeof (rtype_name) * size0 ((array_t *) retarray));
- retarray->offset = 0;
- }
+#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
- abase = a->data;
- bbase = b->data;
- dest = retarray->data;
+/* 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
- if (retarray->dim[0].stride == 0)
- retarray->dim[0].stride = 1;
- if (a->dim[0].stride == 0)
- a->dim[0].stride = 1;
- if (b->dim[0].stride == 0)
- b->dim[0].stride = 1;
+`/* Compiling main function, with selection code for the processor. */
-sinclude(`matmul_asm_'rtype_code`.m4')dnl
+/* Currently, this is i386 only. Adjust for other architectures. */
- 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;
- }
+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);
+ 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 (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
+ matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+ if (matmul_fn == NULL)
{
- axstride = a->dim[0].stride;
- aystride = a->dim[1].stride;
+ matmul_fn = matmul_'rtype_code`_vanilla;
+ if (__builtin_cpu_is ("intel"))
+ {
+ /* Run down the available processors in order of preference. */
+#ifdef HAVE_AVX512F
+ if (__builtin_cpu_supports ("avx512f"))
+ {
+ matmul_fn = matmul_'rtype_code`_avx512f;
+ goto store;
+ }
- count = a->dim[1].ubound + 1 - a->dim[1].lbound;
- xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
- }
+#endif /* HAVE_AVX512F */
- assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
+#ifdef HAVE_AVX2
+ if (__builtin_cpu_supports ("avx2")
+ && __builtin_cpu_supports ("fma"))
+ {
+ matmul_fn = matmul_'rtype_code`_avx2;
+ goto store;
+ }
- 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;
- }
+#endif
- abase = a->data;
- bbase = b->data;
- dest = retarray->data;
+#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
- if (rxstride == 1 && axstride == 1 && bxstride == 1)
- {
- rtype_name *bbase_y;
- rtype_name *dest_y;
- rtype_name *abase_n;
- rtype_name bbase_yn;
-
- if (rystride == ycount)
- memset (dest, 0, (sizeof (rtype_name) * size0((array_t *) retarray)));
- else
- {
- for (y = 0; y < ycount; y++)
- for (x = 0; x < xcount; x++)
- dest[x + y*rystride] = (rtype_name)0;
- }
+ }
+ store:
+ __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
+ }
- for (y = 0; y < ycount; y++)
- {
- bbase_y = bbase + y*bystride;
- dest_y = dest + y*rystride;
- for (n = 0; n < count; n++)
- {
- abase_n = abase + n*aystride;
- bbase_yn = bbase_y[n];
- for (x = 0; x < xcount; x++)
- {
- dest_y[x] += abase_n[x] * bbase_yn;
- }
- }
- }
- }
- else
- {
- 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];
- }
+ (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
}
+#else /* Just the vanilla function. */
+
+'define(`matmul_name',`matmul_'rtype_code)dnl
+define(`target_attribute',`')dnl
+include(matmul_internal.m4)dnl
+`#endif
#endif
+'
projects / thirdparty / gcc.git / blobdiff