diff options
Diffstat (limited to 'libgfortran/generated')
-rw-r--r-- | libgfortran/generated/matmul_c10.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_c16.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_c4.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_c8.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_i16.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_i4.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_i8.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_r10.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_r16.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_r4.c | 47 | ||||
-rw-r--r-- | libgfortran/generated/matmul_r8.c | 47 |
11 files changed, 484 insertions, 33 deletions
diff --git a/libgfortran/generated/matmul_c10.c b/libgfortran/generated/matmul_c10.c index df2cd93c15f..5e3b281245c 100644 --- a/libgfortran/generated/matmul_c10.c +++ b/libgfortran/generated/matmul_c10.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_COMPLEX_10) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_COMPLEX_10 *, const GFC_COMPLEX_10 *, + const int *, const GFC_COMPLEX_10 *, const int *, + const GFC_COMPLEX_10 *, GFC_COMPLEX_10 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_c10 (gfc_array_c10 * const restrict retarray, - gfc_array_c10 * const restrict a, gfc_array_c10 * const restrict b); + gfc_array_c10 * const restrict a, gfc_array_c10 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_c10); void matmul_c10 (gfc_array_c10 * const restrict retarray, - gfc_array_c10 * const restrict a, gfc_array_c10 * const restrict b) + gfc_array_c10 * const restrict a, gfc_array_c10 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_COMPLEX_10 * restrict abase; const GFC_COMPLEX_10 * restrict bbase; @@ -177,6 +193,31 @@ matmul_c10 (gfc_array_c10 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_COMPLEX_10 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_COMPLEX_10 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_c16.c b/libgfortran/generated/matmul_c16.c index 6425eb8d49d..f7301114b37 100644 --- a/libgfortran/generated/matmul_c16.c +++ b/libgfortran/generated/matmul_c16.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_COMPLEX_16) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_COMPLEX_16 *, const GFC_COMPLEX_16 *, + const int *, const GFC_COMPLEX_16 *, const int *, + const GFC_COMPLEX_16 *, GFC_COMPLEX_16 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_c16 (gfc_array_c16 * const restrict retarray, - gfc_array_c16 * const restrict a, gfc_array_c16 * const restrict b); + gfc_array_c16 * const restrict a, gfc_array_c16 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_c16); void matmul_c16 (gfc_array_c16 * const restrict retarray, - gfc_array_c16 * const restrict a, gfc_array_c16 * const restrict b) + gfc_array_c16 * const restrict a, gfc_array_c16 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_COMPLEX_16 * restrict abase; const GFC_COMPLEX_16 * restrict bbase; @@ -177,6 +193,31 @@ matmul_c16 (gfc_array_c16 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_COMPLEX_16 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_COMPLEX_16 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_c4.c b/libgfortran/generated/matmul_c4.c index 2d47a134972..f2984ab48ab 100644 --- a/libgfortran/generated/matmul_c4.c +++ b/libgfortran/generated/matmul_c4.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_COMPLEX_4) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_COMPLEX_4 *, const GFC_COMPLEX_4 *, + const int *, const GFC_COMPLEX_4 *, const int *, + const GFC_COMPLEX_4 *, GFC_COMPLEX_4 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_c4 (gfc_array_c4 * const restrict retarray, - gfc_array_c4 * const restrict a, gfc_array_c4 * const restrict b); + gfc_array_c4 * const restrict a, gfc_array_c4 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_c4); void matmul_c4 (gfc_array_c4 * const restrict retarray, - gfc_array_c4 * const restrict a, gfc_array_c4 * const restrict b) + gfc_array_c4 * const restrict a, gfc_array_c4 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_COMPLEX_4 * restrict abase; const GFC_COMPLEX_4 * restrict bbase; @@ -177,6 +193,31 @@ matmul_c4 (gfc_array_c4 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_COMPLEX_4 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_COMPLEX_4 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_c8.c b/libgfortran/generated/matmul_c8.c index f22719df505..65cc0a52c4b 100644 --- a/libgfortran/generated/matmul_c8.c +++ b/libgfortran/generated/matmul_c8.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_COMPLEX_8) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_COMPLEX_8 *, const GFC_COMPLEX_8 *, + const int *, const GFC_COMPLEX_8 *, const int *, + const GFC_COMPLEX_8 *, GFC_COMPLEX_8 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_c8 (gfc_array_c8 * const restrict retarray, - gfc_array_c8 * const restrict a, gfc_array_c8 * const restrict b); + gfc_array_c8 * const restrict a, gfc_array_c8 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_c8); void matmul_c8 (gfc_array_c8 * const restrict retarray, - gfc_array_c8 * const restrict a, gfc_array_c8 * const restrict b) + gfc_array_c8 * const restrict a, gfc_array_c8 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_COMPLEX_8 * restrict abase; const GFC_COMPLEX_8 * restrict bbase; @@ -177,6 +193,31 @@ matmul_c8 (gfc_array_c8 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_COMPLEX_8 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_COMPLEX_8 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_i16.c b/libgfortran/generated/matmul_i16.c index 73c3fbc108d..a193669d108 100644 --- a/libgfortran/generated/matmul_i16.c +++ b/libgfortran/generated/matmul_i16.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_INTEGER_16) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_INTEGER_16 *, const GFC_INTEGER_16 *, + const int *, const GFC_INTEGER_16 *, const int *, + const GFC_INTEGER_16 *, GFC_INTEGER_16 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_i16 (gfc_array_i16 * const restrict retarray, - gfc_array_i16 * const restrict a, gfc_array_i16 * const restrict b); + gfc_array_i16 * const restrict a, gfc_array_i16 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_i16); void matmul_i16 (gfc_array_i16 * const restrict retarray, - gfc_array_i16 * const restrict a, gfc_array_i16 * const restrict b) + gfc_array_i16 * const restrict a, gfc_array_i16 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_INTEGER_16 * restrict abase; const GFC_INTEGER_16 * restrict bbase; @@ -177,6 +193,31 @@ matmul_i16 (gfc_array_i16 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_INTEGER_16 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_INTEGER_16 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_i4.c b/libgfortran/generated/matmul_i4.c index 63bca0152cd..69b9b487a81 100644 --- a/libgfortran/generated/matmul_i4.c +++ b/libgfortran/generated/matmul_i4.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_INTEGER_4) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_INTEGER_4 *, const GFC_INTEGER_4 *, + const int *, const GFC_INTEGER_4 *, const int *, + const GFC_INTEGER_4 *, GFC_INTEGER_4 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_i4 (gfc_array_i4 * const restrict retarray, - gfc_array_i4 * const restrict a, gfc_array_i4 * const restrict b); + gfc_array_i4 * const restrict a, gfc_array_i4 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_i4); void matmul_i4 (gfc_array_i4 * const restrict retarray, - gfc_array_i4 * const restrict a, gfc_array_i4 * const restrict b) + gfc_array_i4 * const restrict a, gfc_array_i4 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_INTEGER_4 * restrict abase; const GFC_INTEGER_4 * restrict bbase; @@ -177,6 +193,31 @@ matmul_i4 (gfc_array_i4 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_INTEGER_4 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_INTEGER_4 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_i8.c b/libgfortran/generated/matmul_i8.c index caaf9e8f976..23a87a904f7 100644 --- a/libgfortran/generated/matmul_i8.c +++ b/libgfortran/generated/matmul_i8.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_INTEGER_8) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_INTEGER_8 *, const GFC_INTEGER_8 *, + const int *, const GFC_INTEGER_8 *, const int *, + const GFC_INTEGER_8 *, GFC_INTEGER_8 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_i8 (gfc_array_i8 * const restrict retarray, - gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b); + gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_i8); void matmul_i8 (gfc_array_i8 * const restrict retarray, - gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b) + gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_INTEGER_8 * restrict abase; const GFC_INTEGER_8 * restrict bbase; @@ -177,6 +193,31 @@ matmul_i8 (gfc_array_i8 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_INTEGER_8 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_INTEGER_8 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_r10.c b/libgfortran/generated/matmul_r10.c index 8fa1d6d9e49..e4dfd74ef03 100644 --- a/libgfortran/generated/matmul_r10.c +++ b/libgfortran/generated/matmul_r10.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_REAL_10) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_REAL_10 *, const GFC_REAL_10 *, + const int *, const GFC_REAL_10 *, const int *, + const GFC_REAL_10 *, GFC_REAL_10 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_r10 (gfc_array_r10 * const restrict retarray, - gfc_array_r10 * const restrict a, gfc_array_r10 * const restrict b); + gfc_array_r10 * const restrict a, gfc_array_r10 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_r10); void matmul_r10 (gfc_array_r10 * const restrict retarray, - gfc_array_r10 * const restrict a, gfc_array_r10 * const restrict b) + gfc_array_r10 * const restrict a, gfc_array_r10 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_REAL_10 * restrict abase; const GFC_REAL_10 * restrict bbase; @@ -177,6 +193,31 @@ matmul_r10 (gfc_array_r10 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_REAL_10 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_REAL_10 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_r16.c b/libgfortran/generated/matmul_r16.c index 0f61b038168..ec760f2d3d8 100644 --- a/libgfortran/generated/matmul_r16.c +++ b/libgfortran/generated/matmul_r16.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_REAL_16) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_REAL_16 *, const GFC_REAL_16 *, + const int *, const GFC_REAL_16 *, const int *, + const GFC_REAL_16 *, GFC_REAL_16 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_r16 (gfc_array_r16 * const restrict retarray, - gfc_array_r16 * const restrict a, gfc_array_r16 * const restrict b); + gfc_array_r16 * const restrict a, gfc_array_r16 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_r16); void matmul_r16 (gfc_array_r16 * const restrict retarray, - gfc_array_r16 * const restrict a, gfc_array_r16 * const restrict b) + gfc_array_r16 * const restrict a, gfc_array_r16 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_REAL_16 * restrict abase; const GFC_REAL_16 * restrict bbase; @@ -177,6 +193,31 @@ matmul_r16 (gfc_array_r16 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_REAL_16 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_REAL_16 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_r4.c b/libgfortran/generated/matmul_r4.c index d684dd2905c..cf2f45fb125 100644 --- a/libgfortran/generated/matmul_r4.c +++ b/libgfortran/generated/matmul_r4.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_REAL_4) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_REAL_4 *, const GFC_REAL_4 *, + const int *, const GFC_REAL_4 *, const int *, + const GFC_REAL_4 *, GFC_REAL_4 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_r4 (gfc_array_r4 * const restrict retarray, - gfc_array_r4 * const restrict a, gfc_array_r4 * const restrict b); + gfc_array_r4 * const restrict a, gfc_array_r4 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_r4); void matmul_r4 (gfc_array_r4 * const restrict retarray, - gfc_array_r4 * const restrict a, gfc_array_r4 * const restrict b) + gfc_array_r4 * const restrict a, gfc_array_r4 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_REAL_4 * restrict abase; const GFC_REAL_4 * restrict bbase; @@ -177,6 +193,31 @@ matmul_r4 (gfc_array_r4 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_REAL_4 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_REAL_4 * restrict bbase_y; diff --git a/libgfortran/generated/matmul_r8.c b/libgfortran/generated/matmul_r8.c index 41726bce2a5..c746f6c3519 100644 --- a/libgfortran/generated/matmul_r8.c +++ b/libgfortran/generated/matmul_r8.c @@ -36,6 +36,16 @@ Boston, MA 02110-1301, USA. */ #if defined (HAVE_GFC_REAL_8) +/* 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 + matrices. */ + +typedef void (*blas_call)(const char *, const char *, const int *, const int *, + const int *, const GFC_REAL_8 *, const GFC_REAL_8 *, + const int *, const GFC_REAL_8 *, const int *, + const GFC_REAL_8 *, GFC_REAL_8 *, const int *, + int, int); + /* 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: @@ -56,18 +66,24 @@ Boston, MA 02110-1301, USA. */ DO I=1,M S = 0 DO K=1,COUNT - S = S+A(I,K)+B(K,J) + S = S+A(I,K)*B(K,J) C(I,J) = S ENDIF */ +/* 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_r8 (gfc_array_r8 * const restrict retarray, - gfc_array_r8 * const restrict a, gfc_array_r8 * const restrict b); + gfc_array_r8 * const restrict a, gfc_array_r8 * const restrict b, int try_blas, + int blas_limit, blas_call gemm); export_proto(matmul_r8); void matmul_r8 (gfc_array_r8 * const restrict retarray, - gfc_array_r8 * const restrict a, gfc_array_r8 * const restrict b) + gfc_array_r8 * const restrict a, gfc_array_r8 * const restrict b, int try_blas, + int blas_limit, blas_call gemm) { const GFC_REAL_8 * restrict abase; const GFC_REAL_8 * restrict bbase; @@ -177,6 +193,31 @@ matmul_r8 (gfc_array_r8 * const restrict retarray, bbase = b->data; dest = retarray->data; + + /* Now that everything is set up, we're performing the multiplication + itself. */ + +#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x))) + + 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 GFC_REAL_8 one = 1, zero = 0; + const int lda = (axstride == 1) ? aystride : axstride, + ldb = (bxstride == 1) ? bystride : bxstride; + + 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; + } + } + if (rxstride == 1 && axstride == 1 && bxstride == 1) { const GFC_REAL_8 * restrict bbase_y; |