]> git.ipfire.org Git - thirdparty/gcc.git/blobdiff - libgfortran/generated/matmul_c8.c
Update copyright years.
[thirdparty/gcc.git] / libgfortran / generated / matmul_c8.c
index aafa3be182124ce444e6195367de57b880e4ca09..8c19b49deaecffbc29c92b193778d00b399bb074 100644 (file)
@@ -1,5 +1,5 @@
 /* Implementation of the MATMUL intrinsic
-   Copyright (C) 2002-2017 Free Software Foundation, Inc.
+   Copyright (C) 2002-2020 Free Software Foundation, Inc.
    Contributed by Paul Brook <paul@nowt.org>
 
 This file is part of the GNU Fortran runtime library (libgfortran).
@@ -144,8 +144,8 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -153,8 +153,8 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else
@@ -162,17 +162,15 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 1:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
 
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 2:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 2 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
     }
@@ -213,7 +211,9 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
   if (count != GFC_DESCRIPTOR_EXTENT(b,0))
     {
       if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
-       runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
+       runtime_error ("Incorrect extent in argument B in MATMUL intrinsic "
+                      "in dimension 1: is %ld, should be %ld",
+                      (long int) GFC_DESCRIPTOR_EXTENT(b,0), (long int) count);
     }
 
   if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -222,9 +222,9 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
       bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
 
       /* bystride should never be used for 1-dimensional b.
-        in case it is we want it to cause a segfault, rather than
-        an incorrect result. */
-      bystride = 0xDEADBEEF;
+         The value is only used for calculation of the
+         memory by the buffer.  */
+      bystride = 256;
       ycount = 1;
     }
   else
@@ -258,7 +258,18 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
       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,
+         const char *transa, *transb;
+         if (try_blas & 2)
+           transa = "C";
+         else
+           transa = axstride == 1 ? "N" : "T";
+
+         if (try_blas & 4)
+           transb = "C";
+         else
+           transb = bxstride == 1 ? "N" : "T";
+
+         gemm (transa, transb , &m,
                &n, &k, &one, abase, &lda, bbase, &ldb, &zero, dest,
                &ldc, 1, 1);
          return;
@@ -290,6 +301,7 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
                 f13, f14, f23, f24, f33, f34, f43, f44;
       index_type i, j, l, ii, jj, ll;
       index_type isec, jsec, lsec, uisec, ujsec, ulsec;
+      GFC_COMPLEX_8 *t1;
 
       a = abase;
       b = bbase;
@@ -306,25 +318,27 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
       b_offset = 1 + b_dim1;
       b -= b_offset;
 
+      /* Empty c first.  */
+      for (j=1; j<=n; j++)
+       for (i=1; i<=m; i++)
+         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+
       /* Early exit if possible */
       if (m == 0 || n == 0 || k == 0)
        return;
 
       /* Adjust size of t1 to what is needed.  */
-      index_type t1_dim;
-      t1_dim = (a_dim1-1) * 256 + b_dim1;
+      index_type t1_dim, a_sz;
+      if (aystride == 1)
+        a_sz = rystride;
+      else
+        a_sz = a_dim1;
+
+      t1_dim = a_sz * 256 + b_dim1;
       if (t1_dim > 65536)
        t1_dim = 65536;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wvla"
-      GFC_COMPLEX_8 t1[t1_dim]; /* was [256][256] */
-#pragma GCC diagnostic pop
-
-      /* Empty c first.  */
-      for (j=1; j<=n; j++)
-       for (i=1; i<=m; i++)
-         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+      t1 = malloc (t1_dim * sizeof(GFC_COMPLEX_8));
 
       /* Start turning the crank. */
       i1 = n;
@@ -535,6 +549,7 @@ matmul_c8_avx (gfc_array_c8 * const restrict retarray,
                }
            }
        }
+      free(t1);
       return;
     }
   else if (rxstride == 1 && aystride == 1 && bxstride == 1)
@@ -697,8 +712,8 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -706,8 +721,8 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else
@@ -715,17 +730,15 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 1:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
 
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 2:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 2 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
     }
@@ -766,7 +779,9 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
   if (count != GFC_DESCRIPTOR_EXTENT(b,0))
     {
       if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
-       runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
+       runtime_error ("Incorrect extent in argument B in MATMUL intrinsic "
+                      "in dimension 1: is %ld, should be %ld",
+                      (long int) GFC_DESCRIPTOR_EXTENT(b,0), (long int) count);
     }
 
   if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -775,9 +790,9 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
       bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
 
       /* bystride should never be used for 1-dimensional b.
-        in case it is we want it to cause a segfault, rather than
-        an incorrect result. */
-      bystride = 0xDEADBEEF;
+         The value is only used for calculation of the
+         memory by the buffer.  */
+      bystride = 256;
       ycount = 1;
     }
   else
