diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog
index ec72c6daa84bac6c63591b4ec432dde2473d2acd..ec689fd7285f385f8d36565013031050e5a8624a 100644
--- a/libgfortran/ChangeLog
+++ b/libgfortran/ChangeLog
@@ -1,3 +1,22 @@
+2017-03-02  Thomas Koenig  <tkoenig@gcc.gnu.org>
+	    Jakub Jelinek  <jakub@redhat.com>
+
+	* m4/matmul.m4 (matmul_'rtype_code`):  Avoid
+	race condition on storing function pointer.
+        * generated/matmul_c10.c: Regenerated.
+        * generated/matmul_c16.c: Regenerated.
+        * generated/matmul_c4.c: Regenerated.
+        * generated/matmul_c8.c: Regenerated.
+        * generated/matmul_i1.c: Regenerated.
+        * generated/matmul_i16.c: Regenerated.
+        * generated/matmul_i2.c: Regenerated.
+        * generated/matmul_i4.c: Regenerated.
+        * generated/matmul_i8.c: Regenerated.
+        * generated/matmul_r10.c: Regenerated.
+        * generated/matmul_r16.c: Regenerated.
+        * generated/matmul_r4.c: Regenerated.
+        * generated/matmul_r8.c: Regenerated.
+
 2017-03-02  Thomas Koenig  <tkoenig@gcc.gnu.org>
 
 	PR fortran/78379
diff --git a/libgfortran/generated/matmul_c10.c b/libgfortran/generated/matmul_c10.c
index b333a844ea555c641692f19bd8f6fcf7ba7a3391..c8e7a6c37fbddd520c8fc635e213984e1bf289dd 100644
--- a/libgfortran/generated/matmul_c10.c
+++ b/libgfortran/generated/matmul_c10.c
@@ -2256,19 +2256,24 @@ void matmul_c10 (gfc_array_c10 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_c10 * const restrict retarray, 
 	gfc_array_c10 * const restrict a, gfc_array_c10 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_c10 * const restrict retarray, 
+	gfc_array_c10 * const restrict a, gfc_array_c10 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_c10_vanilla;
+      matmul_fn = matmul_c10_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_c10_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_c10_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_c10 (gfc_array_c10 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_c10_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_c10_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_c10 (gfc_array_c10 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_c10_avx;
-	      goto tailcall;
+              matmul_fn = matmul_c10_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_c16.c b/libgfortran/generated/matmul_c16.c
index 0ef66c0f4a2c167b45ad9d4e4fe01a7841d9afdc..adbf1c81a7ae4b984efa3ee65b9c2e660acc1395 100644
--- a/libgfortran/generated/matmul_c16.c
+++ b/libgfortran/generated/matmul_c16.c
@@ -2256,19 +2256,24 @@ void matmul_c16 (gfc_array_c16 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_c16 * const restrict retarray, 
 	gfc_array_c16 * const restrict a, gfc_array_c16 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_c16 * const restrict retarray, 
+	gfc_array_c16 * const restrict a, gfc_array_c16 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_c16_vanilla;
+      matmul_fn = matmul_c16_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_c16_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_c16_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_c16 (gfc_array_c16 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_c16_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_c16_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_c16 (gfc_array_c16 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_c16_avx;
-	      goto tailcall;
+              matmul_fn = matmul_c16_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_c4.c b/libgfortran/generated/matmul_c4.c
index b30320b37ae15a684c58aa11e29c9a84ba9c57ca..ff0e76a9275ed4f71848a76dcac35aa8a45cd4cd 100644
--- a/libgfortran/generated/matmul_c4.c
+++ b/libgfortran/generated/matmul_c4.c
@@ -2256,19 +2256,24 @@ void matmul_c4 (gfc_array_c4 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_c4 * const restrict retarray, 
 	gfc_array_c4 * const restrict a, gfc_array_c4 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_c4 * const restrict retarray, 
+	gfc_array_c4 * const restrict a, gfc_array_c4 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_c4_vanilla;
+      matmul_fn = matmul_c4_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_c4_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_c4_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_c4 (gfc_array_c4 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_c4_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_c4_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_c4 (gfc_array_c4 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_c4_avx;
-	      goto tailcall;
+              matmul_fn = matmul_c4_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_c8.c b/libgfortran/generated/matmul_c8.c
index 75b4680c3a0a77a07e54769545f1e59b12f43304..483ebf114d718733534dfd598ee51e2980d96f2e 100644
--- a/libgfortran/generated/matmul_c8.c
+++ b/libgfortran/generated/matmul_c8.c
@@ -2256,19 +2256,24 @@ void matmul_c8 (gfc_array_c8 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_c8 * const restrict retarray, 
 	gfc_array_c8 * const restrict a, gfc_array_c8 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_c8 * const restrict retarray, 
+	gfc_array_c8 * const restrict a, gfc_array_c8 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_c8_vanilla;
+      matmul_fn = matmul_c8_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_c8_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_c8_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_c8 (gfc_array_c8 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_c8_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_c8_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_c8 (gfc_array_c8 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_c8_avx;
-	      goto tailcall;
+              matmul_fn = matmul_c8_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_i1.c b/libgfortran/generated/matmul_i1.c
index 924826338d8d909b125f4e0f4c0da1a6ce95f3e0..51cfcc05360728ebf1bfc9f95edfe607158ff223 100644
--- a/libgfortran/generated/matmul_i1.c
+++ b/libgfortran/generated/matmul_i1.c
@@ -2256,19 +2256,24 @@ void matmul_i1 (gfc_array_i1 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_i1 * const restrict retarray, 
 	gfc_array_i1 * const restrict a, gfc_array_i1 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_i1 * const restrict retarray, 
+	gfc_array_i1 * const restrict a, gfc_array_i1 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_i1_vanilla;
+      matmul_fn = matmul_i1_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_i1_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_i1_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_i1 (gfc_array_i1 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_i1_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_i1_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_i1 (gfc_array_i1 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_i1_avx;
-	      goto tailcall;
+              matmul_fn = matmul_i1_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_i16.c b/libgfortran/generated/matmul_i16.c
index f10540ed48a1daca20d61bf8869bb464bf9527a4..90a422f42900a914cb52a55c72c1eb010c75abb9 100644
--- a/libgfortran/generated/matmul_i16.c
+++ b/libgfortran/generated/matmul_i16.c
@@ -2256,19 +2256,24 @@ void matmul_i16 (gfc_array_i16 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_i16 * const restrict retarray, 
 	gfc_array_i16 * const restrict a, gfc_array_i16 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_i16 * const restrict retarray, 
+	gfc_array_i16 * const restrict a, gfc_array_i16 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_i16_vanilla;
+      matmul_fn = matmul_i16_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_i16_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_i16_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_i16 (gfc_array_i16 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_i16_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_i16_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_i16 (gfc_array_i16 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_i16_avx;
-	      goto tailcall;
+              matmul_fn = matmul_i16_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_i2.c b/libgfortran/generated/matmul_i2.c
index 55ad5c614e66eee2825f1be0212a9b1b30f6cbe9..58be84fc9c82b43d46d798e7c8b06a18aa193dab 100644
--- a/libgfortran/generated/matmul_i2.c
+++ b/libgfortran/generated/matmul_i2.c
@@ -2256,19 +2256,24 @@ void matmul_i2 (gfc_array_i2 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_i2 * const restrict retarray, 
 	gfc_array_i2 * const restrict a, gfc_array_i2 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_i2 * const restrict retarray, 
+	gfc_array_i2 * const restrict a, gfc_array_i2 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_i2_vanilla;
+      matmul_fn = matmul_i2_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_i2_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_i2_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_i2 (gfc_array_i2 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_i2_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_i2_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_i2 (gfc_array_i2 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_i2_avx;
-	      goto tailcall;
+              matmul_fn = matmul_i2_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_i4.c b/libgfortran/generated/matmul_i4.c
index 97b4a5b6aa0889383561aa0547d5aa8cb2432538..1e724d54ffafa1e8c78cf18502a67bd275668110 100644
--- a/libgfortran/generated/matmul_i4.c
+++ b/libgfortran/generated/matmul_i4.c
@@ -2256,19 +2256,24 @@ void matmul_i4 (gfc_array_i4 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_i4 * const restrict retarray, 
 	gfc_array_i4 * const restrict a, gfc_array_i4 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_i4 * const restrict retarray, 
+	gfc_array_i4 * const restrict a, gfc_array_i4 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_i4_vanilla;
+      matmul_fn = matmul_i4_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_i4_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_i4_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_i4 (gfc_array_i4 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_i4_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_i4_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_i4 (gfc_array_i4 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_i4_avx;
-	      goto tailcall;
+              matmul_fn = matmul_i4_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_i8.c b/libgfortran/generated/matmul_i8.c
index ae78ecfccb68af32deaab2aa38915d3b33d10c49..f5ed84fa2885fa98d90035c4dbe934ba9ed9e346 100644
--- a/libgfortran/generated/matmul_i8.c
+++ b/libgfortran/generated/matmul_i8.c
@@ -2256,19 +2256,24 @@ void matmul_i8 (gfc_array_i8 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_i8 * const restrict retarray, 
 	gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_i8 * const restrict retarray, 
+	gfc_array_i8 * const restrict a, gfc_array_i8 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_i8_vanilla;
+      matmul_fn = matmul_i8_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_i8_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_i8_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_i8 (gfc_array_i8 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_i8_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_i8_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_i8 (gfc_array_i8 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_i8_avx;
-	      goto tailcall;
+              matmul_fn = matmul_i8_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_r10.c b/libgfortran/generated/matmul_r10.c
index 11d059198c92b0f1dfb0483145d040f5567bf2b9..8f76faceabff6aa988b6f343dc3b7430e015989b 100644
--- a/libgfortran/generated/matmul_r10.c
+++ b/libgfortran/generated/matmul_r10.c
@@ -2256,19 +2256,24 @@ void matmul_r10 (gfc_array_r10 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_r10 * const restrict retarray, 
 	gfc_array_r10 * const restrict a, gfc_array_r10 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_r10 * const restrict retarray, 
+	gfc_array_r10 * const restrict a, gfc_array_r10 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_r10_vanilla;
+      matmul_fn = matmul_r10_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_r10_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_r10_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_r10 (gfc_array_r10 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_r10_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_r10_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_r10 (gfc_array_r10 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_r10_avx;
-	      goto tailcall;
+              matmul_fn = matmul_r10_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_r16.c b/libgfortran/generated/matmul_r16.c
index 73e7c9877adf09c4ff13f593aa44e0c478266f49..6a7687fdc33731a304dcdd9926c8264783338e3f 100644
--- a/libgfortran/generated/matmul_r16.c
+++ b/libgfortran/generated/matmul_r16.c
@@ -2256,19 +2256,24 @@ void matmul_r16 (gfc_array_r16 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_r16 * const restrict retarray, 
 	gfc_array_r16 * const restrict a, gfc_array_r16 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_r16 * const restrict retarray, 
+	gfc_array_r16 * const restrict a, gfc_array_r16 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_r16_vanilla;
+      matmul_fn = matmul_r16_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_r16_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_r16_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_r16 (gfc_array_r16 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_r16_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_r16_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_r16 (gfc_array_r16 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_r16_avx;
-	      goto tailcall;
+              matmul_fn = matmul_r16_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_r4.c b/libgfortran/generated/matmul_r4.c
index ac7306fe035f8830c4475b85aba2e00355c8c443..dbb31b05c3bc9d1139cbfc6f4e93b641e376878f 100644
--- a/libgfortran/generated/matmul_r4.c
+++ b/libgfortran/generated/matmul_r4.c
@@ -2256,19 +2256,24 @@ void matmul_r4 (gfc_array_r4 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_r4 * const restrict retarray, 
 	gfc_array_r4 * const restrict a, gfc_array_r4 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_r4 * const restrict retarray, 
+	gfc_array_r4 * const restrict a, gfc_array_r4 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_r4_vanilla;
+      matmul_fn = matmul_r4_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_r4_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_r4_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_r4 (gfc_array_r4 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_r4_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_r4_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_r4 (gfc_array_r4 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_r4_avx;
-	      goto tailcall;
+              matmul_fn = matmul_r4_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/generated/matmul_r8.c b/libgfortran/generated/matmul_r8.c
index 8d2e784de60753f209e94c4c7c34211c5e56d1d4..23efed22e2ab91d9b164a3d1fba8a8eeb62988ec 100644
--- a/libgfortran/generated/matmul_r8.c
+++ b/libgfortran/generated/matmul_r8.c
@@ -2256,19 +2256,24 @@ void matmul_r8 (gfc_array_r8 * const restrict retarray,
 {
   static void (*matmul_p) (gfc_array_r8 * const restrict retarray, 
 	gfc_array_r8 * const restrict a, gfc_array_r8 * const restrict b, int try_blas,
-	int blas_limit, blas_call gemm) = NULL;
+	int blas_limit, blas_call gemm);
+
+  void (*matmul_fn) (gfc_array_r8 * const restrict retarray, 
+	gfc_array_r8 * const restrict a, gfc_array_r8 * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
 
-  if (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_r8_vanilla;
+      matmul_fn = matmul_r8_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_r8_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_r8_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -2277,8 +2282,8 @@ void matmul_r8 (gfc_array_r8 * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_r8_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_r8_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -2286,15 +2291,16 @@ void matmul_r8 (gfc_array_r8 * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_r8_avx;
-	      goto tailcall;
+              matmul_fn = matmul_r8_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */
diff --git a/libgfortran/m4/matmul.m4 b/libgfortran/m4/matmul.m4
index 812a7e7e5714f2c741d051d1dfd486f61baa5327..7976fda8bb4610c8169a142fcb5c444a361e7bd5 100644
--- a/libgfortran/m4/matmul.m4
+++ b/libgfortran/m4/matmul.m4
@@ -121,19 +121,24 @@ void matmul_'rtype_code` ('rtype` * const restrict retarray,
 {
   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) = NULL;
+	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 (matmul_p == NULL)
+  matmul_fn = __atomic_load_n (&matmul_p, __ATOMIC_RELAXED);
+  if (matmul_fn == NULL)
     {
-      matmul_p = matmul_'rtype_code`_vanilla;
+      matmul_fn = matmul_'rtype_code`_vanilla;
       if (__cpu_model.__cpu_vendor == VENDOR_INTEL)
 	{
           /* Run down the available processors in order of preference.  */
 #ifdef HAVE_AVX512F
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX512F))
 	    {
-	      matmul_p = matmul_'rtype_code`_avx512f;
-	      goto tailcall;
+	      matmul_fn = matmul_'rtype_code`_avx512f;
+	      goto store;
 	    }
 
 #endif  /* HAVE_AVX512F */
@@ -142,8 +147,8 @@ void matmul_'rtype_code` ('rtype` * const restrict retarray,
       	  if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX2))
 	     && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
 	    {
-	      matmul_p = matmul_'rtype_code`_avx2;
-	      goto tailcall;
+	      matmul_fn = matmul_'rtype_code`_avx2;
+	      goto store;
 	    }
 
 #endif
@@ -151,15 +156,16 @@ void matmul_'rtype_code` ('rtype` * const restrict retarray,
 #ifdef HAVE_AVX
       	  if (__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
  	    {
-              matmul_p = matmul_'rtype_code`_avx;
-	      goto tailcall;
+              matmul_fn = matmul_'rtype_code`_avx;
+	      goto store;
 	    }
 #endif  /* HAVE_AVX */
         }
+   store:
+      __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
 
-tailcall:
-   (*matmul_p) (retarray, a, b, try_blas, blas_limit, gemm);
+   (*matmul_fn) (retarray, a, b, try_blas, blas_limit, gemm);
 }
 
 #else  /* Just the vanilla function.  */