diff options
Diffstat (limited to 'libgfortran')
| -rw-r--r-- | libgfortran/ChangeLog | 16 | ||||
| -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 | ||||
| -rw-r--r-- | libgfortran/m4/matmul.m4 | 47 |
13 files changed, 544 insertions, 36 deletions
diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog index a9e70825f80..063c62519d3 100644 --- a/libgfortran/ChangeLog +++ b/libgfortran/ChangeLog @@ -1,3 +1,19 @@ +2006-10-22 Francois-Xavier Coudert <coudert@clipper.ens.fr> + + PR fortran/26025 + * m4/matmul.m4: Add possible call to gemm routine. + * generated/matmul_r8.c: Regenerate. + * generated/matmul_r16.c: Regenerate. + * generated/matmul_c8.c: Regenerate. + * generated/matmul_i8.c: Regenerate. + * generated/matmul_c16.c: Regenerate. + * generated/matmul_r10.c: Regenerate. + * generated/matmul_r4.c: Regenerate. + * generated/matmul_c10.c: Regenerate. + * generated/matmul_c4.c: Regenerate. + * generated/matmul_i4.c: Regenerate. + * generated/matmul_i16.c: Regenerate. + 2006-10-21 Steven G. Kargl <kargl@gcc.gnu.org> * runtime/error.c: Add errno.h 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; diff --git a/libgfortran/m4/matmul.m4 b/libgfortran/m4/matmul.m4 index 3678c639f2a..ef2f0fb88dc 100644 --- a/libgfortran/m4/matmul.m4 +++ b/libgfortran/m4/matmul.m4 @@ -37,6 +37,16 @@ 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 + matrices. */ + +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); + /* 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: @@ -57,18 +67,24 @@ include(iparm.m4)dnl 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_`'rtype_code (rtype * const restrict retarray, - rtype * const restrict a, rtype * const restrict b); + rtype * const restrict a, rtype * const restrict b, int try_blas, + 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) + 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; @@ -179,6 +195,31 @@ sinclude(`matmul_asm_'rtype_code`.m4')dnl 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 rtype_name 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 rtype_name * restrict bbase_y; |