@@ -811,7 +826,18 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
       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,
+         const char *transa, *transb;
+         if (try_blas & 2)
+           transa = "C";
+         else
+           transa = axstride == 1 ? "N" : "T";
+
+         if (try_blas & 4)
+           transb = "C";
+         else
+           transb = bxstride == 1 ? "N" : "T";
+
+         gemm (transa, transb , &m,
                &n, &k, &one, abase, &lda, bbase, &ldb, &zero, dest,
                &ldc, 1, 1);
          return;
@@ -843,6 +869,7 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
                 f13, f14, f23, f24, f33, f34, f43, f44;
       index_type i, j, l, ii, jj, ll;
       index_type isec, jsec, lsec, uisec, ujsec, ulsec;
+      GFC_COMPLEX_8 *t1;
 
       a = abase;
       b = bbase;
@@ -859,25 +886,27 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
       b_offset = 1 + b_dim1;
       b -= b_offset;
 
+      /* Empty c first.  */
+      for (j=1; j<=n; j++)
+       for (i=1; i<=m; i++)
+         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+
       /* Early exit if possible */
       if (m == 0 || n == 0 || k == 0)
        return;
 
       /* Adjust size of t1 to what is needed.  */
-      index_type t1_dim;
-      t1_dim = (a_dim1-1) * 256 + b_dim1;
+      index_type t1_dim, a_sz;
+      if (aystride == 1)
+        a_sz = rystride;
+      else
+        a_sz = a_dim1;
+
+      t1_dim = a_sz * 256 + b_dim1;
       if (t1_dim > 65536)
        t1_dim = 65536;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wvla"
-      GFC_COMPLEX_8 t1[t1_dim]; /* was [256][256] */
-#pragma GCC diagnostic pop
-
-      /* Empty c first.  */
-      for (j=1; j<=n; j++)
-       for (i=1; i<=m; i++)
-         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+      t1 = malloc (t1_dim * sizeof(GFC_COMPLEX_8));
 
       /* Start turning the crank. */
       i1 = n;
@@ -1088,6 +1117,7 @@ matmul_c8_avx2 (gfc_array_c8 * const restrict retarray,
                }
            }
        }
+      free(t1);
       return;
     }
   else if (rxstride == 1 && aystride == 1 && bxstride == 1)
@@ -1250,8 +1280,8 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -1259,8 +1289,8 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else
@@ -1268,17 +1298,15 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 1:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
 
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 2:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 2 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
     }
@@ -1319,7 +1347,9 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
   if (count != GFC_DESCRIPTOR_EXTENT(b,0))
     {
       if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
-       runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
+       runtime_error ("Incorrect extent in argument B in MATMUL intrinsic "
+                      "in dimension 1: is %ld, should be %ld",
+                      (long int) GFC_DESCRIPTOR_EXTENT(b,0), (long int) count);
     }
 
   if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -1328,9 +1358,9 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
       bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
 
       /* bystride should never be used for 1-dimensional b.
-        in case it is we want it to cause a segfault, rather than
-        an incorrect result. */
-      bystride = 0xDEADBEEF;
+         The value is only used for calculation of the
+         memory by the buffer.  */
+      bystride = 256;
       ycount = 1;
     }
   else
@@ -1364,7 +1394,18 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
       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,
+         const char *transa, *transb;
+         if (try_blas & 2)
+           transa = "C";
+         else
+           transa = axstride == 1 ? "N" : "T";
+
+         if (try_blas & 4)
+           transb = "C";
+         else
+           transb = bxstride == 1 ? "N" : "T";
+
+         gemm (transa, transb , &m,
                &n, &k, &one, abase, &lda, bbase, &ldb, &zero, dest,
                &ldc, 1, 1);
          return;
@@ -1396,6 +1437,7 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
                 f13, f14, f23, f24, f33, f34, f43, f44;
       index_type i, j, l, ii, jj, ll;
       index_type isec, jsec, lsec, uisec, ujsec, ulsec;
+      GFC_COMPLEX_8 *t1;
 
       a = abase;
       b = bbase;
@@ -1412,25 +1454,27 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
       b_offset = 1 + b_dim1;
       b -= b_offset;
 
+      /* Empty c first.  */
+      for (j=1; j<=n; j++)
+       for (i=1; i<=m; i++)
+         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+
       /* Early exit if possible */
       if (m == 0 || n == 0 || k == 0)
        return;
 
       /* Adjust size of t1 to what is needed.  */
-      index_type t1_dim;
-      t1_dim = (a_dim1-1) * 256 + b_dim1;
+      index_type t1_dim, a_sz;
+      if (aystride == 1)
+        a_sz = rystride;
+      else
+        a_sz = a_dim1;
+
+      t1_dim = a_sz * 256 + b_dim1;
       if (t1_dim > 65536)
        t1_dim = 65536;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wvla"
