diff options
Diffstat (limited to 'libgfortran/generated/matmul_r4.c')
| -rw-r--r-- | libgfortran/generated/matmul_r4.c | 47 |
1 files changed, 44 insertions, 3 deletions
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; |