From 097bb2ca8d8b0413d19a8b9be042926c42c4049c Mon Sep 17 00:00:00 2001
From: Angelika Schwarz <17718454+angsch@users.noreply.github.com>
Date: Sun, 23 Jul 2023 22:24:53 +0200
Subject: [PATCH] Fix loop over block columns of C

If the column count of C exceeds 1000, the code
encounters a segmentation fault. The block size
is defined #define nb 1000 and used to allocate
the buffers. However, n instead of nb is
passed.

* add the missing loop over block columns of C
* adjust the leading dimensions to select
  whatever the row count is
---
 src/HowToOptimizeGemm/parameters.h |  6 +++---
 src/MMult_4x4_15.c                 | 20 +++++++++++---------
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/src/HowToOptimizeGemm/parameters.h b/src/HowToOptimizeGemm/parameters.h
index f26f626..c139c6f 100644
--- a/src/HowToOptimizeGemm/parameters.h
+++ b/src/HowToOptimizeGemm/parameters.h
@@ -35,6 +35,6 @@ leading dimension of the array that stores matrix X.  If LDX=-1
 then the leading dimension is set to the row dimension of matrix X.
 */
 
-#define LDA 1000
-#define LDB 1000
-#define LDC 1000
+#define LDA -1
+#define LDB -1
+#define LDC -1
diff --git a/src/MMult_4x4_15.c b/src/MMult_4x4_15.c
index 7336df8..b183e74 100644
--- a/src/MMult_4x4_15.c
+++ b/src/MMult_4x4_15.c
@@ -22,15 +22,17 @@ void MY_MMult( int m, int n, int k, double *a, int lda,
                                     double *b, int ldb,
                                     double *c, int ldc )
 {
-  int i, p, pb, ib;
-
-  /* This time, we compute a mc x n block of C by a call to the InnerKernel */
-
-  for ( p=0; p<k; p+=kc ){
-    pb = min( k-p, kc );
-    for ( i=0; i<m; i+=mc ){
-      ib = min( m-i, mc );
-      InnerKernel( ib, n, pb, &A( i,p ), lda, &B(p, 0 ), ldb, &C( i,0 ), ldc, i==0 );
+  int i, p, pb, ib, j, jb;
+
+  /* This time, we compute a mc x nb block of C by a call to the InnerKernel */
+  for ( j=0; j<n; j+=nb ){
+    jb = min( n-j, nb );
+    for ( p=0; p<k; p+=kc ){
+      pb = min( k-p, kc );
+      for ( i=0; i<m; i+=mc ){
+        ib = min( m-i, mc );
+        InnerKernel( ib, jb, pb, &A( i,p ), lda, &B(p,j ), ldb, &C( i,j ), ldc, i==0 );
+      }
     }
   }
 }