-      GFC_COMPLEX_8 t1[t1_dim]; /* was [256][256] */
-#pragma GCC diagnostic pop
-
-      /* Empty c first.  */
-      for (j=1; j<=n; j++)
-       for (i=1; i<=m; i++)
-         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+      t1 = malloc (t1_dim * sizeof(GFC_COMPLEX_8));
 
       /* Start turning the crank. */
       i1 = n;
@@ -1641,6 +1685,7 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
                }
            }
        }
+      free(t1);
       return;
     }
   else if (rxstride == 1 && aystride == 1 && bxstride == 1)
@@ -1737,6 +1782,24 @@ matmul_c8_avx512f (gfc_array_c8 * const restrict retarray,
 
 #endif  /* HAVE_AVX512F */
 
+/* AMD-specifix funtions with AVX128 and FMA3/FMA4.  */
+
+#if defined(HAVE_AVX) && defined(HAVE_FMA3) && defined(HAVE_AVX128)
+void
+matmul_c8_avx128_fma3 (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) __attribute__((__target__("avx,fma")));
+internal_proto(matmul_c8_avx128_fma3);
+#endif
+
+#if defined(HAVE_AVX) && defined(HAVE_FMA4) && defined(HAVE_AVX128)
+void
+matmul_c8_avx128_fma4 (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) __attribute__((__target__("avx,fma4")));
+internal_proto(matmul_c8_avx128_fma4);
+#endif
+
 /* Function to fall back to if there is no special processor-specific version.  */
 static void
 matmul_c8_vanilla (gfc_array_c8 * const restrict retarray, 
@@ -1799,8 +1862,8 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -1808,8 +1871,8 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else
@@ -1817,17 +1880,15 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 1:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
 
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 2:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 2 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
     }
@@ -1868,7 +1929,9 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
   if (count != GFC_DESCRIPTOR_EXTENT(b,0))
     {
       if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
-       runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
+       runtime_error ("Incorrect extent in argument B in MATMUL intrinsic "
+                      "in dimension 1: is %ld, should be %ld",
+                      (long int) GFC_DESCRIPTOR_EXTENT(b,0), (long int) count);
     }
 
   if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -1877,9 +1940,9 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
       bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
 
       /* bystride should never be used for 1-dimensional b.
-        in case it is we want it to cause a segfault, rather than
-        an incorrect result. */
-      bystride = 0xDEADBEEF;
+         The value is only used for calculation of the
+         memory by the buffer.  */
+      bystride = 256;
       ycount = 1;
     }
   else
@@ -1913,7 +1976,18 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
       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,
+         const char *transa, *transb;
+         if (try_blas & 2)
+           transa = "C";
+         else
+           transa = axstride == 1 ? "N" : "T";
+
+         if (try_blas & 4)
+           transb = "C";
+         else
+           transb = bxstride == 1 ? "N" : "T";
+
+         gemm (transa, transb , &m,
                &n, &k, &one, abase, &lda, bbase, &ldb, &zero, dest,
                &ldc, 1, 1);
          return;
@@ -1945,6 +2019,7 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
                 f13, f14, f23, f24, f33, f34, f43, f44;
       index_type i, j, l, ii, jj, ll;
       index_type isec, jsec, lsec, uisec, ujsec, ulsec;
+      GFC_COMPLEX_8 *t1;
 
       a = abase;
       b = bbase;
@@ -1961,25 +2036,27 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
       b_offset = 1 + b_dim1;
       b -= b_offset;
 
+      /* Empty c first.  */
+      for (j=1; j<=n; j++)
+       for (i=1; i<=m; i++)
+         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+
       /* Early exit if possible */
       if (m == 0 || n == 0 || k == 0)
        return;
 
       /* Adjust size of t1 to what is needed.  */
-      index_type t1_dim;
-      t1_dim = (a_dim1-1) * 256 + b_dim1;
+      index_type t1_dim, a_sz;
+      if (aystride == 1)
+        a_sz = rystride;
+      else
+        a_sz = a_dim1;
+
+      t1_dim = a_sz * 256 + b_dim1;
       if (t1_dim > 65536)
        t1_dim = 65536;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wvla"
-      GFC_COMPLEX_8 t1[t1_dim]; /* was [256][256] */
-#pragma GCC diagnostic pop
-
-      /* Empty c first.  */
-      for (j=1; j<=n; j++)
-       for (i=1; i<=m; i++)
-         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+      t1 = malloc (t1_dim * sizeof(GFC_COMPLEX_8));
 
       /* Start turning the crank. */
       i1 = n;
@@ -2190,6 +2267,7 @@ matmul_c8_vanilla (gfc_array_c8 * const restrict retarray,
                }
            }
        }
+      free(t1);
       return;
     }
   else if (rxstride == 1 && aystride == 1 && bxstride == 1)
@@ -2336,6 +2414,26 @@ void matmul_c8 (gfc_array_c8 * const restrict retarray,
            }
 #endif  /* HAVE_AVX */
         }
+    else if (__cpu_model.__cpu_vendor == VENDOR_AMD)
+      {
+#if defined(HAVE_AVX) && defined(HAVE_FMA3) && defined(HAVE_AVX128)
+        if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
+           && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA)))
+         {
+            matmul_fn = matmul_c8_avx128_fma3;
+           goto store;
+         }
+#endif
+#if defined(HAVE_AVX) && defined(HAVE_FMA4) && defined(HAVE_AVX128)
+        if ((__cpu_model.__cpu_features[0] & (1 << FEATURE_AVX))
+            && (__cpu_model.__cpu_features[0] & (1 << FEATURE_FMA4)))
+         {
+            matmul_fn = matmul_c8_avx128_fma4;
+           goto store;
+         }
+#endif
+
+      }
    store:
       __atomic_store_n (&matmul_p, matmul_fn, __ATOMIC_RELAXED);
    }
@@ -2406,8 +2504,8 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -2415,8 +2513,8 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic: is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
       else
@@ -2424,17 +2522,15 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
          arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 1:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 1 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
 
          arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
          ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
          if (arg_extent != ret_extent)
-           runtime_error ("Incorrect extent in return array in"
-                          " MATMUL intrinsic for dimension 2:"
-                          " is %ld, should be %ld",
+           runtime_error ("Array bound mismatch for dimension 2 of "
+                          "array (%ld/%ld) ",
                           (long int) ret_extent, (long int) arg_extent);
        }
     }
@@ -2475,7 +2571,9 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
   if (count != GFC_DESCRIPTOR_EXTENT(b,0))
     {
       if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
-       runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
+       runtime_error ("Incorrect extent in argument B in MATMUL intrinsic "
+                      "in dimension 1: is %ld, should be %ld",
+                      (long int) GFC_DESCRIPTOR_EXTENT(b,0), (long int) count);
     }
 
   if (GFC_DESCRIPTOR_RANK (b) == 1)
@@ -2484,9 +2582,9 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
       bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
 
       /* bystride should never be used for 1-dimensional b.
-        in case it is we want it to cause a segfault, rather than
-        an incorrect result. */
-      bystride = 0xDEADBEEF;
+         The value is only used for calculation of the
+         memory by the buffer.  */
+      bystride = 256;
       ycount = 1;
     }
   else
@@ -2520,7 +2618,18 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
       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,
+         const char *transa, *transb;
+         if (try_blas & 2)
+           transa = "C";
+         else
+           transa = axstride == 1 ? "N" : "T";
+
+         if (try_blas & 4)
+           transb = "C";
+         else
+           transb = bxstride == 1 ? "N" : "T";
+
+         gemm (transa, transb , &m,
                &n, &k, &one, abase, &lda, bbase, &ldb, &zero, dest,
                &ldc, 1, 1);
          return;
@@ -2552,6 +2661,7 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
                 f13, f14, f23, f24, f33, f34, f43, f44;
       index_type i, j, l, ii, jj, ll;
       index_type isec, jsec, lsec, uisec, ujsec, ulsec;
+      GFC_COMPLEX_8 *t1;
 
       a = abase;
       b = bbase;
@@ -2568,25 +2678,27 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
       b_offset = 1 + b_dim1;
       b -= b_offset;
 
+      /* Empty c first.  */
+      for (j=1; j<=n; j++)
+       for (i=1; i<=m; i++)
+         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+
       /* Early exit if possible */
       if (m == 0 || n == 0 || k == 0)
        return;
 
       /* Adjust size of t1 to what is needed.  */
-      index_type t1_dim;
-      t1_dim = (a_dim1-1) * 256 + b_dim1;
+      index_type t1_dim, a_sz;
+      if (aystride == 1)
+        a_sz = rystride;
+      else
+        a_sz = a_dim1;
+
+      t1_dim = a_sz * 256 + b_dim1;
       if (t1_dim > 65536)
        t1_dim = 65536;
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wvla"
-      GFC_COMPLEX_8 t1[t1_dim]; /* was [256][256] */
-#pragma GCC diagnostic pop
-
-      /* Empty c first.  */
-      for (j=1; j<=n; j++)
-       for (i=1; i<=m; i++)
-         c[i + j * c_dim1] = (GFC_COMPLEX_8)0;
+      t1 = malloc (t1_dim * sizeof(GFC_COMPLEX_8));
 
       /* Start turning the crank. */
       i1 = n;
@@ -2797,6 +2909,7 @@ matmul_c8 (gfc_array_c8 * const restrict retarray,
                }
            }
        }
+      free(t1);
       return;
     }
   else if (rxstride == 1 && aystride == 1 && bxstride == 1)