Merge pull request #5990 from jjerphan/wasm128-gemv

Add WASM SIMD128 SGEMV/DGEMV kernels
2026-08-22 12:10:59 +02:00 · 2026-08-21 16:18:59 +02:00 · 2026-08-21 06:25:29 +02:00 · 2026-08-21 01:09:35 +02:00 · 2026-08-21 00:57:38 +02:00 · 2026-08-20 22:32:08 +02:00
171 changed files with 20205 additions and 1853 deletions
@@ -593,7 +593,7 @@ jobs:
      - name: Build OpenBLAS
        run: |
          cd build
-          cmake --build . --target dgemm_thread_safety dgemm_thread_safety_mixed dgemv_thread_safety
+          cmake --build . --target dgemm_thread_safety dgemm_thread_safety_mixed dgemm_thread_safety_shutdown dgemv_thread_safety

      - name: Show ccache status
        continue-on-error: true
@@ -611,7 +611,7 @@ jobs:
        run: |
          cd build
          export PATH="$PWD/lib:$PATH"
-          OPENBLAS_NUM_THREADS=8 OMP_NUM_THREADS=16 ctest -R 'dgemm_thread_safety|dgemm_thread_safety_mixed|dgemv_thread_safety' --output-on-failure
+          OPENBLAS_NUM_THREADS=8 OMP_NUM_THREADS=16 ctest -R 'dgemm_thread_safety|dgemm_thread_safety_mixed|dgemm_thread_safety_shutdown|dgemv_thread_safety' --output-on-failure


  cross_build:
@@ -51,7 +51,7 @@ jobs:
      - name: Install APT deps
        run: |
          sudo apt-get update && \
-          sudo apt-get install autoconf automake autotools-dev ninja-build make ccache qemu-user-static \
+          sudo apt-get install autoconf automake autotools-dev ninja-build make ccache \
                               gcc-14-loongarch64-linux-gnu g++-14-loongarch64-linux-gnu gfortran-14-loongarch64-linux-gnu

      - name: Compilation cache
@@ -70,11 +70,18 @@ jobs:
          echo "compression = true" >> ~/.ccache/ccache.conf
          ccache -s

-      - name: Disable utest dsdot:dsdot_n_1
+      - name: Install QEMU 10.2.1
        run: |
-          echo -n > utest/test_dsdot.c
-          echo "Due to the current version of qemu causing utest cases to fail,"
-          echo "the utest dsdot:dsdot_n_1 have been temporarily disabled."
+          qemu_deb="${RUNNER_TEMP}/qemu-user_10.2.1+ds-1ubuntu3_amd64.deb"
+          qemu_root="${RUNNER_TEMP}/qemu-user-10.2.1"
+          wget --output-document="${qemu_deb}" \
+            https://archive.ubuntu.com/ubuntu/pool/universe/q/qemu/qemu-user_10.2.1+ds-1ubuntu3_amd64.deb
+          echo "e6bdaac070232a44b8657cad84ac0853d84b7167b55af82370af040291f55bd8  ${qemu_deb}" \
+            | sha256sum --check -
+          dpkg-deb --extract "${qemu_deb}" "${qemu_root}"
+          sudo install -m 0755 "${qemu_root}/usr/bin/qemu-loongarch64" \
+            /usr/local/bin/qemu-loongarch64-static
+          qemu-loongarch64-static --version

      - name: Build OpenBLAS
        run: |
@@ -92,12 +92,6 @@ jobs:
          echo "compression = true" >> ~/.ccache/ccache.conf
          ccache -s

-      - name: Disable utest dsdot:dsdot_n_1
-        run: |
-          echo -n > utest/test_dsdot.c
-          echo "Due to the qemu versions 7.2 causing utest cases to fail,"
-          echo "the utest dsdot:dsdot_n_1 have been temporarily disabled."
-
      - name: Build OpenBLAS
        run: make CC='ccache clang --target=loongarch64-linux-gnu --sysroot=/opt/loongson-gnu-toolchain-8.3-x86_64-loongarch64-linux-gnu-rc1.3/loongarch64-linux-gnu/sysroot/ -static' FC='ccache loongarch64-linux-gnu-gfortran -static' HOSTCC='ccache clang' CROSS_SUFFIX=llvm-  NO_SHARED=1 ${{ matrix.opts }} -j$(nproc)

@@ -56,7 +56,7 @@ jobs:
        with:
          repository: qemu/qemu
          path: qemu
-          ref: ae35f033b874c627d81d51070187fbf55f0bf1a7
+          ref: 3e0bcba1ca7d6607ca49a988d165f052a3a53323

      - name: build qemu
        run: |
@@ -283,4 +283,8 @@ hheei <hheei@users.noreply.github.com>
 * Aadityansha Verma <https://github.com/aadityansha06>
  * [2026-07-14] Add independent transpose support for C in GEADD (sgeadd/dgeadd/cgeadd/zgeadd).

+* Vincent Lovero <https://github.com/vlovero>
+  * [2026-08-11] ARM v9.2 SME GEMM kernels for Apple M

+* Hugo Meiland <hugo@meiland.nl>
+  * [2026-08-09] Add Cortex-A72 DGEMM 6x8 microkernel and blocking
@@ -158,6 +158,9 @@ x86_64:

 wasm:
 - added optimized kernels for STRSM and DTRSM
+ - added a 4x4 WASM SIMD128 GEMM microkernel for SGEMM/DGEMM with matching
+   4-wide packing and TRMM, using relaxed SIMD madd in the inner loop
+ - gated leftover WASM TRSM call counters behind OPENBLAS_WASM_TRSM_PROFILE

 ====================================================================
 Version 0.3.32
@@ -59,7 +59,7 @@ endif
 endif

 ifeq ($(CORE), ARMV9SME)
-CCOMMON_OPT += -march=armv9-a+sve2+sme
+CCOMMON_OPT += -march=armv9-a+sve2+sme+sme-f64f64
 FCOMMON_OPT += -march=armv9-a+sve2
 ifdef OS_WINDOWS
 ifeq ($(C_COMPILER), CLANG)
@@ -310,7 +310,7 @@ endif

 ifeq ($(CORE), VORTEXM4)
 ifneq ($(C_COMPILER), GCC)
-CCOMMON_OPT += -march=armv8.4-a+sme
+CCOMMON_OPT += -march=armv8.4-a+sme+sme-f64f64
 #ifneq ($(APPLECLANG),1)
 #override LDFLAGS += -lclang_rt_builtins-aarch64
 #endif
@@ -285,7 +285,7 @@ GEMM_GEMV_FORWARD = 1
 SBGEMM_GEMV_FORWARD = 1
 BGEMM_GEMV_FORWARD = 1
 endif
-ifeq ($(ARCH), riscv)
+ifeq ($(ARCH), riscv64)
 GEMM_GEMV_FORWARD = 1
 endif
 ifeq ($(ARCH), power)
@@ -1 +1 @@
-CCOMMON_OPT += -msimd128
+CCOMMON_OPT += -msimd128 -mrelaxed-simd
@@ -165,7 +165,7 @@ jobs:

 - job: OSX_OpenMP
  pool:
-     vmImage: 'macOS-14'
+     vmImage: 'macOS-15'
  steps:   
  - script: |
      brew update
@@ -175,7 +175,7 @@ jobs:
     
 - job: OSX_GCC_Nothreads
  pool:
-     vmImage: 'macOS-14'
+     vmImage: 'macOS-15'
  steps:   
  - script: |
      brew update
@@ -190,6 +190,7 @@ jobs:
      make CC=gcc-15 FC=gfortran-15

 - job: OSX_LLVM_flangnew
+  timeoutInMinutes: 100
  pool:
     vmImage: 'macOS-latest'
  variables:
@@ -203,6 +204,7 @@ jobs:
      make TARGET=NEHALEM CC=/usr/local/opt/llvm/bin/clang FC=/usr/local/opt/flang/bin/flang NO_SHARED=1

 - job: OSX_OpenMP_Clang
+  timeoutInMinutes: 100
  pool:
     vmImage: 'macOS-latest'
  variables:
@@ -216,6 +218,7 @@ jobs:
      make TARGET=CORE2 USE_OPENMP=1 DYNAMIC_ARCH=1 DYNAMIC_LIST='NEHALEM HASWELL SKYLAKEX' CC=/usr/local/opt/llvm/bin/clang NOFORTRAN=1

 - job: OSX_OpenMP_Clang_cmake
+  timeoutInMinutes: 100
  pool:
     vmImage: 'macOS-latest'
  variables:
@@ -233,8 +236,9 @@ jobs:
      ctest
      
 - job: OSX_dynarch_cmake
+  timeoutInMinutes: 100
  pool:
-     vmImage: 'macOS-14'
+     vmImage: 'macOS-15'
  variables:
     LD_LIBRARY_PATH: /usr/local/opt/llvm/lib
     LIBRARY_PATH: /usr/local/opt/llvm/lib
@@ -283,7 +287,7 @@ jobs:
 
 - job: OSX_NDK_ARMV7
  pool:
-     vmImage: 'macOS-14'
+     vmImage: 'macOS-15'
  steps:   
  - script: | 
      brew update
@@ -293,19 +297,19 @@ jobs:

 - job: OSX_IOS_ARMV8
  pool:
-     vmImage: 'macOS-14'
+     vmImage: 'macOS-15'
  variables:
-     CC: /Applications/Xcode_16.2.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang
-     CFLAGS: -O2 -Wno-macro-redefined -isysroot /Applications/Xcode_16.2.app/Contents/Developer/Platforms/iPhoneOS.platform/Developer/SDKs/iPhoneOS18.2.sdk -arch arm64 -miphoneos-version-min=10.0
+     CC: /Applications/Xcode_16.4.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang
+     CFLAGS: -O2 -Wno-macro-redefined -isysroot /Applications/Xcode_16.4.app/Contents/Developer/Platforms/iPhoneOS.platform/Developer/SDKs/iPhoneOS18.5.sdk -arch arm64 -miphoneos-version-min=10.0
  steps:
  - script: |
-     ls /Applications/Xcode_16.2.app/Contents/Developer/Platforms/
-     ls /Applications/Xcode_16.2.app/Contents/Developer/Platforms/iPhoneOS.platform/Developer/SDKs
+     ls /Applications/Xcode_16.4.app/Contents/Developer/Platforms/
+     ls /Applications/Xcode_16.4.app/Contents/Developer/Platforms/iPhoneOS.platform/Developer/SDKs
     make TARGET=ARMV8 DYNAMIC_ARCH=1 NUM_THREADS=32 HOSTCC=clang NOFORTRAN=1

 - job: OSX_IOS_ARMV7
  pool:
-     vmImage: 'macOS-14'
+     vmImage: 'macOS-15'
  variables:
     CC: /Applications/Xcode_16.2.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang
     CFLAGS: -O2 -mno-thumb -Wno-macro-redefined -isysroot /Applications/Xcode_16.2.app/Contents/Developer/Platforms/iPhoneOS.platform/Developer/SDKs/iPhoneOS18.2.sdk -arch armv7 -miphoneos-version-min=5.1
@@ -315,15 +319,15 @@ jobs:

 - job: OSX_xbuild_DYNAMIC_ARM64
  pool:
-     vmImage: 'macOS-14'
+     vmImage: 'macOS-15'
  variables:
-     CC: /Applications/Xcode_16.2.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang
-     CFLAGS: -O2 -Wno-macro-redefined -isysroot /Applications/Xcode_16.2.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX15.2.sdk -arch arm64
+     CC: /Applications/Xcode_16.4.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang
+     CFLAGS: -O2 -Wno-macro-redefined -isysroot /Applications/Xcode_16.4.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX15.5.sdk -arch arm64
  steps:
  - script: |
-     ls /Applications/Xcode_16.2.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs
-     /Applications/Xcode_12.2.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang -arch arm64 --print-supported-cpus
-     /Applications/Xcode_16.2.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang --version
+     ls /Applications/Xcode_16.4.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs
+     /Applications/Xcode_16.4.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang -arch arm64 --print-supported-cpus
+     /Applications/Xcode_16.4.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clang --version
     make TARGET=ARMV8 DYNAMIC_ARCH=1 NUM_THREADS=32 HOSTCC=clang NOFORTRAN=1

 - job: ALPINE_MUSL
@@ -95,10 +95,15 @@ else
 GOTO_HFLOAT_TARGETS=
 endif

+ifeq ($(USE_OPENMP), 1)
+SMALLSCALING=smallscaling
+endif
+
 ifeq ($(OSNAME), WINNT)

 goto :: slinpack.goto dlinpack.goto clinpack.goto zlinpack.goto \
       scholesky.goto dcholesky.goto ccholesky.goto zcholesky.goto \
+       cblas_sgemm.goto cblas_dgemm.goto cblas_cgemm.goto cblas_zgemm.goto \
       sgemm.goto dgemm.goto cgemm.goto zgemm.goto \
       strmm.goto dtrmm.goto ctrmm.goto ztrmm.goto \
       strsm.goto dtrsm.goto ctrsm.goto ztrsm.goto \
@@ -268,6 +273,7 @@ mkl :: slinpack.mkl dlinpack.mkl clinpack.mkl zlinpack.mkl \
 else

 goto :: sgemm.goto dgemm.goto cgemm.goto zgemm.goto \
+       cblas_sgemm.goto cblas_dgemm.goto cblas_cgemm.goto cblas_zgemm.goto \
       strmm.goto dtrmm.goto ctrmm.goto ztrmm.goto \
       strsm.goto dtrsm.goto ctrsm.goto ztrsm.goto \
       sspr.goto dspr.goto \
@@ -301,7 +307,7 @@ goto :: sgemm.goto dgemm.goto cgemm.goto zgemm.goto \
       stpsv.goto dtpsv.goto ctpsv.goto ztpsv.goto \
       strsv.goto dtrsv.goto ctrsv.goto ztrsv.goto \
       ssymm.goto dsymm.goto csymm.goto zsymm.goto \
-       smallscaling \
+       $(SMALLSCALING) \
       isamax.goto idamax.goto icamax.goto izamax.goto \
       ismax.goto idmax.goto \
       isamin.goto idamin.goto icamin.goto izamin.goto \
@@ -681,6 +687,18 @@ endif
 sgemm.goto : sgemm.$(SUFFIX) ../$(LIBNAME)
 	$(CC) $(CFLAGS) -o $(@F) $^ $(CEXTRALIB) $(EXTRALIB) $(FEXTRALIB) -lm

+cblas_sgemm.goto : cblas_sgemm.$(SUFFIX) ../$(LIBNAME)
+	$(CC) $(CFLAGS) -o $(@F) $^ $(CEXTRALIB) $(EXTRALIB) $(FEXTRALIB) -lm
+
+cblas_dgemm.goto : cblas_dgemm.$(SUFFIX) ../$(LIBNAME)
+	$(CC) $(CFLAGS) -o $(@F) $^ $(CEXTRALIB) $(EXTRALIB) $(FEXTRALIB) -lm
+
+cblas_cgemm.goto : cblas_cgemm.$(SUFFIX) ../$(LIBNAME)
+	$(CC) $(CFLAGS) -o $(@F) $^ $(CEXTRALIB) $(EXTRALIB) $(FEXTRALIB) -lm
+
+cblas_zgemm.goto : cblas_zgemm.$(SUFFIX) ../$(LIBNAME)
+	$(CC) $(CFLAGS) -o $(@F) $^ $(CEXTRALIB) $(EXTRALIB) $(FEXTRALIB) -lm
+
 sgemm.acml : sgemm.$(SUFFIX)
 	-$(CC) $(CFLAGS) -o $(@F) $^ $(LIBACML) $(CEXTRALIB) $(EXTRALIB) $(FEXTRALIB)

@@ -3027,6 +3045,18 @@ cgemm.$(SUFFIX) : gemm.c
 zgemm.$(SUFFIX) : gemm.c
 	$(CC) $(CFLAGS) -c -DCOMPLEX -DDOUBLE -o $(@F) $^

+cblas_sgemm.$(SUFFIX) : cblasgemm.c
+	$(CC) $(CFLAGS) -c -UCOMPLEX -UDOUBLE -o $(@F) $^
+
+cblas_dgemm.$(SUFFIX) : cblasgemm.c
+	$(CC) $(CFLAGS) -c -UCOMPLEX -DDOUBLE -o $(@F) $^
+
+cblas_cgemm.$(SUFFIX) : cblasgemm.c
+	$(CC) $(CFLAGS) -c -DCOMPLEX -UDOUBLE -o $(@F) $^
+
+cblas_zgemm.$(SUFFIX) : cblasgemm.c
+	$(CC) $(CFLAGS) -c -DCOMPLEX -DDOUBLE -o $(@F) $^
+
 ssymm.$(SUFFIX) : symm.c
 	$(CC) $(CFLAGS) -c -UCOMPLEX -UDOUBLE -o $(@F) $^

@@ -3533,7 +3563,11 @@ zomatcopy.$(SUFFIX) : omatcopy.c


 smallscaling: smallscaling.c ../$(LIBNAME)
+ifeq ($(C_COMPILER), GCC)
 	$(CC) $(CFLAGS) -o $(@F) $^ $(EXTRALIB) -fopenmp -lm -lpthread
+else
+	$(CC) $(CFLAGS) -o $(@F) $^ $(EXTRALIB) -quak -openmp -lm -lpthread
+endif

 clean ::
 	@rm -f *.goto *.mkl *.acml *.atlas *.veclib *.essl smallscaling
@@ -0,0 +1,233 @@
+/***************************************************************************
+Copyright (c) 2014, The OpenBLAS Project
+All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in
+the documentation and/or other materials provided with the
+distribution.
+3. Neither the name of the OpenBLAS project nor the names of
+its contributors may be used to endorse or promote products
+derived from this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*****************************************************************************/
+
+#include "bench.h"
+#include "cblas.h"
+#undef GEMM
+
+#ifndef COMPLEX
+
+#ifdef DOUBLE
+#define GEMM   cblas_dgemm
+#elif defined(BFLOAT16) && defined(BGEMM)
+#define GEMM   cblas_bgemm
+#elif defined(BFLOAT16)
+#define GEMM   cblas_sbgemm
+#undef IFLOAT
+#define IFLOAT bfloat16
+#elif defined(HFLOAT16)
+#define GEMM   cblas_shgemm
+#undef IFLOAT
+#define IFLOAT hfloat16
+#else
+#define GEMM   cblas_sgemm
+#undef IFLOAT
+#define IFLOAT float
+#endif
+
+#else
+
+#ifdef DOUBLE
+#define GEMM   cblas_zgemm
+#else
+#define GEMM   cblas_cgemm
+#endif
+
+#endif
+
+int main(int argc, char *argv[]){
+
+  IFLOAT *a, *b;
+  //IFLOAT *aa, *bb;
+  FLOAT *c;
+  //FLOAT *cc;
+#ifdef BGEMM
+  blasint one=1;
+  blasint two=2;
+  float alpha_in[] = {1.0, 0.0};
+  float beta_in[] = {0.0, 0.0};
+  FLOAT alpha[2], beta[2];
+  sbstobf16_(&two, alpha_in, &one, alpha, &one);
+  sbstobf16_(&two, beta_in, &one, beta, &one);
+#else
+#ifdef COMPLEX
+  FLOAT alpha[] = {1.0, 0.0};
+  FLOAT beta [] = {0.0, 0.0};
+#else
+  FLOAT alpha = 1.0;
+  FLOAT beta = 0.0;
+#endif
+#endif
+  CBLAS_TRANSPOSE transa = CblasNoTrans;
+  CBLAS_TRANSPOSE transb = CblasNoTrans;
+ char transac, transbc;
+  blasint m, n, k, i, j, lda, ldb, ldc;
+  int loops = 1;
+  int has_param_m = 0;
+  int has_param_n = 0;
+  int has_param_k = 0;
+  int has_param_lda = 0;
+  int has_param_ldb = 0;
+  char *p;
+//blasint sme=0;
+  int from =   1;
+  int to   = 200;
+  int step =   1;
+
+  double time1, timeg;
+
+  argc--;argv++;
+
+  if (argc > 0) { from = atol(*argv);            argc--; argv++; }
+  if (argc > 0) { to   = MAX(atol(*argv), from); argc--; argv++; }
+  if (argc > 0) { step = atol(*argv);            argc--; argv++; }
+
+  if ((p = getenv("OPENBLAS_TRANS"))) {
+  transa=(*p=='N') ? CblasNoTrans : CblasTrans;
+  transb=(*p=='N') ? CblasNoTrans : CblasTrans;
+  }
+  if ((p = getenv("OPENBLAS_TRANSA"))) {
+  transa=(*p=='N') ? CblasNoTrans : CblasTrans;
+  }
+  if ((p = getenv("OPENBLAS_TRANSB"))) {
+  transb=(*p=='N') ? CblasNoTrans : CblasTrans;
+  }
+
+  transac=(transa==CblasNoTrans) ? 'N' : 'T';
+  transbc=(transb==CblasNoTrans) ? 'N' : 'T';
+  fprintf(stderr, "From : %3d  To : %3d Step=%d : Transa=%c : Transb=%c\n", from, to, step, transac, transbc);
+
+  p = getenv("OPENBLAS_LOOPS");
+  if ( p != NULL ) {
+    loops = atoi(p);
+  }
+
+  if ((p = getenv("OPENBLAS_PARAM_M"))) {
+    m = atoi(p);
+    has_param_m=1;
+  } else {
+    m = to;
+  }
+  if ((p = getenv("OPENBLAS_PARAM_N"))) {
+    n = atoi(p);
+    has_param_n=1;
+  } else {
+    n = to;
+  }
+  if ((p = getenv("OPENBLAS_PARAM_K"))) {
+    k = atoi(p);
+    has_param_k=1;
+  } else {
+    k = to;
+  }
+  if ((p = getenv("OPENBLAS_PARAM_LDA"))) {
+    lda = atoi(p);
+    has_param_lda=1;
+  } 
+  if ((p = getenv("OPENBLAS_PARAM_LDB"))) {
+    ldb = atoi(p);
+    has_param_ldb=1;
+  } 
+
+  if (( a = (IFLOAT *)malloc(sizeof(IFLOAT) * m * k * COMPSIZE)) == NULL) {
+    fprintf(stderr,"Out of Memory!!\n");exit(1);
+  }
+  if (( b = (IFLOAT *)malloc(sizeof(IFLOAT) * k * n * COMPSIZE)) == NULL) {
+    fprintf(stderr,"Out of Memory!!\n");exit(1);
+  }
+  if (( c = (FLOAT *)malloc(sizeof(FLOAT) * m * n * COMPSIZE)) == NULL) {
+    fprintf(stderr,"Out of Memory!!\n");exit(1);
+  }
+  //if (( aa = (IFLOAT *)malloc(sizeof(IFLOAT) * m * k * COMPSIZE)) == NULL) {
+  //  fprintf(stderr,"Out of Memory!!\n");exit(1);
+  //}
+  //if (( bb = (IFLOAT *)malloc(sizeof(IFLOAT) * k * n * COMPSIZE)) == NULL) {
+  //  fprintf(stderr,"Out of Memory!!\n");exit(1);
+  //}
+  //if (( cc = (FLOAT *)malloc(sizeof(FLOAT) * m * n * COMPSIZE)) == NULL) {
+  //  fprintf(stderr,"Out of Memory!!\n");exit(1);
+  //}
+
+#ifdef __linux
+  srandom(getpid());
+#endif
+
+  for (i = 0; i < m * k * COMPSIZE; i++) {
+    a[i] = ((IFLOAT) rand() / (IFLOAT) RAND_MAX) - 0.5;
+  //  aa[i]=a[i];
+  }
+  for (i = 0; i < k * n * COMPSIZE; i++) {
+    b[i] = ((IFLOAT) rand() / (IFLOAT) RAND_MAX) - 0.5;
+  //  bb[i]=b[i];
+  }
+  for (i = 0; i < m * n * COMPSIZE; i++) {
+    c[i] = ((FLOAT) rand() / (FLOAT) RAND_MAX) - 0.5;
+  //  cc[i]=c[i];
+  }
+
+  fprintf(stderr, "          SIZE                   Flops             Time\n");
+
+  for (i = from; i <= to; i += step) {
+    
+    timeg=0;
+
+    if (!has_param_m) { m = i; }
+    if (!has_param_n) { n = i; }
+    if (!has_param_k) { k = i; }
+    
+    if (!has_param_lda) {
+      if (transa == CblasNoTrans) { lda = k; }
+      else { lda = m; }
+    }
+    if (!has_param_ldb) {
+      if (transb == CblasNoTrans) { ldb = n; }
+      else { ldb = k; }
+    }
+    ldc = n;
+
+    fprintf(stderr, " M=%4d, N=%4d, K=%4d : ", (int)m, (int)n, (int)k);
+    begin();
+
+    for (j=0; j<loops; j++) {
+      GEMM (CblasRowMajor,transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
+    }
+
+ // for (ii = 0; ii < m * n * COMPSIZE; ii++) if (fabsf(c[ii]-cc[ii])>1.5e-5){fprintf(stderr,"mismatch %d %f !=  %f: %g\n",ii,c[ii],cc[ii],fabsf(c[ii]-cc[ii]));}
+    end();
+    time1 = getsec();
+
+    timeg = time1/loops;
+    fprintf(stderr,
+	    " %10.2f MFlops %10.6f sec\n",
+	    COMPSIZE * COMPSIZE * 2. * (double)k * (double)m * (double)n / timeg * 1.e-6, time1);
+    
+  }
+
+  return 0;
+}
+
+// void main(int argc, char *argv[]) __attribute__((weak, alias("MAIN__")));
@@ -1,3 +1,4 @@
+//#pragma clang optimize off
 /***************************************************************************
 Copyright (c) 2014, The OpenBLAS Project
 All rights reserved.
@@ -45,6 +46,7 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define IFLOAT hfloat16
 #else
 #define GEMM   BLASFUNC(sgemm)
+#undef IFLOAT
 #define IFLOAT float
 #endif

@@ -186,7 +188,7 @@ int main(int argc, char *argv[]){

    timeg = time1/loops;
    fprintf(stderr,
-	    " %10.2f MFlops %10.6f sec\n",
+	    " %10.2lf MFlops %10.6f sec\n",
 	    COMPSIZE * COMPSIZE * 2. * (double)k * (double)m * (double)n / timeg * 1.e-6, time1);
    
  }
@@ -315,7 +315,7 @@ if (${CORE} STREQUAL ARMV9SME)
    if (${CMAKE_C_COMPILER_ID} STREQUAL "NVHPC" AND NOT NO_SVE)
      set (CCOMMON_OPT  "${CCOMMON_OPT} -tp=host")
    else ()
-      set (CCOMMON_OPT "${CCOMMON_OPT} -march=armv9-a+sme")
+      set (CCOMMON_OPT "${CCOMMON_OPT} -march=armv9-a+sme+sme-f64f64")
      if (CMAKE_SYSTEM_NAME STREQUAL "Windows" AND CMAKE_C_COMPILER_ID MATCHES "Clang")
        set (CCOMMON_OPT "${CCOMMON_OPT} -mllvm --aarch64-stack-hazard-size=0")
      endif ()
@@ -329,7 +329,7 @@ if (${CORE} STREQUAL VORTEXM4)
 	    set (CCOMMON_OPT  "${CCOMMON_OPT} -tp=host")
    else ()
 	  if (${CMAKE_C_COMPILER_ID} STREQUAL "AppleClang")
-        set (CCOMMON_OPT "${CCOMMON_OPT} -march=armv8.4-a+sme -mcpu=apple-m4")
+        set (CCOMMON_OPT "${CCOMMON_OPT} -march=armv8.4-a+sme+sme-f64f64 -mcpu=apple-m4")
 	  else ()
 	    set (CCOMMON_OPT "${CCOMMON_OPT} -march=armv8.4-a -mcpu=apple-m4")
      endif ()
@@ -1313,7 +1313,7 @@ endif ()
    set(ZGEMM_UNROLL_M 2)
    set(ZGEMM_UNROLL_N 4)
    set(SYMV_P 16)
-  elseif ("${TCORE}" STREQUAL "ARMV8SVE" OR "${TCORE}" STREQUAL "CORTEXA510" OR "${TCORE}" STREQUAL "CORTEXX2" OR "${TCORE}" STREQUAL "ARMV9")
+  elseif ("${TCORE}" STREQUAL "ARMV8SVE" OR "${TCORE}" STREQUAL "CORTEXA510" OR "${TCORE}" STREQUAL "CORTEXX2" OR "${TCORE}" STREQUAL "ARMV9SME")
    file(APPEND ${TARGET_CONF_TEMP}
      "#define L1_DATA_SIZE\t32768\n"
      "#define L1_DATA_LINESIZE\t64\n"
@@ -147,6 +147,9 @@ if (POWER)
    set(HAVE_GAS 0)
  elseif (CMAKE_ASM_COMPILER_ID STREQUAL "Clang")
    set(CCOMMON_OPT "${CCOMMON_OPT} -fno-integrated-as")
+	if (${CMAKE_SYSTEM_NAME} STREQUAL "FreeBSD")
+      set(CCOMMON_OPT "${CCOMMON_OPT} -Wa,-mpwr10")
+    endif ()
    set(HAVE_GAS 0)
  endif ()
  set(GETARCH_FLAGS "${GETARCH_FLAGS} -DHAVE_GAS=${HAVE_GAS}")
@@ -366,13 +369,13 @@ if (${TARGET} STREQUAL NEOVERSEV1)
    endif()
  endif()
  if (${TARGET} STREQUAL ARMV9SME)
-      set (KERNEL_DEFINITIONS  "${KERNEL_DEFINITIONS} -march=armv9-a+sme -O3")
+      set (KERNEL_DEFINITIONS  "${KERNEL_DEFINITIONS} -march=armv9-a+sme+sme-f64f64 -O3")
 	  if (${CMAKE_SYSTEM_NAME} STREQUAL Windows AND ${CMAKE_C_COMPILER_ID} MATCHES "Clang")
 	  set (KERNEL_DEFINITIONS  "${KERNEL_DEFINITIONS} -mllvm --aarch64-stack-hazard-size=0")
 	  endif()
  endif()
  if (${TARGET} STREQUAL VORTEXM4)
-      set (KERNEL_DEFINITIONS  "${KERNEL_DEFINITIONS} -march=armv8.4-a+sme -O3")
+      set (KERNEL_DEFINITIONS  "${KERNEL_DEFINITIONS} -march=armv8.4-a+sme+sme-f64f64 -O3")
 	  if (${CMAKE_SYSTEM_NAME} STREQUAL Windows AND ${CMAKE_C_COMPILER_ID} MATCHES "Clang")
 	  set (KERNEL_DEFINITIONS  "${KERNEL_DEFINITIONS} -mllvm --aarch64-stack-hazard-size=0")
 	  endif()
@@ -145,6 +145,16 @@ static __inline int blas_quickdivide(unsigned int x, unsigned int y){
 	.ent REALNAME; \
 REALNAME:

+#define PROLOGUE_EXPORT \
+	.arch ev6; \
+	.set noat; \
+	.set noreorder; \
+.text; \
+	.align 5; \
+	.globl REALNAME; \
+	.ent REALNAME; \
+REALNAME:
+
 #ifdef PROFILE
 #define PROFCODE \
 	ldgp	$gp, 0($27); \
@@ -119,6 +119,7 @@
 #endif

 #define	CGEMM_BETA		cgemm_beta
+#define	SME_CGEMM_KERNEL	sme_cgemm_kernel

 #define	CGEMM_KERNEL_N		cgemm_kernel_n
 #define	CGEMM_KERNEL_L		cgemm_kernel_l
@@ -326,6 +327,7 @@
 #define	CTRSM_ILTNCOPY		gotoblas -> ctrsm_iltncopy

 #define	CGEMM_BETA		gotoblas -> cgemm_beta
+#define	SME_CGEMM_KERNEL	gotoblas -> sme_cgemm_kernel
 #define	CGEMM_KERNEL_N		gotoblas -> cgemm_kernel_n
 #define	CGEMM_KERNEL_L		gotoblas -> cgemm_kernel_l
 #define	CGEMM_KERNEL_R		gotoblas -> cgemm_kernel_r
@@ -114,6 +114,7 @@
 #define	DGEMM_BETA		dgemm_beta

 #define	DGEMM_KERNEL		dgemm_kernel
+#define SME_DGEMM_KERNEL	sme_dgemm_kernel

 #define	DTRMM_KERNEL_LN		dtrmm_kernel_LN
 #define	DTRMM_KERNEL_LT		dtrmm_kernel_LT
@@ -246,6 +247,7 @@

 #define	DGEMM_BETA		gotoblas -> dgemm_beta
 #define	DGEMM_KERNEL		gotoblas -> dgemm_kernel
+#define	SME_DGEMM_KERNEL	gotoblas -> sme_dgemm_kernel

 #define	DTRMM_KERNEL_LN		gotoblas -> dtrmm_kernel_LN
 #define	DTRMM_KERNEL_LT		gotoblas -> dtrmm_kernel_LT
@@ -371,6 +371,13 @@ static __inline int blas_quickdivide(unsigned int x, unsigned int y){
 	.proc REALNAME; \
 REALNAME:

+#define PROLOGUE_EXPORT \
+	.explicit; \
+	.text; \
+	.align 128; \
+	.global REALNAME; \
+	.proc REALNAME; \
+REALNAME:

 #ifdef PROFILE
 #define PROFCODE \
@@ -135,6 +135,31 @@ void ssyr2k_direct_alpha_betaLT(BLASLONG N, BLASLONG K,
 	float beta,
 	float * R, BLASLONG strideR);

+void sme_sgemm_kernel(char*, char*, BLASLONG M, BLASLONG N, BLASLONG K,
+	float * alpha,
+	float * A, BLASLONG ldA,
+	float * B, BLASLONG ldB,
+	float * beta,
+	float * C, BLASLONG ldC);
+void sme_dgemm_kernel(const char*, const char*, const BLASLONG M, const BLASLONG N, const BLASLONG K,
+	const double * alpha,
+	const double * A, const BLASLONG ldA,
+	const double * B, const BLASLONG ldB,
+	const double * beta,
+	double * C, const BLASLONG ldC);
+void sme_cgemm_kernel(const char*, const char*, const BLASLONG M, const BLASLONG N, const BLASLONG K,
+	const float alpha_r, const float alpha_i,
+	const float * A, const BLASLONG ldA,
+	const float * B, const BLASLONG ldB,
+	const float beta_r, const float beta_i,
+	float * C, const BLASLONG ldC);
+void sme_zgemm_kernel(const char*, const char*, const BLASLONG M, const BLASLONG N, const BLASLONG K,
+	const double alpha_r, const double alpha_i,
+	const double * A, const BLASLONG ldA,
+	const double * B, const BLASLONG ldB,
+	const double beta_r, const double beta_i,
+	double * C, const BLASLONG ldC);
+
 int sgemm_direct_performant(BLASLONG M, BLASLONG N, BLASLONG K);

 int shgemm_beta(BLASLONG, BLASLONG, BLASLONG, float,
@@ -276,6 +276,7 @@ int (*shgemv_t) (BLASLONG, BLASLONG, float, hfloat16 *, BLASLONG, hfloat16 *, BL
  void (*ssyr2k_direct_alpha_betaUT) (BLASLONG, BLASLONG, float, float *, BLASLONG, float *, BLASLONG, float, float *, BLASLONG);
  void (*ssyr2k_direct_alpha_betaLN) (BLASLONG, BLASLONG, float, float *, BLASLONG, float *, BLASLONG, float, float *, BLASLONG);
  void (*ssyr2k_direct_alpha_betaLT) (BLASLONG, BLASLONG, float, float *, BLASLONG, float *, BLASLONG, float, float *, BLASLONG);
+  void (*sme_sgemm_kernel) (char*, char*, BLASLONG, BLASLONG, BLASLONG, float*, float *, BLASLONG , float *, BLASLONG ,float*, float *, BLASLONG);
 #endif

  
@@ -401,6 +402,9 @@ int (*shgemv_t) (BLASLONG, BLASLONG, float, hfloat16 *, BLASLONG, hfloat16 *, BL
  int    (*dsymv_U) (BLASLONG, BLASLONG, double,  double  *, BLASLONG, double  *, BLASLONG, double  *, BLASLONG, double *);
 #endif
 #if (BUILD_DOUBLE==1) || (BUILD_COMPLEX16==1)
+#ifdef ARCH_ARM64
+  void (*sme_dgemm_kernel) (const char*, const char*, const BLASLONG, const BLASLONG, const BLASLONG, const double*, const double *, const BLASLONG , const double *, const BLASLONG ,const double*, double *, const BLASLONG);
+#endif
  int    (*dgemm_kernel   )(BLASLONG, BLASLONG, BLASLONG, double, double *, double *, double *, BLASLONG);
  int    (*dgemm_beta     )(BLASLONG, BLASLONG, BLASLONG, double, double *, BLASLONG, double *, BLASLONG, double  *, BLASLONG);

@@ -616,6 +620,9 @@ int (*shgemv_t) (BLASLONG, BLASLONG, float, hfloat16 *, BLASLONG, hfloat16 *, BL
  int    (*chemv_M) (BLASLONG, BLASLONG, float,  float, float  *, BLASLONG, float  *, BLASLONG, float  *, BLASLONG, float *);
  int    (*chemv_V) (BLASLONG, BLASLONG, float,  float, float  *, BLASLONG, float  *, BLASLONG, float  *, BLASLONG, float *);

+#ifdef ARCH_ARM64
+  void (*sme_cgemm_kernel) (const char*, const char*, const BLASLONG, const BLASLONG, const BLASLONG, const float, const float, const float *, const BLASLONG , const float *, const BLASLONG, const float, const float, float *, const BLASLONG);
+#endif
  int    (*cgemm_kernel_n )(BLASLONG, BLASLONG, BLASLONG, float, float, float *, float *, float *, BLASLONG);
  int    (*cgemm_kernel_l )(BLASLONG, BLASLONG, BLASLONG, float, float, float *, float *, float *, BLASLONG);
  int    (*cgemm_kernel_r )(BLASLONG, BLASLONG, BLASLONG, float, float, float *, float *, float *, BLASLONG);
@@ -826,6 +833,10 @@ int (*shgemv_t) (BLASLONG, BLASLONG, float, hfloat16 *, BLASLONG, hfloat16 *, BL
  int    (*zhemv_M) (BLASLONG, BLASLONG, double,  double, double  *, BLASLONG, double  *, BLASLONG, double  *, BLASLONG, double *);
  int    (*zhemv_V) (BLASLONG, BLASLONG, double,  double, double  *, BLASLONG, double  *, BLASLONG, double  *, BLASLONG, double *);

+#ifdef ARCH_ARM64
+  void (*sme_zgemm_kernel) (const char*, const char*, const BLASLONG, const BLASLONG, const BLASLONG, const double, const double, const double *, const BLASLONG , const double *, const BLASLONG, const double, const double, double *, const BLASLONG);
+#endif
+
  int    (*zgemm_kernel_n )(BLASLONG, BLASLONG, BLASLONG, double, double, double *, double *, double *, BLASLONG);
  int    (*zgemm_kernel_l )(BLASLONG, BLASLONG, BLASLONG, double, double, double *, double *, double *, BLASLONG);
  int    (*zgemm_kernel_r )(BLASLONG, BLASLONG, BLASLONG, double, double, double *, double *, double *, BLASLONG);
@@ -668,6 +668,17 @@ REALNAME:;\
 	.csect .text[PR],5;\
 .REALNAME:

+#define PROLOGUE_EXPORT \
+	.machine "any";\
+	.toc;\
+	.globl .REALNAME;\
+	.globl REALNAME;\
+	.csect REALNAME[DS],3;\
+REALNAME:;\
+	.long .REALNAME, TOC[tc0], 0;\
+	.csect .text[PR],5;\
+.REALNAME:
+
 #define EPILOGUE \
 _section_.text:;\
 	.csect .data[RW],4;\
@@ -686,6 +697,17 @@ REALNAME:;\
 	.csect .text[PR], 5;\
 .REALNAME:

+#define PROLOGUE_EXPORT \
+	.machine "any";\
+	.toc;\
+	.globl .REALNAME;\
+	.globl REALNAME;\
+	.csect REALNAME[DS],3;\
+REALNAME:;\
+	.llong .REALNAME, TOC[tc0], 0;\
+	.csect .text[PR], 5;\
+.REALNAME:
+
 #define EPILOGUE \
 _section_.text:;\
 	.csect .data[RW],4;\
@@ -705,6 +727,16 @@ _section_.text:;\
 	.text
 	.align 4
 	.globl REALNAME
+REALNAME:
+	.endmacro
+
+	.macro PROLOGUE_EXPORT
+	.section __TEXT,__text,regular,pure_instructions
+	.section __TEXT,__picsymbolstub1,symbol_stubs,pure_instructions,32
+	.machine ppc
+	.text
+	.align 4
+	.globl REALNAME
 REALNAME:
 	.endmacro
 #else
@@ -715,6 +747,16 @@ REALNAME:
 	.text
 	.align 4
 	.globl REALNAME
+REALNAME:
+	.endmacro
+
+	.macro PROLOGUE_EXPORT
+	.section __TEXT,__text,regular,pure_instructions
+	.section __TEXT,__picsymbolstub1,symbol_stubs,pure_instructions,32
+	.machine ppc64
+	.text
+	.align 4
+	.globl REALNAME
 REALNAME:
 	.endmacro
 #endif
@@ -76,6 +76,7 @@
 #define	SGEMM_ITCOPY		sgemm_itcopy
 #endif

+#define SME_SGEMM_KERNEL	sme_sgemm_kernel
 #define	STRMM_OUNUCOPY		strmm_ounucopy
 #define	STRMM_OUNNCOPY		strmm_ounncopy
 #define	STRMM_OUTUCOPY		strmm_outucopy
@@ -248,6 +249,7 @@
 #define  SSYR2K_DIRECT_ALPHA_BETA_UT     gotoblas -> ssyr2k_direct_alpha_betaUT
 #define  SSYR2K_DIRECT_ALPHA_BETA_LN     gotoblas -> ssyr2k_direct_alpha_betaLN
 #define  SSYR2K_DIRECT_ALPHA_BETA_LT     gotoblas -> ssyr2k_direct_alpha_betaLT
+#define SME_SGEMM_KERNEL	gotoblas -> sme_sgemm_kernel
 #endif

 #define	SGEMM_ONCOPY		gotoblas -> sgemm_oncopy
@@ -211,6 +211,14 @@ static __inline int blas_quickdivide(blasint x, blasint y){
 	.proc	07; \
 REALNAME:;

+#define PROLOGUE_EXPORT \
+	.section	".text"; \
+	.align 32; \
+	.global REALNAME;\
+	.type	REALNAME, #function; \
+	.proc	07; \
+REALNAME:;
+
 #if defined(__linux__) && defined(__ELF__)
 #define GNUSTACK .section .note.GNU-stack,"",@progbits
 #else
@@ -106,7 +106,15 @@ typedef struct blas_queue {
  blas_arg_t *args;
  void *range_m;
  void *range_n;
+
+  /*
+   * sa and sb are caller-owned inputs. worker_sb reports the sb workspace
+   * used by the threading backend for the most recent invocation.
+   * Keeping these roles separate prevents a reused queue from treating a
+   * released workspace as input to its next invocation.
+   */
  void *sa, *sb;
+  void *worker_sb;

  struct blas_queue *next;

@@ -148,6 +156,7 @@ static __inline int num_cpu_avail(int level) {
     `omp_get_max_threads()` below (which is to get a default
     in case the user hasn't made an explicit choice). */
  if (blas_is_num_threads_set_explicitly) {
+	  if (omp_in_parallel()) return 1;
    return blas_cpu_number;
  }

@@ -183,6 +192,7 @@ static __inline void blas_queue_init(blas_queue_t *queue){

  queue -> sa    = NULL;
  queue -> sb    = NULL;
+  queue -> worker_sb = NULL;
  queue-> next  = NULL;
 }

@@ -119,6 +119,7 @@
 #endif

 #define	ZGEMM_BETA		zgemm_beta
+#define	SME_ZGEMM_KERNEL	sme_zgemm_kernel

 #define	ZGEMM_KERNEL_N		zgemm_kernel_n
 #define	ZGEMM_KERNEL_L		zgemm_kernel_l
@@ -326,6 +327,7 @@
 #define	ZTRSM_ILTNCOPY		gotoblas -> ztrsm_iltncopy

 #define	ZGEMM_BETA		gotoblas -> zgemm_beta
+#define	SME_ZGEMM_KERNEL	gotoblas -> sme_zgemm_kernel
 #define	ZGEMM_KERNEL_N		gotoblas -> zgemm_kernel_n
 #define	ZGEMM_KERNEL_L		gotoblas -> zgemm_kernel_l
 #define	ZGEMM_KERNEL_R		gotoblas -> zgemm_kernel_r
@@ -19,6 +19,7 @@ endif()
 set(CPP_THREAD_SAFETY_DGEMM_ARGS "" CACHE STRING "Arguments passed to the DGEMM thread safety test")
 set(CPP_THREAD_SAFETY_DGEMM_MIXED_ARGS "" CACHE STRING "Arguments passed to the mixed DGEMM thread safety test")
 set(CPP_THREAD_SAFETY_DGEMV_ARGS "" CACHE STRING "Arguments passed to the DGEMV thread safety test")
+set(CPP_THREAD_SAFETY_SHUTDOWN_ARGS "" CACHE STRING "Arguments passed to the DGEMM shutdown safety test")

 if (CPP_THREAD_SAFETY_TEST)
  message(STATUS "building thread safety test")
@@ -29,6 +30,16 @@ if (CPP_THREAD_SAFETY_TEST)
  add_executable(dgemm_thread_safety_mixed dgemm_thread_safety_mixed.cpp)
  target_link_libraries(dgemm_thread_safety_mixed ${CPP_THREAD_SAFETY_LIBS})
  add_test(NAME dgemm_thread_safety_mixed COMMAND ${CMAKE_CURRENT_BINARY_DIR}/dgemm_thread_safety_mixed ${CPP_THREAD_SAFETY_DGEMM_MIXED_ARGS})
+
+  if (WIN32)
+    add_executable(dgemm_thread_safety_shutdown dgemm_thread_safety_shutdown.cpp)
+    target_link_libraries(dgemm_thread_safety_shutdown ${CPP_THREAD_SAFETY_LIBS})
+    add_test(NAME dgemm_thread_safety_shutdown COMMAND ${CMAKE_CURRENT_BINARY_DIR}/dgemm_thread_safety_shutdown ${CPP_THREAD_SAFETY_SHUTDOWN_ARGS})
+    set_tests_properties(dgemm_thread_safety_shutdown PROPERTIES TIMEOUT 900)
+  endif()
+  if (USE_THREAD AND (USE_OPENMP OR (NOT WIN32 AND NOT CYGWIN)))
+    add_test(NAME dgemm_thread_safety_mixed_callback COMMAND ${CMAKE_CURRENT_BINARY_DIR}/dgemm_thread_safety_mixed ${CPP_THREAD_SAFETY_DGEMM_MIXED_ARGS} --callback)
+  endif()
 endif()


@@ -3,6 +3,26 @@ include $(TOPDIR)/Makefile.system

 all :: dgemv_tester dgemm_tester dgemm_mixed_tester

+ifeq ($(OSNAME), WINNT)
+all :: dgemm_shutdown_tester
+endif
+
+CALLBACK_THREAD_TEST =
+ifeq ($(SMP),1)
+ifeq ($(USE_OPENMP),1)
+CALLBACK_THREAD_TEST = 1
+else
+ifneq ($(OSNAME),$(filter $(OSNAME),WINNT CYGWIN_NT))
+CALLBACK_THREAD_TEST = 1
+endif
+endif
+endif
+
+ifeq ($(CALLBACK_THREAD_TEST),1)
+.PHONY : dgemm_mixed_callback_tester
+all :: dgemm_mixed_callback_tester
+endif
+
 dgemv_tester :
 	$(CXX) $(COMMON_OPT) -Wall -Wextra -Wshadow -std=c++11 dgemv_thread_safety.cpp ../$(LIBNAME) $(EXTRALIB) $(FEXTRALIB) -o dgemv_tester
 	./dgemv_tester
@@ -15,5 +35,14 @@ dgemm_mixed_tester : dgemm_tester
 	$(CXX) $(COMMON_OPT) -Wall -Wextra -Wshadow -std=c++11 dgemm_thread_safety_mixed.cpp ../$(LIBNAME) $(EXTRALIB) $(FEXTRALIB) -o dgemm_mixed_tester
 	./dgemm_mixed_tester

+dgemm_shutdown_tester : dgemm_mixed_tester
+	$(CXX) $(COMMON_OPT) -Wall -Wextra -Wshadow -std=c++11 dgemm_thread_safety_shutdown.cpp ../$(LIBNAME) $(EXTRALIB) $(FEXTRALIB) -o dgemm_shutdown_tester
+	./dgemm_shutdown_tester
+
+ifeq ($(CALLBACK_THREAD_TEST),1)
+dgemm_mixed_callback_tester : dgemm_mixed_tester
+	./dgemm_mixed_tester --callback
+endif
+
 clean ::
-	rm -f dgemv_tester dgemm_tester dgemm_mixed_tester
+	rm -f dgemv_tester dgemm_tester dgemm_mixed_tester dgemm_shutdown_tester
@@ -14,6 +14,21 @@
 #endif
 #include "cpp_thread_safety_common.h"

+std::atomic<uint32_t> callbackInvocations(0);
+
+void thread_callback(int sync, openblas_dojob_callback doJob, int numJobs,
+                     size_t jobDataElementSize, void* jobData, int doJobData){
+	(void)sync;
+	callbackInvocations.fetch_add(1, std::memory_order_relaxed);
+	std::vector<std::thread> workers;
+	workers.reserve(numJobs);
+	char* jobs = static_cast<char*>(jobData);
+	for(int i=0; i<numJobs; i++)
+		workers.emplace_back(doJob, i, jobs + i * jobDataElementSize, doJobData);
+	for(auto& worker : workers)
+		worker.join();
+}
+
 void compute_dgemm_pair(std::vector<double>& transA, std::vector<double>& noTransA, std::vector<double>& B, double* firstOutput, double* secondOutput, const blasint randomMatSize, const bool sameVariant){
 	cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, randomMatSize, 2, 2, 1.0, &transA[0], randomMatSize, &B[0], 2, 0.0, firstOutput, 2);
 	if (sameVariant)
@@ -48,26 +63,31 @@ int main(int argc, char* argv[]){
 	uint32_t numTestRounds = 200;
 	uint32_t maxHwThreads = GetMaxHwThreads();
 	bool sameVariant = false;
+	bool useCallback = false;

 	if (maxHwThreads < numConcurrentThreads)
 		numConcurrentThreads = maxHwThreads;

-	if (argc != 1 && argc != 4 && argc != 5){
-		std::cout<<"ERROR: expected zero arguments, or: <M> <threads> <rounds> [sameVariant]"<<std::endl;
+	std::vector<std::string> positionalArgs;
+	for (int i = 1; i < argc; i++){
+		std::cout<<argv[i]<<std::endl;
+		if (std::string(argv[i]) == "--callback")
+			useCallback = true;
+		else
+			positionalArgs.push_back(argv[i]);
+	}
+
+	if (!positionalArgs.empty() && positionalArgs.size() != 3 && positionalArgs.size() != 4){
+		std::cout<<"ERROR: expected: [<M> <threads> <rounds> [sameVariant]] [--callback]"<<std::endl;
 		return 1;
 	}

-	if(argc == 4 || argc == 5){
-		std::vector<std::string> cliArgs;
-		for (int i = 1; i < argc; i++){
-			cliArgs.push_back(argv[i]);
-			std::cout<<argv[i]<<std::endl;
-		}
-		randomMatSize = std::stoul(cliArgs[0]);
-		numConcurrentThreads = std::stoul(cliArgs[1]);
-		numTestRounds = std::stoul(cliArgs[2]);
-		if (argc == 5)
-			sameVariant = std::stoul(cliArgs[3]) != 0;
+	if(!positionalArgs.empty()){
+		randomMatSize = std::stoul(positionalArgs[0]);
+		numConcurrentThreads = std::stoul(positionalArgs[1]);
+		numTestRounds = std::stoul(positionalArgs[2]);
+		if (positionalArgs.size() == 4)
+			sameVariant = std::stoul(positionalArgs[3]) != 0;
 	}

 	FailIfThreadsAreZero(numConcurrentThreads);
@@ -92,6 +112,8 @@ int main(int argc, char* argv[]){
 	std::cout<<"Number of testing rounds : "<<numTestRounds<<'\n';
 	std::cout<<"Second DGEMM uses "<<(sameVariant ? "the same transpose variant" : "a different transpose variant")<<'\n';
 	std::cout<<"OpenBLAS internal threads : "<<openblas_get_num_threads()<<'\n';
+	if (useCallback)
+		std::cout<<"Thread execution backend : caller callback\n";
 	std::cout<<"This test will need "<<(static_cast<uint64_t>(matrixElements) * 2 * 8 + static_cast<uint64_t>(outputElements) * (2 + 2 * numConcurrentThreads) * 8)/static_cast<double>(1024*1024)<<" MiB of RAM\n"<<std::endl;

 	std::cout<<"Filling matrices with deterministic values..."<<std::flush;
@@ -110,6 +132,9 @@ int main(int argc, char* argv[]){
 	compute_dgemm_pair(transA, noTransA, B, &referenceFirst[0], &referenceSecond[0], randomMatSize, sameVariant);
 	std::cout<<"done\n";

+	if (useCallback)
+		openblas_set_threads_callback_function(thread_callback);
+
 	std::cout<<"Testing mixed CBLAS DGEMM thread safety\n";
 	std::cout<<"Launching "<<numConcurrentThreads<<" worker threads..."<<std::flush;
 	for(uint32_t i=0; i<numConcurrentThreads; i++){
@@ -128,6 +153,15 @@ int main(int argc, char* argv[]){
 	}
 	std::cout<<"done\n";

+	if (useCallback) {
+		const uint32_t invocations = callbackInvocations.load();
+		std::cout<<"Thread callback invocations: "<<invocations<<std::endl;
+		if (invocations == 0) {
+			std::cout<<"Thread callback was not invoked!"<<std::endl;
+			return 1;
+		}
+	}
+
 	std::cout<<"Mixed DGEMM mismatches: "<<mismatches<<std::endl;
 	if (mismatches != 0) {
 		std::cout<<"Mixed CBLAS DGEMM thread safety test FAILED!"<<std::endl;
@@ -0,0 +1,165 @@
+/* Stress test for library shutdown racing with in-flight BLAS calls
+ * (https://github.com/OpenMathLib/OpenBLAS/issues/5954).
+ *
+ * Windows only. On POSIX, exit() runs the library destructor while worker
+ * threads are still computing into OpenBLAS-owned buffers, which no amount of
+ * locking inside blas_shutdown can make safe, so there is nothing to assert
+ * there; CMakeLists.txt only registers this test on WIN32.
+ *
+ * The parent re-executes itself as short-lived children and checks that each
+ * one terminates cleanly, turning shutdown-path crashes and deadlocks into
+ * ordinary test failures. Each child (--child-storm N) starts N callers that
+ * allocate their matrices and park on a gate, releases them so they all enter
+ * their first dgemm at once, and exits a millisecond later while that
+ * allocation storm is still in flight.
+ *
+ * N must exceed NUM_BUFFERS = MAX(50, NUM_THREADS * 2 * NUM_PARALLEL) for the
+ * build under test; below that every slot is already mapped and the race is
+ * unreachable.
+ */
+#include <atomic>
+#include <chrono>
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <iostream>
+#include <string>
+#include <thread>
+#include <vector>
+
+#ifdef OPENBLAS_USE_GENERATED_CBLAS_H
+#include "generated/cblas.h"
+#else
+#include "../cblas.h"
+#endif
+
+#include <windows.h>
+
+namespace {
+
+const blasint stormM = 200, stormK = 120, stormN = 90; /* the gh-5954 shape */
+const blasint poolDim = 320; /* above the multithreading threshold, so the pool spins up */
+const uint32_t defaultStormCallers = 128;
+const uint32_t stormDelayMs = 3; /* gate to sweep; at 0 the sweep beats the allocations */
+const int stormBlasThreads = 4;
+const int stormTimeoutSec = 15;
+const int numStormChildren = 40;
+
+std::atomic<uint32_t> parked(0); /* callers built and waiting on the gate */
+std::atomic<bool> gate(false);
+
+void fillOperands(std::vector<double>& A, std::vector<double>& B) {
+	for (size_t i = 0; i < A.size(); i++) A[i] = (i % 1000) / 1000.0;
+	for (size_t i = 0; i < B.size(); i++) B[i] = (i % 997) / 997.0;
+}
+
+void dgemmOnce(blasint m, blasint k, blasint n) {
+	std::vector<double> A(m * k), B(k * n), C(m * n);
+	fillOperands(A, B);
+	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k,
+	            1.0, A.data(), m, B.data(), k, 0.1, C.data(), m);
+}
+
+/* Allocate before parking, so that when the gate opens nothing stands between
+   the thread and its first dgemm. */
+void gatedWorker(blasint m, blasint k, blasint n) {
+	std::vector<double> A(m * k), B(k * n), C(m * n);
+	fillOperands(A, B);
+	parked.fetch_add(1, std::memory_order_release);
+	while (!gate.load(std::memory_order_acquire)) std::this_thread::yield();
+	for (;;)
+		cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k,
+		            1.0, A.data(), m, B.data(), k, 0.1, C.data(), m);
+}
+
+int ChildStorm(uint32_t nCallers) {
+	SetErrorMode(SEM_FAILCRITICALERRORS | SEM_NOGPFAULTERRORBOX);
+	openblas_set_num_threads(stormBlasThreads);
+
+	/* Build the OpenBLAS worker pool first, so the storm is buffer allocation
+	   and not pool startup. */
+	dgemmOnce(poolDim, poolDim, poolDim);
+
+	for (uint32_t i = 0; i < nCallers; i++)
+		std::thread(gatedWorker, stormM, stormK, stormN).detach();
+	for (int ms = 0; parked.load(std::memory_order_acquire) < nCallers && ms < 10000; ms++)
+		std::this_thread::sleep_for(std::chrono::milliseconds(1));
+
+	gate.store(true, std::memory_order_release);
+	std::this_thread::sleep_for(std::chrono::milliseconds(stormDelayMs));
+	std::exit(0);
+}
+
+/* Returns 0 if the child exited cleanly, nonzero otherwise; fills outcome. */
+int RunChild(const std::string& args, int timeoutSec, std::string& outcome) {
+	char exe[MAX_PATH];
+	if (GetModuleFileNameA(NULL, exe, MAX_PATH) == 0) {
+		outcome = "GetModuleFileName failed";
+		return 1;
+	}
+	std::string cmd = "\"" + std::string(exe) + "\" " + args;
+
+	STARTUPINFOA si;
+	PROCESS_INFORMATION pi;
+	ZeroMemory(&si, sizeof(si));
+	si.cb = sizeof(si);
+	ZeroMemory(&pi, sizeof(pi));
+	if (!CreateProcessA(NULL, &cmd[0], NULL, NULL, FALSE, 0, NULL, NULL, &si, &pi)) {
+		outcome = "CreateProcess failed";
+		return 1;
+	}
+
+	int ret = 1;
+	char buf[64];
+	if (WaitForSingleObject(pi.hProcess, timeoutSec * 1000) != WAIT_OBJECT_0) {
+		TerminateProcess(pi.hProcess, 1);
+		WaitForSingleObject(pi.hProcess, 5000);
+		snprintf(buf, sizeof(buf), "HANG (killed after %ds)", timeoutSec);
+	} else {
+		DWORD code = 1;
+		GetExitCodeProcess(pi.hProcess, &code);
+		if (code == 0) {
+			snprintf(buf, sizeof(buf), "clean exit");
+			ret = 0;
+		} else {
+			snprintf(buf, sizeof(buf), "CRASH (exit code 0x%08lX)", (unsigned long)code);
+		}
+	}
+	outcome = buf;
+
+	CloseHandle(pi.hThread);
+	CloseHandle(pi.hProcess);
+	return ret;
+}
+
+} // namespace
+
+int main(int argc, char* argv[]) {
+	if (argc >= 3 && std::strcmp(argv[1], "--child-storm") == 0)
+		return ChildStorm(uint32_t(std::atoi(argv[2])));
+	SetErrorMode(SEM_FAILCRITICALERRORS | SEM_NOGPFAULTERRORBOX);
+
+	uint32_t callers = defaultStormCallers;
+	if (argc >= 2) {
+		int n = std::atoi(argv[1]);
+		if (n > 0) callers = uint32_t(n);
+	}
+
+	int failures = 0;
+	std::cout << "Testing process exit during an allocation storm (" << callers << " callers)"
+	          << std::endl;
+	for (int i = 0; i < numStormChildren; i++) {
+		std::string outcome;
+		failures += RunChild("--child-storm " + std::to_string(callers), stormTimeoutSec, outcome);
+		std::cout << "  storm child " << i << ": " << outcome << std::endl;
+	}
+
+	if (failures) {
+		std::cout << "CBLAS DGEMM shutdown safety test FAILED! (" << failures
+		          << " child processes)" << std::endl;
+		return 1;
+	}
+	std::cout << "CBLAS DGEMM shutdown safety test PASSED!" << std::endl;
+	return 0;
+}
@@ -349,7 +349,7 @@ int CNAME(BLASLONG n, BLASLONG k, FLOAT *alpha, FLOAT *a, BLASLONG lda, FLOAT *x
 #else
 	    ONE, ZERO,
 #endif
-	    (FLOAT*)(queue[i].sb), 1, buffer, 1, NULL, 0);
+	    (FLOAT*)(queue[i].worker_sb), 1, buffer, 1, NULL, 0);
  }

  AXPYU_K(n, 0, 0,
@@ -79,10 +79,32 @@ static int inner_small_matrix_thread(blas_arg_t *args, BLASLONG *range_m, BLASLO
 int CNAME(blas_arg_t * args_array, BLASLONG nums){
  XFLOAT *buffer;
  XFLOAT *sa, *sb;
+  blasint info;
  int nthreads=1;
  int (*routine)(blas_arg_t *, void *, void *, XFLOAT *, XFLOAT *, BLASLONG);
  int i=0, /*j,*/ current_nums;

+#ifndef COMPLEX
+#ifdef XDOUBLE
+#define ERROR_NAME "QGEMM_BATCH "
+#elif defined(DOUBLE)
+#define ERROR_NAME "DGEMM_BATCH "
+#elif defined(BFLOAT16)
+#define ERROR_NAME "SBGEMM_BATCH "
+#else
+#define ERROR_NAME "SGEMM_BATCH "
+#endif
+#else
+#ifdef XDOUBLE
+#define ERROR_NAME "XGEMM_BATCH "
+#elif defined(DOUBLE)
+#define ERROR_NAME "ZGEMM_BATCH "
+#else
+#define ERROR_NAME "CGEMM_BATCH "
+#endif
+#endif
+
+
 #ifdef SMP
  blas_queue_t * queue=NULL;
 #endif
@@ -90,6 +112,12 @@ int CNAME(blas_arg_t * args_array, BLASLONG nums){
  if(nums <=0 ) return 0;

  buffer = (XFLOAT *)blas_memory_alloc(0);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return(1); 
+  }
+
  sa = (XFLOAT *)((BLASLONG)buffer +GEMM_OFFSET_A);
  sb = (XFLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
  
@@ -119,6 +147,9 @@ int CNAME(blas_arg_t * args_array, BLASLONG nums){
    queue=(blas_queue_t *)malloc((nums+1) * sizeof(blas_queue_t));
    if(queue == NULL){
      openblas_warning(0, "memory alloc failed!\n");
+      if (buffer) blas_memory_free(buffer);
+      info = -999;
+      BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
      return(1);
    }
    for(i=0; i<nums; i++){
@@ -37,6 +37,11 @@
 /* or implied, of The University of Texas at Austin.                 */
 /*********************************************************************/

+#if defined(__loongarch__) && defined(LA464) && defined(DOUBLE)
+#undef GEMM_P
+#define GEMM_P 96
+#endif
+
 #ifndef CACHE_LINE_SIZE
 #define CACHE_LINE_SIZE 8
 #endif
@@ -335,7 +335,7 @@ static void legacy_exec(void *func, int mode, blas_arg_t *args, void *sb){
      }
 }

-#if defined(OS_LINUX) && !defined(NO_AFFINITY) 
+#if defined(OS_LINUX) && !defined(NO_AFFINITY)
 int gotoblas_set_affinity(int);
 int gotoblas_set_affinity2(int);
 int get_node(void);
@@ -345,7 +345,7 @@ static int increased_threads = 0;

 #ifdef OS_LINUX
 extern int openblas_get_num_threads(void);  
-
+OPENBLAS_EXPORT
 int openblas_setaffinity(int thread_idx, size_t cpusetsize, cpu_set_t* cpu_set) {
  const int active_threads = openblas_get_num_threads();

@@ -360,6 +360,7 @@ int openblas_setaffinity(int thread_idx, size_t cpusetsize, cpu_set_t* cpu_set)

  return pthread_setaffinity_np(thread, cpusetsize, cpu_set);
 }
+OPENBLAS_EXPORT
 int openblas_getaffinity(int thread_idx, size_t cpusetsize, cpu_set_t* cpu_set) {
  const int active_threads = openblas_get_num_threads();

@@ -836,6 +837,7 @@ if (openblas_threads_callback_) {
 #endif

  routine = (int (*)(blas_arg_t *, void *, void *, double *, double *, BLASLONG))queue -> routine;
+  queue->worker_sb = queue->sb;

  if (queue -> mode & BLAS_LEGACY) {
    legacy_exec(routine, queue -> mode, queue -> args, queue -> sb);
@@ -1139,9 +1141,10 @@ if (!(queue -> mode & BLAS_COMPLEX)){
        /* Other types in future */
    }
 }
-queue->sb=sb;
    }

+    queue->worker_sb = sb;
+
 #ifdef MONITOR
 main_status[cpu] = MAIN_RUNNING2;
 #endif
@@ -118,17 +118,19 @@ void goto_set_num_threads(int num_threads) {
 #endif

 }
+OPENBLAS_EXPORT
 void openblas_set_num_threads(int num_threads) {
 	blas_is_num_threads_set_explicitly = 1;
 	goto_set_num_threads(num_threads);
 }

 #ifdef OS_LINUX
-
+OPENBLAS_EXPORT
 int openblas_setaffinity(int thread_idx, size_t cpusetsize, cpu_set_t* cpu_set) {
  fprintf(stderr,"OpenBLAS: use OpenMP environment variables for setting cpu affinity\n");
  return -1;
 }
+OPENBLAS_EXPORT
 int openblas_getaffinity(int thread_idx, size_t cpusetsize, cpu_set_t* cpu_set) {
  fprintf(stderr,"OpenBLAS: use OpenMP environment variables for querying cpu affinity\n");
  return -1;
@@ -327,7 +329,6 @@ static void exec_threads(int thread_num, blas_queue_t *queue, int buf_index){

    if (sa == NULL) {
      sa = (void *)((BLASLONG)buffer + GEMM_OFFSET_A);
-      queue->sa=sa;
    }

    if (sb == NULL) {
@@ -376,10 +377,11 @@ fprintf(stderr,"UNHANDLED COMPLEX\n");
          /* Other types in future */
 	  }
      }
-      queue->sb=sb;
    }
  }

+  queue->worker_sb = sb;
+
  if (queue -> mode & BLAS_LEGACY) {
    legacy_exec(queue -> routine, queue -> mode, queue -> args, sb);
  } else
@@ -311,9 +311,10 @@ static DWORD WINAPI blas_thread_server(void *arg){
            /* Other types in future */
 	    }
 	}
-	queue->sb=sb;
      }

+      queue->worker_sb = sb;
+
 #ifdef MONITOR
      main_status[cpu] = MAIN_RUNNING2;
 #endif
@@ -483,6 +484,7 @@ int exec_blas(BLASLONG num, blas_queue_t *queue){
  if ((num > 1) && queue -> next) exec_blas_async(1, queue -> next);

  routine = queue -> routine;
+  queue->worker_sb = queue->sb;

  if (queue -> mode & BLAS_LEGACY) {
    legacy_exec(routine, queue -> mode, queue -> args, queue -> sb);
@@ -171,6 +171,8 @@ extern gotoblas_t  gotoblas_CORTEXA53;
 #else
 extern gotoblas_t  gotoblas_CORTEXA57;
 #endif
+/* Keep A72 out of the default DYNAMIC_CORE set (see #4389); opt in via
+ * DYNAMIC_LIST=CORTEXA72. Default DYNAMIC_ARCH continues to alias A57. */
 #define gotoblas_CORTEXA72 gotoblas_CORTEXA57
 #define gotoblas_CORTEXA73 gotoblas_CORTEXA57
 #define gotoblas_FALKOR gotoblas_CORTEXA57
@@ -2099,7 +2099,7 @@ int openblas_get_num_threads(void) {

 struct release_t {
  void *address;
-  void (*func)(struct release_t *);
+  void (* _Atomic func)(struct release_t *);
  long attr;
 };

@@ -2123,6 +2123,37 @@ static pthread_spinlock_t alloc_lock = 0;
 static BLASULONG  alloc_lock = 0UL;
 #endif

+static void blas_release_register(void *address, void (*func)(struct release_t *), long attr) {
+
+  struct release_t *release;
+  int rpos;
+
+#if (defined(SMP) || defined(USE_LOCKING)) && !defined(USE_OPENMP)
+  LOCK_COMMAND(&alloc_lock);
+#endif
+#if defined(HAVE_C11) && !defined(__cplusplus)
+  rpos = atomic_fetch_add(&release_pos, 1);
+#elif defined(__GNUC__)
+  rpos = __sync_fetch_and_add(&release_pos, 1);
+#elif defined(OS_WINDOWS)
+  rpos = InterlockedIncrement((LONG volatile *)&release_pos) - 1;
+#else
+  rpos = release_pos++;
+#endif
+  if (likely(rpos < NUM_BUFFERS)) {
+    release = &release_info[rpos];
+  } else {
+    release = &new_release_info[rpos - NUM_BUFFERS];
+  }
+  release->address = address;
+  release->attr    = attr;
+  WMB;
+  release->func    = func;
+#if (defined(SMP) || defined(USE_LOCKING)) && !defined(USE_OPENMP)
+  UNLOCK_COMMAND(&alloc_lock);
+#endif
+}
+
 #ifdef ALLOC_MMAP

 static void alloc_mmap_free(struct release_t *release){
@@ -2154,20 +2185,7 @@ static void *alloc_mmap(void *address){
  }

  if (map_address != (void *)-1) {
-#if (defined(SMP) || defined(USE_LOCKING)) && !defined(USE_OPENMP)
-    LOCK_COMMAND(&alloc_lock);
-#endif
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].func    = alloc_mmap_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_mmap_free;
-    }
-#if (defined(SMP) || defined(USE_LOCKING)) && !defined(USE_OPENMP)
-    UNLOCK_COMMAND(&alloc_lock);
-#endif
+    blas_release_register(map_address, alloc_mmap_free, 0);
  } else {
 #ifdef DEBUG
        int errsv=errno;
@@ -2323,20 +2341,7 @@ static void *alloc_mmap(void *address){
 #endif

  if (map_address != (void *)-1) {
-#if (defined(SMP) || defined(USE_LOCKING)) && !defined(USE_OPENMP)
-    LOCK_COMMAND(&alloc_lock);
-#endif
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].func    = alloc_mmap_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_mmap_free;
-    }
-#if (defined(SMP) || defined(USE_LOCKING)) && !defined(USE_OPENMP)
-    UNLOCK_COMMAND(&alloc_lock);
-#endif
+    blas_release_register(map_address, alloc_mmap_free, 0);
  }

  return map_address;
@@ -2364,14 +2369,7 @@ static void *alloc_malloc(void *address){
  if (map_address == (void *)NULL) map_address = (void *)-1;

  if (map_address != (void *)-1) {
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].func    = alloc_malloc_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_malloc_free;
-    }
+    blas_release_register(map_address, alloc_malloc_free, 0);
  }

  return map_address;
@@ -2403,14 +2401,7 @@ static void *alloc_qalloc(void *address){
  if (map_address == (void *)NULL) map_address = (void *)-1;

  if (map_address != (void *)-1) {
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].func    = alloc_qalloc_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_qalloc_free;
-    }
+    blas_release_register(map_address, alloc_qalloc_free, 0);
  }

  return (void *)(((BLASULONG)map_address + FIXED_PAGESIZE - 1) & ~(FIXED_PAGESIZE - 1));
@@ -2437,14 +2428,7 @@ static void *alloc_windows(void *address){
  if (map_address == (void *)NULL) map_address = (void *)-1;

  if (map_address != (void *)-1) {
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].func    = alloc_windows_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_windows_free;
-    }
+    blas_release_register(map_address, alloc_windows_free, 0);
  }

  return map_address;
@@ -2486,16 +2470,7 @@ static void *alloc_devicedirver(void *address){
                     fd, 0);

  if (map_address != (void *)-1) {
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].attr    = fd;
-    release_info[rpos].func    = alloc_devicedirver_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].attr    = fd;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_devicedirver_free;
-    }
+    blas_release_register(map_address, alloc_devicedirver_free, fd);
  }

  return map_address;
@@ -2530,16 +2505,7 @@ static void *alloc_shm(void *address){

    shmctl(shmid, IPC_RMID, 0);

-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].attr    = shmid;
-    release_info[rpos].func    = alloc_shm_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].attr    = shmid;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_shm_free;
-    }
+    blas_release_register(map_address, alloc_shm_free, shmid);
  }

  return map_address;
@@ -2647,14 +2613,7 @@ fprintf(stderr,"alloc_hugetlb got called\n");
 #endif

  if (map_address != (void *)-1){
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].func    = alloc_hugetlb_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_hugetlb_free;
-    }
+    blas_release_register(map_address, alloc_hugetlb_free, 0);
  }

  return map_address;
@@ -2699,16 +2658,7 @@ static void *alloc_hugetlbfile(void *address){
                     fd, 0);

  if (map_address != (void *)-1) {
-    int rpos = release_pos++;
-    if (likely(rpos < NUM_BUFFERS)) {
-    release_info[rpos].address = map_address;
-    release_info[rpos].attr    = fd;
-    release_info[rpos].func    = alloc_hugetlbfile_free;
-    } else {
-    new_release_info[rpos-NUM_BUFFERS].address = map_address;
-    new_release_info[rpos-NUM_BUFFERS].attr    = fd;
-    new_release_info[rpos-NUM_BUFFERS].func    = alloc_hugetlbfile_free;
-    }
+    blas_release_register(map_address, alloc_hugetlbfile_free, fd);
  }

  return map_address;
@@ -3056,7 +3006,9 @@ void *blas_memory_alloc(int procpos){
 #endif
  memory_overflowed=1;
  MB;
-  new_release_info = (struct release_t*) malloc(NEW_BUFFERS * sizeof(struct release_t));
+  /* zeroed so blas_shutdown sees NULL func in slots that were reserved but
+     never published */
+  new_release_info = (struct release_t*) calloc(NEW_BUFFERS, sizeof(struct release_t));
  newmemory = (struct newmemstruct*) malloc(NEW_BUFFERS * sizeof(struct newmemstruct));
  for (i = 0; i < NEW_BUFFERS; i++) {
  newmemory[i].addr   = (void *)0;
@@ -3240,9 +3192,41 @@ void blas_memory_free_nolock(void * map_address) {
  free(map_address);
 }

+#if defined(OS_WINDOWS) && !defined(OS_CYGWIN_NT)
+/* During process termination Windows has already killed every other thread,
+   possibly while one held alloc_lock or a blas server lock, so any cleanup
+   here can only deadlock or crash; the OS reclaims the memory anyway.
+
+   Nothing in the SDK reports this from a destructor. DllMain's lpReserved
+   distinguishes the two cases, but outside MSVC gotoblas_quit runs from the
+   CRT's fini array via __attribute__((destructor)) and never sees
+   it. RtlDllShutdownInProgress is documented under Win32 Dev Notes but
+   deliberately absent from the SDK headers, so callers declare it themselves
+   If it cannot be resolved we fall back to the previous behaviour. */
+typedef BOOLEAN (WINAPI *rtl_dll_shutdown_in_progress_t)(VOID);
+static rtl_dll_shutdown_in_progress_t rtl_dll_shutdown_in_progress = NULL;
+
+/* Resolved at init, not on the way out: GetModuleHandle takes LdrpSnapsLock,
+   which ExitProcess does not release before it kills the other threads. */
+static void blas_shutdown_check_init(void) {
+  HMODULE ntdll = GetModuleHandleA("ntdll.dll");
+  if (!ntdll) return;
+  rtl_dll_shutdown_in_progress = (rtl_dll_shutdown_in_progress_t)(void *)
+      GetProcAddress(ntdll, "RtlDllShutdownInProgress");
+}
+
+static int blas_process_is_terminating(void) {
+  return rtl_dll_shutdown_in_progress && rtl_dll_shutdown_in_progress();
+}
+#endif
+
 void blas_shutdown(void){

-  int pos;
+  int pos, release_count;
+
+#if defined(OS_WINDOWS) && !defined(OS_CYGWIN_NT)
+  if (blas_process_is_terminating()) return;
+#endif

 #ifdef SMP
  BLASFUNC(blas_thread_shutdown)();
@@ -3250,12 +3234,18 @@ void blas_shutdown(void){

  LOCK_COMMAND(&alloc_lock);

-  for (pos = 0; pos < release_pos; pos ++) {
-    if (likely(pos < NUM_BUFFERS))
-    release_info[pos].func(&release_info[pos]);
-    else
-    new_release_info[pos-NUM_BUFFERS].func(&new_release_info[pos-NUM_BUFFERS]);
+  release_count = release_pos;
+  for (pos = 0; pos < release_count; pos ++) {
+    struct release_t *release = likely(pos < NUM_BUFFERS) ?
+        &release_info[pos] : &new_release_info[pos-NUM_BUFFERS];
+    void (*func)(struct release_t *) = release->func;
+    RMB;
+    if (func == NULL) continue; /* reserved but never published: owner died mid-allocation */
+    func(release);
+    release->func = NULL;
+    release->address = NULL;
  }
+  release_pos = 0;

 #ifdef SEEK_ADDRESS
  base_address      = 0UL;
@@ -3282,6 +3272,8 @@ void blas_shutdown(void){
    }
    free((void*)newmemory);
    newmemory = NULL;
+    free(new_release_info);
+    new_release_info = NULL;
    memory_overflowed = 0;
  }

@@ -3406,6 +3398,10 @@ void CONSTRUCTOR gotoblas_init(void) {

  if (gotoblas_initialized) return;

+#if defined(OS_WINDOWS) && !defined(OS_CYGWIN_NT)
+  blas_shutdown_check_init();
+#endif
+
 #ifdef SMP
  openblas_fork_handler();
 #endif
@@ -226,6 +226,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0) y -= (leny-1)*incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (m * n < 250000 || kl+ku < 15 )
@@ -268,7 +268,6 @@ void NAME(char *TRANSA, char *TRANSB,

  int transa, transb, nrowa, nrowb;
  blasint info;
-  int order = -1;

  char transA, transB;
  IFLOAT *buffer;
@@ -317,7 +316,7 @@ void NAME(char *TRANSA, char *TRANSB,

  args.alpha = (void *)alpha;
  args.beta  = (void *)beta;
-
+  
  transA = *TRANSA;
  transB = *TRANSB;

@@ -346,7 +345,6 @@ void NAME(char *TRANSA, char *TRANSB,
  if (transB == 'R') transb = 2;
  if (transB == 'C') transb = 3;
 #endif
-
  nrowa = args.m;
  if (transa & 1) nrowa = args.k;
  nrowb = args.k;
@@ -562,22 +560,80 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
 || strcmp(gotoblas_corename(), "vortexm4") == 0
 #endif
 )
-// if (support_sme1())
 #endif
+
  if (order == CblasRowMajor && k==lda && n==ldb && n==ldc && beta == 0 && alpha == 1.0 && TransA == CblasNoTrans && TransB == CblasNoTrans && SGEMM_DIRECT_PERFORMANT(m,n,k)) {
        SGEMM_DIRECT(m, n, k, a, lda, b, ldb, c, ldc);
        return;
  }
-else
- if (order == CblasRowMajor && k==lda && n==ldb && n==ldc && TransA == CblasNoTrans && TransB == CblasNoTrans && SGEMM_DIRECT_PERFORMANT(m,n,k)) {
+else if (order == CblasRowMajor && k==lda && n==ldb && n==ldc && TransA == CblasNoTrans && TransB == CblasNoTrans && SGEMM_DIRECT_PERFORMANT(m,n,k)) {
        SGEMM_DIRECT_ALPHA_BETA(m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
        return;
  }

-#endif
-#endif
+#endif //defined arm64
+#endif //defined complex

+
+#endif //defined CBLAS
+
+#if !defined(BFLOAT16)  && !defined(HFLOAT16)
+#if defined(ARCH_ARM64) && (defined(USE_SGEMM_KERNEL_DIRECT)||defined(DYNAMIC_ARCH))
+#if defined(DYNAMIC_ARCH)
+if (strcmp(gotoblas_corename(), "armv9sme") == 0
+#if defined(__clang__)
+ || strcmp(gotoblas_corename(), "vortexm4") == 0
 #endif
+)
+#endif //defined dynarch
+{
+char* TA,*TB;
+      if (transa & 1) 
+        TA = "T";
+	else 
+	TA= "N";
+      if (transb & 1) 
+        TB = "T";
+	else 
+	TB= "N";
+#ifndef COMPLEX
+      if (transa == 3) 
+	TA= "T";
+      if (transb == 3) 
+	TB= "T";
+FLOAT* al=(FLOAT*)args.alpha;
+FLOAT* be=(FLOAT*)args.beta;
+#ifndef DOUBLE
+        SME_SGEMM_KERNEL(TA,TB, args.m, args.n, args.k, al, args.a, args.lda, args.b, args.ldb, be, args.c, args.ldc);
+#else
+        SME_DGEMM_KERNEL(TA,TB, args.m, args.n, args.k, al, args.a, args.lda, args.b, args.ldb, be, args.c, args.ldc);
+#endif
+#else
+      if (transa == 2) 
+	TA= "R";
+      if (transb == 2) 
+	TB= "R";
+      if (transa == 3) 
+	TA= "C";
+      if (transb == 3) 
+	TB= "C";
+FLOAT* al=(FLOAT*)args.alpha;
+FLOAT* be=(FLOAT*)args.beta;
+#ifndef DOUBLE
+float _Complex c_al={al[0],al[1]};
+float _Complex c_be={be[0],be[1]};
+        SME_CGEMM_KERNEL(TA,TB, args.m, args.n, args.k, al[0],al[1], args.a, args.lda, args.b, args.ldb, be[0],be[1], args.c, args.ldc);
+#else
+double _Complex c_al={al[0],al[1]};
+double _Complex c_be={be[0],be[1]};
+        SME_ZGEMM_KERNEL(TA,TB, args.m, args.n, args.k, al[0],al[1], args.a, args.lda, args.b, args.ldb, be[0],be[1], args.c, args.ldc);
+#endif
+#endif
+	return;
+}
+#endif //defined arm64
+
+#endif //defined b/hfloat16 

 #if defined(__linux__) && defined(__x86_64__) && defined(BFLOAT16)
 #if defined(DYNAMIC_ARCH)
@@ -679,6 +735,7 @@ else
 #if USE_SMALL_MATRIX_OPT
 #if !defined(COMPLEX)
  if(GEMM_SMALL_MATRIX_PERMIT(transa, transb, args.m, args.n, args.k, *(FLOAT *)(args.alpha), *(FLOAT *)(args.beta))){
+
 	  if(*(FLOAT *)(args.beta) == 0.0){
 		(GEMM_SMALL_KERNEL_B0((transb << 2) | transa))(args.m, args.n, args.k, args.a, args.lda, *(FLOAT *)(args.alpha), args.b, args.ldb, args.c, args.ldc);
 	  }else{
@@ -699,7 +756,6 @@ else
 #endif

  buffer = (XFLOAT *)blas_memory_alloc(0);
-
 //For LOONGARCH64, applying an offset to the buffer is essential
 //for minimizing cache conflicts and optimizing performance.
 #if defined(ARCH_LOONGARCH64) && !defined(NO_AFFINITY)
@@ -353,18 +353,18 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE *  transa_array, enum CB
 #if !defined(COMPLEX)
      if(*(FLOAT *)(group_beta) == 0.0){
 	group_mode=mode | BLAS_SMALL_B0_OPT;
-	group_small_matrix_opt_routine=(void *)(gemm_small_kernel_b0[(group_transb<<2)|group_transa]);
+	group_small_matrix_opt_routine=SMALL_KERNEL_ADDR(gemm_small_kernel_b0, ((group_transb<<2)|group_transa));
      }else{
 	group_mode=mode | BLAS_SMALL_OPT;
-	group_small_matrix_opt_routine=(void *)(gemm_small_kernel[(group_transb<<2)|group_transa]);
+	group_small_matrix_opt_routine=SMALL_KERNEL_ADDR(gemm_small_kernel, ((group_transb<<2)|group_transa));
      }
 #else
      if(((FLOAT *)(group_beta))[0] == 0.0 && ((FLOAT *)(group_beta))[1] == 0.0){
 	group_mode=mode | BLAS_SMALL_B0_OPT;
-	group_small_matrix_opt_routine=(void *)(zgemm_small_kernel_b0[(group_transb<<2)|group_transa]);
+	group_small_matrix_opt_routine=SMALL_KERNEL_ADDR(zgemm_small_kernel_b0, ((group_transb<<2)|group_transa));
      }else{
 	group_mode=mode | BLAS_SMALL_OPT;
-	group_small_matrix_opt_routine=(void *)(zgemm_small_kernel[(group_transb<<2)|group_transa]);
+	group_small_matrix_opt_routine=SMALL_KERNEL_ADDR(zgemm_small_kernel, ((group_transb<<2)|group_transa));
      }

 #endif
@@ -229,13 +229,16 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE transa, enum CBLAS_TRANS
    group_transb = -1;
    info   = 0;
   
-#if defined(CBLAS) 
+#if defined(CBLAS)
    if (order == CblasColMajor) {
      group_m = m;
      group_n = n;
      group_k = k;

-      
+      group_lda = lda;
+      group_ldb = ldb;
+      group_ldc = ldc;
+
      if (transa == CblasNoTrans)     group_transa = 0;
      if (transa == CblasTrans)       group_transa = 1;
 #ifndef COMPLEX
@@ -74,7 +74,7 @@ int NAME(blasint *k, blasint *n, blasint *n1, FLOAT *d,

 #ifdef SMP
  int nthreads = 1;
-  if (n >= 64) nthreads = num_cpu_avail(4);
+  if (nval >= 64) nthreads = num_cpu_avail(4);

  if (nthreads == 1) {
 #endif
@@ -198,6 +198,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0 ) y -= (n - 1) * incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMPTEST
  nthreads = num_cpu_avail(2);
@@ -182,6 +182,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0 ) y -= (n - 1) * incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMPTEST
  nthreads = num_cpu_avail(2);
@@ -192,6 +192,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -196,6 +196,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0 ) y -= (n - 1) * incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -405,6 +405,11 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
  FUNCTION_PROFILE_START();

  buffer = (FLOAT *)blas_memory_alloc(0);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  sa = (FLOAT *)((BLASLONG)buffer + GEMM_OFFSET_A);
  sb = (FLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
@@ -180,6 +180,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, blasint n, FLOAT alpha,
  if (incy < 0 ) y -= (n - 1) * incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (n <200)
@@ -195,6 +195,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, blasint n, FLOAT alpha,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -197,6 +197,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, blasint n, FLOAT alpha,
  if (incy < 0 ) y -= (n - 1) * incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -393,6 +393,11 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
  FUNCTION_PROFILE_START();

  buffer = (FLOAT *)blas_memory_alloc(0);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  sa = (FLOAT *)((BLASLONG)buffer + GEMM_OFFSET_A);
  sb = (FLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
@@ -376,6 +376,11 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
  FUNCTION_PROFILE_START();

  buffer = (FLOAT *)blas_memory_alloc(0);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  sa = (FLOAT *)((BLASLONG)buffer + GEMM_OFFSET_A);
  sb = (FLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
@@ -223,6 +223,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -202,6 +202,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  (tbsv[(trans<<2) | (uplo<<1) | unit])(n, k, a, lda, x, incx, buffer);

@@ -221,6 +221,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -193,6 +193,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  (tpsv[(trans<<2) | (uplo<<1) | unit])(n, a, x, incx, buffer);

@@ -218,6 +218,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -397,6 +397,11 @@ if (strcmp(gotoblas_corename(), "armv9sme") == 0
  FUNCTION_PROFILE_START();

  buffer = (FLOAT *)blas_memory_alloc(0);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  sa = (FLOAT *)((BLASLONG)buffer + GEMM_OFFSET_A);
  sb = (FLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
@@ -197,6 +197,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  (trsv[(trans<<2) | (uplo<<1) | unit])(n, a, lda, x, incx, buffer);

@@ -251,6 +251,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0) y -= (leny - 1) * incy * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (m * n  < 125000 || ku + kl < 15) 
@@ -204,6 +204,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0 ) y -= (n - 1) * incy * COMPSIZE;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMPBUG
  nthreads = num_cpu_avail(2);
@@ -199,6 +199,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, blasint n, void *VALPHA
  if (incy < 0 ) y -= (n - 1) * incy * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (n<MULTI_THREAD_MINIMAL) {
@@ -177,6 +177,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, blasint n, FLOAT alpha,
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (n < 100)
@@ -186,6 +186,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, blasint n, void *VALPHA
  if (incy < 0 ) y -= (n - 1) * incy * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (n < 100)
@@ -194,6 +194,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0 ) y -= (n - 1) * incy * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -175,6 +175,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (n < 100)
@@ -187,6 +187,11 @@ void CNAME(enum CBLAS_ORDER order,
  if (incy < 0 ) y -= (n - 1) * incy * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  if (n < 100)
@@ -44,140 +44,145 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {

  FUNCTION_PROFILE_START();

-  if (db_r == ZERO && db_i == ZERO) {
-    *C        = ONE;
-    *(S  + 0) = ZERO;
-    *(S  + 1) = ZERO;
-    return;
-  }
-
-  long double safmax = 1./safmin;
-#if defined DOUBLE
-  long double rtmax = safmax /DBL_EPSILON;
-#else
-  long double rtmax = safmax /FLT_EPSILON;
-#endif
-		*(S1 + 0) = *(DB + 0);
-		*(S1 + 1) = *(DB + 1) *-1;
-	if (da_r == ZERO && da_i == ZERO) {
-	    *C = ZERO;
-	    if (db_r == ZERO) {
-		    (*DA) = fabsl(db_i);
-		*S = *S1 /(*DA);
-		*(S+1) = *(S1+1) /(*DA);
-		return;
-	    } else if ( db_i == ZERO) {
-		    *DA = fabsl(db_r);
-		*S = *S1 /(*DA);
-		*(S+1) = *(S1+1) /(*DA);
-		return;
-	    } else {
-	        long double g1 = MAX( fabsl(db_r), fabsl(db_i));
-	        rtmax =sqrt(safmax/2.);
-	        if (g1 > rtmin && g1 < rtmax) { // unscaled
-		    d = sqrt(adb);
-		    *S = *S1 /d;
-		    *(S+1) = *(S1+1) /d;
-		    *DA = d ;
-		    *(DA+1) = ZERO;
-		    return;
-    	        } else { // scaled algorithm
-		  long double u = MIN ( safmax, MAX ( safmin, g1));
-		    FLOAT gs_r = db_r/u;
-		    FLOAT gs_i = db_i/u;
-		    d = sqrt ( gs_r*gs_r + gs_i*gs_i);
-		    *S = gs_r / d;
-		    *(S + 1) = (gs_i * -1) / d;
-		    *DA = d * u;
-		    *(DA+1) = ZERO;
-		    return;
-	        }
-	    }
-	} else {
-	       FLOAT f1 = MAX ( fabsl(da_r), fabsl(da_i));
-	       FLOAT g1 = MAX ( fabsl(db_r), fabsl(db_i));
-	       rtmax = sqrt(safmax / 4.);
-	       if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax) { //unscaled
-		    long double h = ada + adb;
-	   	    double adahsq = sqrt(ada * h);
-		    if (ada >= h *safmin) {
-			*C = sqrt(ada/h);
-			*R = *DA / *C;
-			*(R+1) = *(DA+1) / *C;
-			rtmax *= 2.;
-			if ( ada > rtmin && h < rtmax) { // no risk of intermediate overflow
-				*S = *S1 * (*DA / adahsq) - *(S1+1)* (*(DA+1)/adahsq);
-				*(S+1) = *S1 * (*(DA+1) / adahsq) + *(S1+1) * (*DA/adahsq);
-			} else {
-				*S = *S1 * (*R/h) - *(S1+1) * (*(R+1)/h);
-				*(S+1) = *S1 * (*(R+1)/h) + *(S1+1) * (*(R)/h);
-			}
-	    	    } else {
-		        *C = ada / adahsq;
-		        if (*C >= safmin) {
-			    *R = *DA / *C;
-			    *(R+1) = *(DA+1) / *C;
-			} else {
-			    *R = *DA * (h / adahsq);
-			    *(R+1) = *(DA+1) * (h / adahsq);
-			}
-		        *S = *S1 * ada / adahsq;
-		    	*(S+1) = *(S1+1) * ada / adahsq;
-		    }
-		    *DA=*R;
-		    *(DA+1)=*(R+1);
-		    return;
-	        } else {	// scaled
-		   FLOAT fs_r, fs_i, gs_r, gs_i;
-		long double v,w,f2,g2,h;
-		long double u = MIN ( safmax, MAX ( safmin, MAX(f1,g1)));
-		    gs_r = db_r/u;
-		    gs_i = db_i/u;
-		    g2 = sqrt ( gs_r*gs_r + gs_i*gs_i);
-		    if (f1 /u < rtmin) {
-			v = MIN (safmax, MAX (safmin, f1));
-			w = v / u;
-			fs_r = *DA/ v;
-			fs_i = *(DA+1) / v;
-		        f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
-		        h = f2 * w * w + g2;		
-		    } else { // use same scaling for both
-			w = 1.;
-			fs_r = *DA/ u;
-			fs_i = *(DA+1) / u;
-		        f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
-			h = f2 + g2;
-		    }
-		    if ( f2 >= h * safmin) {
-			    *C = sqrt ( f2 / h );
-			    *DA = fs_r / *C;
-			    *(DA+1) = fs_i / *C;
-			    rtmax *= 2;
-			    if ( f2 > rtmin && h < rtmax) {
-				    *S = gs_r * (fs_r /sqrt(f2*h)) - gs_i * (fs_i / sqrt(f2*h));
-				    *(S+1) = gs_r * (fs_i /sqrt(f2*h)) + gs_i * -1. * (fs_r / sqrt(f2*h));
-			    } else {
-				    *S = gs_r * (*DA/h) - gs_i * (*(DA+1) / h);
-				    *(S+1) = gs_r * (*(DA+1) /h) + gs_i * -1. * (*DA / h);
-			    }
-		    } else { // intermediates might overflow
-			d = sqrt ( f2 * h);
-			*C = f2 /d;
-			if (*C >= safmin) {
-				*DA = fs_r / *C;
-				*(DA+1) = fs_i / *C;
-			} else {
-				*DA = fs_r * (h / d);
-				*(DA+1) = fs_i / (h / d);
-			}
-			*S = gs_r * (fs_r /d) - gs_i * (fs_i / d);
-			*(S+1) = gs_r * (fs_i /d) + gs_i * -1. * (fs_r / d);
-	            }
-		    *C *= w;
-		    *DA *= u;
-		    *(DA+1) *= u;
-		    return;
+  do {
+		if (db_r == ZERO && db_i == ZERO) {
+			*C        = ONE;
+			*(S  + 0) = ZERO;
+			*(S  + 1) = ZERO;
+			break;
 		}
-	}
+
+		long double safmax = 1./safmin;
+	#if defined DOUBLE
+		long double rtmax = safmax /DBL_EPSILON;
+	#else
+		long double rtmax = safmax /FLT_EPSILON;
+	#endif
+			*(S1 + 0) = *(DB + 0);
+			*(S1 + 1) = *(DB + 1) *-1;
+		if (da_r == ZERO && da_i == ZERO) {
+				*C = ZERO;
+				if (db_r == ZERO) {
+					(*DA) = fabsl(db_i);
+			*S = *S1 /(*DA);
+			*(S+1) = *(S1+1) /(*DA);
+			break;
+				} else if ( db_i == ZERO) {
+					*DA = fabsl(db_r);
+			*S = *S1 /(*DA);
+			*(S+1) = *(S1+1) /(*DA);
+			break;
+				} else {
+						long double g1 = MAX( fabsl(db_r), fabsl(db_i));
+						rtmax =sqrt(safmax/2.);
+						if (g1 > rtmin && g1 < rtmax) { // unscaled
+					d = sqrt(adb);
+					*S = *S1 /d;
+					*(S+1) = *(S1+1) /d;
+					*DA = d ;
+					*(DA+1) = ZERO;
+					break;
+								} else { // scaled algorithm
+				long double u = MIN ( safmax, MAX ( safmin, g1));
+					FLOAT gs_r = db_r/u;
+					FLOAT gs_i = db_i/u;
+					d = sqrt ( gs_r*gs_r + gs_i*gs_i);
+					*S = gs_r / d;
+					*(S + 1) = (gs_i * -1) / d;
+					*DA = d * u;
+					*(DA+1) = ZERO;
+					break;
+						}
+				}
+		} else {
+					FLOAT f1 = MAX ( fabsl(da_r), fabsl(da_i));
+					FLOAT g1 = MAX ( fabsl(db_r), fabsl(db_i));
+					rtmax = sqrt(safmax / 4.);
+					if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax) { //unscaled
+					long double h = ada + adb;
+						double adahsq = sqrt(ada * h);
+					if (ada >= h *safmin) {
+				*C = sqrt(ada/h);
+				*R = *DA / *C;
+				*(R+1) = *(DA+1) / *C;
+				rtmax *= 2.;
+				if ( ada > rtmin && h < rtmax) { // no risk of intermediate overflow
+					*S = *S1 * (*DA / adahsq) - *(S1+1)* (*(DA+1)/adahsq);
+					*(S+1) = *S1 * (*(DA+1) / adahsq) + *(S1+1) * (*DA/adahsq);
+				} else {
+					*S = *S1 * (*R/h) - *(S1+1) * (*(R+1)/h);
+					*(S+1) = *S1 * (*(R+1)/h) + *(S1+1) * (*(R)/h);
+				}
+							} else {
+							*C = ada / adahsq;
+							if (*C >= safmin) {
+						*R = *DA / *C;
+						*(R+1) = *(DA+1) / *C;
+				} else {
+						*R = *DA * (h / adahsq);
+						*(R+1) = *(DA+1) * (h / adahsq);
+				}
+							*S = *S1 * ada / adahsq;
+						*(S+1) = *(S1+1) * ada / adahsq;
+					}
+					*DA=*R;
+					*(DA+1)=*(R+1);
+					break;
+						} else {	// scaled
+				FLOAT fs_r, fs_i, gs_r, gs_i;
+			long double v,w,f2,g2,h;
+			long double u = MIN ( safmax, MAX ( safmin, MAX(f1,g1)));
+					gs_r = db_r/u;
+					gs_i = db_i/u;
+					g2 = sqrt ( gs_r*gs_r + gs_i*gs_i);
+					if (f1 /u < rtmin) {
+				v = MIN (safmax, MAX (safmin, f1));
+				w = v / u;
+				fs_r = *DA/ v;
+				fs_i = *(DA+1) / v;
+							f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
+							h = f2 * w * w + g2;		
+					} else { // use same scaling for both
+				w = 1.;
+				fs_r = *DA/ u;
+				fs_i = *(DA+1) / u;
+							f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
+				h = f2 + g2;
+					}
+					if ( f2 >= h * safmin) {
+						*C = sqrt ( f2 / h );
+						*DA = fs_r / *C;
+						*(DA+1) = fs_i / *C;
+						rtmax *= 2;
+						if ( f2 > rtmin && h < rtmax) {
+							*S = gs_r * (fs_r /sqrt(f2*h)) - gs_i * (fs_i / sqrt(f2*h));
+							*(S+1) = gs_r * (fs_i /sqrt(f2*h)) + gs_i * -1. * (fs_r / sqrt(f2*h));
+						} else {
+							*S = gs_r * (*DA/h) - gs_i * (*(DA+1) / h);
+							*(S+1) = gs_r * (*(DA+1) /h) + gs_i * -1. * (*DA / h);
+						}
+					} else { // intermediates might overflow
+				d = sqrt ( f2 * h);
+				*C = f2 /d;
+				if (*C >= safmin) {
+					*DA = fs_r / *C;
+					*(DA+1) = fs_i / *C;
+				} else {
+					*DA = fs_r * (h / d);
+					*(DA+1) = fs_i / (h / d);
+				}
+				*S = gs_r * (fs_r /d) - gs_i * (fs_i / d);
+				*(S+1) = gs_r * (fs_i /d) + gs_i * -1. * (fs_r / d);
+								}
+					*C *= w;
+					*DA *= u;
+					*(DA+1) *= u;
+					break;
+			}
+		}
+	} while (0);
+
+	FUNCTION_PROFILE_END(4, 4, 4);
+	IDEBUG_END;
 }
 	
@@ -139,6 +139,11 @@ void NAME(char *UPLO, blasint *N, blasint *K, FLOAT  *ALPHA, FLOAT *a, blasint *
  if (incy < 0 ) c -= (n - 1) * incy * COMPSIZE;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMPTEST
  nthreads = num_cpu_avail(2);
@@ -128,6 +128,11 @@ void NAME(char *UPLO, blasint *N, FLOAT  *ALPHA, FLOAT *a,
  if (incy < 0 ) c -= (n - 1) * incy * COMPSIZE;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMPTEST
  nthreads = num_cpu_avail(2);
@@ -120,6 +120,11 @@ void NAME(char *UPLO, blasint *N, FLOAT  *ALPHA,
  if (incx < 0 ) x -= (n - 1) * incx;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -123,6 +123,11 @@ void NAME(char *UPLO, blasint *N, FLOAT  *ALPHA,
  if (incy < 0 ) y -= (n - 1) * incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -118,6 +118,11 @@ void NAME(char *UPLO, blasint *N, FLOAT  *ALPHA, FLOAT *a, blasint *LDA,
  if (incy < 0 ) c -= (n - 1) * incy * COMPSIZE;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -205,6 +205,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, int n, void* valpha, FL
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -125,6 +125,11 @@ void NAME(char *UPLO, blasint *N, FLOAT  *ALPHA,
  if (incy < 0 ) y -= (n - 1) * incy;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -238,6 +238,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -211,6 +211,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  (tbsv[(trans<<2) | (uplo<<1) | unit])(n, k, a, lda, x, incx, buffer);

@@ -229,6 +229,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

 #ifdef SMP
  nthreads = num_cpu_avail(2);
@@ -202,6 +202,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  (tpsv[(trans<<2) | (uplo<<1) | unit])(n, a, x, incx, buffer);

@@ -208,6 +208,11 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo,
  if (incx < 0 ) x -= (n - 1) * incx * 2;

  buffer = (FLOAT *)blas_memory_alloc(1);
+  if (!buffer) {
+    info = -999;
+    BLASFUNC(xerbla)(ERROR_NAME, &info, sizeof(ERROR_NAME));
+    return;
+  }

  (trsv[(trans<<2) | (uplo<<1) | unit])(n, a, lda, x, incx, buffer);

@@ -290,6 +290,10 @@ function (build_core TARGET_CORE KDIR TSUFFIX KERNEL_DEFINITIONS)
   	  if (HAVE_SME)
 	  GenerateNamedObjects("${KERNELDIR}/${SGEMMDIRECTSMEKERNEL}" "" "gemm_direct_sme1_2VLx2VL" false "" "" false SINGLE)
 	  GenerateNamedObjects("${KERNELDIR}/${SGEMMDIRECTPREKERNEL}" "" "gemm_direct_sme1_preprocess" false "" "" false SINGLE)
+	  GenerateNamedObjects("${KERNELDIR}/sme_sgemm_kernel.c" "" "sme_*gemm_kernel" false "" "" false "SINGLE")
+	  GenerateNamedObjects("${KERNELDIR}/sme_dgemm_kernel.c" "" "sme_*gemm_kernel" false "" "" false "DOUBLE")
+	  GenerateNamedObjects("${KERNELDIR}/sme_cgemm_kernel.c" "" "sme_*gemm_kernel" false "" "" false "COMPLEX")
+	  GenerateNamedObjects("${KERNELDIR}/sme_zgemm_kernel.c" "" "sme_*gemm_kernel" false "" "" false "ZCOMPLEX")
 	  endif ()
   	  endif ()
    endif()
@@ -27,7 +27,7 @@ endif

 ifdef TARGET_CORE
 ifeq ($(TARGET_CORE), ARMV9SME)
- override CFLAGS += -DBUILD_KERNEL -DTABLE_NAME=gotoblas_$(TARGET_CORE) -march=armv9-a+sve2+sme
+ override CFLAGS += -DBUILD_KERNEL -DTABLE_NAME=gotoblas_$(TARGET_CORE) -march=armv9-a+sve2+sme+sme-f64f64
 ifdef OS_WINDOWS
 ifeq ($(C_COMPILER), CLANG)
   override CFLAGS += --aarch64-stack-hazard-size=0
@@ -38,7 +38,7 @@ ifeq ($(TARGET_CORE), VORTEXM4)
  ifeq ($(C_COMPILER), GCC)
    override CFLAGS += -DBUILD_KERNEL -DTABLE_NAME=gotoblas_$(TARGET_CORE) -UHAVE_SME -march=armv8.4-a
  else
-    override CFLAGS += -DBUILD_KERNEL -DTABLE_NAME=gotoblas_$(TARGET_CORE) -march=armv8.4-a+sme
+    override CFLAGS += -DBUILD_KERNEL -DTABLE_NAME=gotoblas_$(TARGET_CORE) -march=armv8.4-a+sme+sme-f64f64
 #   ifneq ($(APPLECLANG),1)
 #     override LDFLAGS += -lclang_rt_builtins-aarch64
 #   endif
@@ -264,7 +264,11 @@ SKERNELOBJS += \
 ifdef USE_SME
 SKERNELOBJS += \
 	sgemm_direct_sme1_2VLx2VL$(TSUFFIX).$(SUFFIX) \
-	sgemm_direct_sme1_preprocess$(TSUFFIX).$(SUFFIX)
+	sgemm_direct_sme1_preprocess$(TSUFFIX).$(SUFFIX) \
+	sme_sgemm_kernel$(TSUFFIX).$(SUFFIX) \
+	sme_dgemm_kernel$(TSUFFIX).$(SUFFIX) \
+	sme_cgemm_kernel$(TSUFFIX).$(SUFFIX) \
+	sme_zgemm_kernel$(TSUFFIX).$(SUFFIX)
 endif
 endif
 endif
@@ -1064,6 +1068,14 @@ $(KDIR)sgemm_direct_sme1_2VLx2VL$(TSUFFIX).$(SUFFIX) :
 	$(CC) $(CFLAGS) -c $(KERNELDIR)/sgemm_direct_sme1_2VLx2VL.S -UDOUBLE -UCOMPLEX  -o $@
 $(KDIR)sgemm_direct_sme1_preprocess$(TSUFFIX).$(SUFFIX) : 
 	$(CC) $(CFLAGS) -c $(KERNELDIR)/sgemm_direct_sme1_preprocess.S -UDOUBLE -UCOMPLEX  -o $@
+$(KDIR)sme_sgemm_kernel$(TSUFFIX).$(SUFFIX) : 
+	$(CC) $(CFLAGS) -c $(KERNELDIR)/sme_sgemm_kernel.c -UDOUBLE -UCOMPLEX  -o $@
+$(KDIR)sme_dgemm_kernel$(TSUFFIX).$(SUFFIX) : 
+	$(CC) $(CFLAGS) -c $(KERNELDIR)/sme_dgemm_kernel.c -DDOUBLE -UCOMPLEX  -o $@
+$(KDIR)sme_cgemm_kernel$(TSUFFIX).$(SUFFIX) : 
+	$(CC) $(CFLAGS) -c $(KERNELDIR)/sme_cgemm_kernel.c -UDOUBLE -DCOMPLEX  -o $@
+$(KDIR)sme_zgemm_kernel$(TSUFFIX).$(SUFFIX) : 
+	$(CC) $(CFLAGS) -c $(KERNELDIR)/sme_zgemm_kernel.c -DDOUBLE -DCOMPLEX  -o $@
 endif
 endif
 endif
@@ -1,3 +1,28 @@
 include $(KERNELDIR)/KERNEL.CORTEXA57

+# Cortex-A72 DGEMM uses a dedicated 6x8 microkernel (see
+# dgemm_kernel_6x8_cortexa72.S). Stock TARGET=CORTEXA72 previously
+# aliased CORTEXA57 (8x4). Keep SGEMM/CGEMM/ZGEMM from A57 for now.
+#
+# Stock generic/gemm_{t,n}copy_6.c packs as 4+2, not contiguous MR=6 —
+# use the panel packers. NR=8 reuses the asm dgemm_{n,t}copy_8.S.
+# No dtrmm_kernel_6x8.S yet; fall back to generic 2x2 for DTRMM.
+# Stock generic/trsm_kernel_*.c assumes power-of-2 UNROLL_M; use the
+# UNROLLM6-aware kernels (same idea as loongarch64 UNROLLN6).
+DGEMMKERNEL    =  dgemm_kernel_6x8_cortexa72.S
+DTRMMKERNEL    =  ../generic/trmmkernel_2x2.c

+DTRSMKERNEL_LN =  trsm_kernel_LN_UNROLLM6.c
+DTRSMKERNEL_LT =  trsm_kernel_LT_UNROLLM6.c
+DTRSMKERNEL_RN =  trsm_kernel_RN_UNROLLM6.c
+DTRSMKERNEL_RT =  trsm_kernel_RT_UNROLLM6.c
+
+DGEMMINCOPY    =  gemm_ncopy_6_panel.c
+DGEMMITCOPY    =  gemm_tcopy_6_panel.c
+DGEMMINCOPYOBJ =  dgemm_incopy$(TSUFFIX).$(SUFFIX)
+DGEMMITCOPYOBJ =  dgemm_itcopy$(TSUFFIX).$(SUFFIX)
+
+DGEMMONCOPY    =  dgemm_ncopy_8.S
+DGEMMOTCOPY    =  dgemm_tcopy_8.S
+DGEMMONCOPYOBJ =  dgemm_oncopy$(TSUFFIX).$(SUFFIX)
+DGEMMOTCOPYOBJ =  dgemm_otcopy$(TSUFFIX).$(SUFFIX)
@@ -168,6 +168,11 @@ SCABS_KERNEL	= ../generic/cabs.c
 DCABS_KERNEL	= ../generic/cabs.c
 QCABS_KERNEL	= ../generic/cabs.c

+SOMATCOPY_CT = ../arm/omatcopy_ct.c
+SOMATCOPY_RT = ../arm/omatcopy_rt.c
+DOMATCOPY_CT = ../arm/omatcopy_ct.c
+DOMATCOPY_RT = ../arm/omatcopy_rt.c
+
 #Dump kernel
 CGEMM3MKERNEL    = ../generic/zgemm3mkernel_dump.c
 ZGEMM3MKERNEL    = ../generic/zgemm3mkernel_dump.c
@@ -0,0 +1,94 @@
+/* True MR=6 no-transpose pack (INCOPY) for OpenBLAS DGEMM 6xN.
+ *
+ * Used when A is transposed (TN/TT/...): INCOPY(m=Kc, n=Mc).
+ * Packs 6 columns at a time into contiguous 6-wide micropanels along m
+ * (same layout as loongarch64/gemm_ncopy_6.prefx.c).
+ *
+ * Stock generic/gemm_ncopy_6.c packs as 4-wide and is wrong for 6x8.
+ */
+#include "common.h"
+
+int CNAME(BLASLONG m, BLASLONG n, FLOAT *a, BLASLONG lda, FLOAT *b)
+{
+  BLASLONG i, j;
+  FLOAT *aoffset, *aoffset1, *aoffset2, *aoffset3, *aoffset4, *aoffset5, *aoffset6;
+  FLOAT *boffset;
+  FLOAT c1, c2, c3, c4, c5, c6;
+
+  aoffset = a;
+  boffset = b;
+
+  j = n / 6;
+  if (j > 0) {
+    do {
+      aoffset1 = aoffset;
+      aoffset2 = aoffset1 + lda;
+      aoffset3 = aoffset2 + lda;
+      aoffset4 = aoffset3 + lda;
+      aoffset5 = aoffset4 + lda;
+      aoffset6 = aoffset5 + lda;
+      aoffset += 6 * lda;
+
+      i = m;
+      if (i > 0) {
+        do {
+          c1 = *(aoffset1);
+          c2 = *(aoffset2);
+          c3 = *(aoffset3);
+          c4 = *(aoffset4);
+          c5 = *(aoffset5);
+          c6 = *(aoffset6);
+          aoffset1++; aoffset2++; aoffset3++;
+          aoffset4++; aoffset5++; aoffset6++;
+          *(boffset + 0) = c1;
+          *(boffset + 1) = c2;
+          *(boffset + 2) = c3;
+          *(boffset + 3) = c4;
+          *(boffset + 4) = c5;
+          *(boffset + 5) = c6;
+          boffset += 6;
+          i--;
+        } while (i > 0);
+      }
+      j--;
+    } while (j > 0);
+  }
+
+  /* n remainder: 4, 2, 1 column groups (matches ukernel N-edge order for A-side) */
+  {
+    BLASLONG nr = n - (n / 6) * 6;
+    if (nr >= 4) {
+      aoffset1 = aoffset;
+      aoffset2 = aoffset1 + lda;
+      aoffset3 = aoffset2 + lda;
+      aoffset4 = aoffset3 + lda;
+      aoffset += 4 * lda;
+      for (i = 0; i < m; i++) {
+        boffset[0] = *aoffset1++;
+        boffset[1] = *aoffset2++;
+        boffset[2] = *aoffset3++;
+        boffset[3] = *aoffset4++;
+        boffset += 4;
+      }
+      nr -= 4;
+    }
+    if (nr >= 2) {
+      aoffset1 = aoffset;
+      aoffset2 = aoffset1 + lda;
+      aoffset += 2 * lda;
+      for (i = 0; i < m; i++) {
+        boffset[0] = *aoffset1++;
+        boffset[1] = *aoffset2++;
+        boffset += 2;
+      }
+      nr -= 2;
+    }
+    if (nr >= 1) {
+      aoffset1 = aoffset;
+      for (i = 0; i < m; i++) {
+        *boffset++ = *aoffset1++;
+      }
+    }
+  }
+  return 0;
+}
@@ -0,0 +1,63 @@
+/* True MR=6 transpose pack (ITCOPY) for OpenBLAS DGEMM 6xN.
+ *
+ * Call signature matches gemm_tcopy_*: (m, n, a, lda, b) where for NN
+ * ICOPY uses ITCOPY(m=Kc, n=Mc). Output is Goto-style micropanels:
+ *   for each panel of 6 rows: Kc contiguous packs of 6 doubles.
+ * Remainders after full panels: 4, then 2, then 1 (kernel cascade).
+ *
+ * Stock generic/gemm_tcopy_6.c is a clone of tcopy_4 (4+2) and does NOT
+ * match a contiguous 6-wide ukernel — do not use it for 6x8.
+ */
+#include "common.h"
+
+int CNAME(BLASLONG m, BLASLONG n, FLOAT *a, BLASLONG lda, FLOAT *b)
+{
+  BLASLONG j, js;
+  FLOAT *b_ptr = b;
+  BLASLONG n6 = n / 6;
+  BLASLONG nr = n - n6 * 6;
+
+  for (js = 0; js < n6; js++) {
+    FLOAT *a_row = a + js * 6;
+    for (j = 0; j < m; j++) {
+      FLOAT *ap = a_row + j * lda;
+      b_ptr[0] = ap[0];
+      b_ptr[1] = ap[1];
+      b_ptr[2] = ap[2];
+      b_ptr[3] = ap[3];
+      b_ptr[4] = ap[4];
+      b_ptr[5] = ap[5];
+      b_ptr += 6;
+    }
+  }
+
+  FLOAT *a_row = a + n6 * 6;
+  if (nr >= 4) {
+    for (j = 0; j < m; j++) {
+      FLOAT *ap = a_row + j * lda;
+      b_ptr[0] = ap[0];
+      b_ptr[1] = ap[1];
+      b_ptr[2] = ap[2];
+      b_ptr[3] = ap[3];
+      b_ptr += 4;
+    }
+    a_row += 4;
+    nr -= 4;
+  }
+  if (nr >= 2) {
+    for (j = 0; j < m; j++) {
+      FLOAT *ap = a_row + j * lda;
+      b_ptr[0] = ap[0];
+      b_ptr[1] = ap[1];
+      b_ptr += 2;
+    }
+    a_row += 2;
+    nr -= 2;
+  }
+  if (nr >= 1) {
+    for (j = 0; j < m; j++) {
+      *b_ptr++ = a_row[j * lda];
+    }
+  }
+  return 0;
+}
@@ -0,0 +1,544 @@
+//#include <algorithm>
+#include <arm_neon.h>
+//#include <cstddef>
+#include <stddef.h>
+#include <stdalign.h>
+#include <stdbool.h>
+#include "common.h"
+#ifndef stdmin
+#define stdmin(a,b)   (a>b? b:a)
+#endif
+typedef float _Complex cfloat;
+static cfloat CMUL(cfloat a, cfloat b,bool conja, bool conjb) {
+float ra=creal(a);
+float rb=creal(b);
+float ia=conja ? -cimag(a) : cimag(a);
+float ib=conjb ? -cimag(b) : cimag(b);
+float r1=ra*rb;
+float r2=ia*ib;
+float r=r1-r2;
+float i=(ra+ia)*(rb+ib)-r1-r2;
+cfloat res={r,i};
+return res;
+}
+#define KERNEL_ALPHA 0
+#define USE_VECTORIZED_PACKING 1
+
+
+static cfloat czero={0.,0.};
+static cfloat cone={1.,0.};
+
+#define MC 256
+#define KC 2048
+#define NC 1024
+
+static inline void cgemm_sme_compute_16x16_tile(blasint current_K, const float *A_ptr, const float *B_ptr, cfloat *C_ptr, size_t ldc, int beta_mode, const cfloat *beta_ptr)
+{
+    size_t ldc_bytes = ldc * sizeof(cfloat);
+
+    asm volatile("smstart\n\t"
+		 "ptrue p0.s\n\t"
+                 "zero {za}\n\t"
+                 "cmp %w[beta_mode], #0\n\t"
+                 "b.eq 19f\n\t"
+
+                 // FIX 1: Correct SVE mnemonics for loading 32-bit floats (Real/Imag)
+                 "cmp %w[beta_mode], #1\n\t"
+                 "b.ne 10f\n\t"
+                 "ld1rw z30.s, p0/z, [%[beta_ptr]]\n\t" // Load beta.real()
+                 "add x15, %[beta_ptr], #4\n\t"         // 4-byte offset for float
+                 "ld1rw z31.s, p0/z, [x15]\n\t"         // Load beta.imag()
+                 "10:\n\t"
+
+                 "mov w12, #0\n\t"
+                 "mov x13, %[c]\n\t"
+                 "100:\n\t"
+                 "ld1w z0.s, p0/z, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "ld1w z1.s, p0/z, [x14]\n\t"
+                 "uzp1 z2.s, z0.s, z1.s\n\t" // z2 = C_re
+                 "uzp2 z3.s, z0.s, z1.s\n\t" // z3 = C_im
+
+                 "cmp %w[beta_mode], #1\n\t"
+                 "b.ne 101f\n\t"
+
+                 // FIX 2: Fully implemented Complex Beta Multiplication using movprfx
+                 "movprfx z4, z2\n\t"
+                 "fmul z4.s, p0/m, z4.s, z30.s\n\t" // z4 = Cre * Bre
+                 "movprfx z5, z3\n\t"
+                 "fmul z5.s, p0/m, z5.s, z31.s\n\t" // z5 = Cim * Bim
+                 "movprfx z6, z4\n\t"
+                 "fsub z6.s, p0/m, z6.s, z5.s\n\t" // z6 = Cre' (Cre*Bre - Cim*Bim)
+
+                 "movprfx z8, z2\n\t"
+                 "fmul z8.s, p0/m, z8.s, z31.s\n\t" // z8 = Cre * Bim
+                 "movprfx z9, z3\n\t"
+                 "fmul z9.s, p0/m, z9.s, z30.s\n\t" // z9 = Cim * Bre
+                 "movprfx z7, z8\n\t"
+                 "fadd z7.s, p0/m, z7.s, z9.s\n\t" // z7 = Cim' (Cre*Bim + Cim*Bre)
+
+                 "mova za0v.s[w12, 0], p0/m, z6.s\n\t"
+                 "mova za2v.s[w12, 0], p0/m, z7.s\n\t"
+                 "b 102f\n\t"
+
+                 "101:\n\t" // Fallback: beta == 1.0
+                 "mova za0v.s[w12, 0], p0/m, z2.s\n\t"
+                 "mova za2v.s[w12, 0], p0/m, z3.s\n\t"
+
+                 "102:\n\t"
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w12, w12, #1\n\t"
+                 "cmp w12, #16\n\t"
+                 "b.ne 100b\n\t"
+
+                 "19:\n\t"
+                 "mov w10, %w[k]\n\t"
+                 "cbz w10, 3f\n\t"
+                 "11:\n\t"
+                 "ld1w z0.s, p0/z, [%[a], #0, mul vl]\n\t"
+                 "ld1w z1.s, p0/z, [%[a], #1, mul vl]\n\t"
+                 "ld1w z2.s, p0/z, [%[b], #0, mul vl]\n\t"
+                 "ld1w z3.s, p0/z, [%[b], #1, mul vl]\n\t"
+                 "fmopa za0.s, p0/m, p0/m, z0.s, z2.s\n\t"
+                 "fmopa za1.s, p0/m, p0/m, z1.s, z3.s\n\t"
+                 "fmopa za2.s, p0/m, p0/m, z0.s, z3.s\n\t"
+                 "fmopa za3.s, p0/m, p0/m, z1.s, z2.s\n\t"
+                 "add %[a], %[a], #128\n\t"
+                 "add %[b], %[b], #128\n\t"
+                 "subs w10, w10, #1\n\t"
+                 "b.ne 11b\n\t"
+
+                 "3:\n\t"
+                 "mov w12, #0\n\t"
+                 "mov x13, %[c]\n\t"
+                 "200:\n\t"
+                 "mova z0.s, p0/m, za0v.s[w12, 0]\n\t"
+                 "mova z1.s, p0/m, za1v.s[w12, 0]\n\t"
+                 "mova z2.s, p0/m, za2v.s[w12, 0]\n\t"
+                 "mova z3.s, p0/m, za3v.s[w12, 0]\n\t"
+
+                 // FIX 3: Non-destructive SVE arithmetic on store
+                 "movprfx z4, z0\n\t"
+                 "fsub z4.s, p0/m, z4.s, z1.s\n\t"
+                 "movprfx z5, z2\n\t"
+                 "fadd z5.s, p0/m, z5.s, z3.s\n\t"
+
+                 "zip1 z6.s, z4.s, z5.s\n\t"
+                 "zip2 z7.s, z4.s, z5.s\n\t"
+                 "st1w z6.s, p0, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "st1w z7.s, p0, [x14]\n\t"
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w12, w12, #1\n\t"
+                 "cmp w12, #16\n\t"
+                 "b.ne 200b\n\t"
+		 "smstop\n\t"
+                 "msr fpsr, xzr\n\t"
+    : [a] "+&r"(A_ptr), [b] "+&r"(B_ptr)
+    : [k] "r"(current_K), [c] "r"(C_ptr), [ldc_bytes] "r"(ldc_bytes), [beta_mode] "r"(beta_mode), [beta_ptr] "r"(beta_ptr)
+    // Updated Clobber list to cover x15 and z8-z9, z30-z31
+    :  "x10", "x12", "x13", "x14", "x15", "cc", "memory",
+       "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7",
+       "v8", "v9", "v10", "v11", "v12", "v13", "v14", "v15",
+       "v16", "v17", "v18", "v19", "v20", "v21", "v22", "v23",
+       "v24", "v25", "v26", "v27", "v28", "v29", "v30", "v31",
+       "z0", "z1", "z2", "z3", "z4", "z5", "z6", "z7",
+       "z8", "z9", "z10", "z11", "z12", "z13", "z14", "z15",
+       "z16", "z17", "z18", "z19", "z20", "z21", "z22", "z23",
+       "z24", "z25", "z26", "z27", "z28", "z29", "z30", "z31", "za",
+       "p0", "p1", "p2", "p3", "p4", "p5", "p6", "p7",
+       "p8", "p9", "p10", "p11", "p12", "p13", "p14", "p15");
+
+
+}
+
+static void cgemm_sme_NN(int M, int N, int K, const cfloat alpha, const cfloat *A, int lda, const cfloat *b, int ldb, const cfloat beta, cfloat *C, int ldc, bool conja, bool conjb)
+{
+
+    if (alpha == czero || K == 0) {
+        if (beta == czero) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = czero;
+                }
+            }
+        } else {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] *= beta;
+                }
+            }
+	}
+        return;
+    }
+
+
+
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC * 2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 15) & ~15;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta==czero) ? 0 : (beta != cone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 16) {
+                float *__restrict B_out = &B_pack[(jj / 16) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    for (int bc = 0; bc < 16; ++bc) {
+                        cfloat val = (jj + bc < current_N) ? b[(j + jj + bc) * ldb + k + kk] : czero;
+                        B_out[kk * 32 + bc] = creal(val);
+                        B_out[kk * 32 + 16 + bc] = conjb ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 15) & ~15;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 16) {
+                    float *__restrict A_out = &A_pack[(ii / 16) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        for (int br = 0; br < 16; ++br) {
+                            cfloat val = (ii + br < current_M) ? CMUL(A[(k + kk) * lda + i + ii + br] , alpha, conja, 0) : czero;
+                            A_out[kk * 32 + br] = creal(val);
+                            A_out[kk * 32 + 16 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 16) {
+                    int current_N_block = stdmin(16, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 16) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 16) {
+                        int current_M_block = stdmin(16, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 16) * panel_stride_A];
+                        cfloat *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 16 && current_N_block == 16) {
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) cfloat C_buffer[256];
+                            if (beta_mode == 0) {
+                                for (int idx = 0; idx < 256; ++idx) {
+                                    C_buffer[idx] =  czero;
+                                }
+                            }
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 16] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_buffer, 16, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 16];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void cgemm_sme_TN(int M, int N, int K, const cfloat alpha, const cfloat *A, int lda, const cfloat *b, int ldb, const cfloat beta, cfloat *C, int ldc, bool conja, bool conjb)
+{
+    if (alpha == czero || K == 0) {
+        if (beta == czero) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = czero;
+                }
+            }
+        } else {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = CMUL(C[i+j*ldc],beta,0,0);
+                }
+            }
+	}
+        return;
+    }
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC * 2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 15) & ~15;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta==czero) ? 0 : (beta != cone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 16) {
+                float *__restrict B_out = &B_pack[(jj / 16) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    for (int bc = 0; bc < 16; ++bc) {
+                        cfloat val = (jj + bc < current_N) ? b[(j + jj + bc) * ldb + k + kk] : czero;
+                        B_out[kk * 32 + bc] = creal(val);
+                        B_out[kk * 32 + 16 + bc] = conjb ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 15) & ~15;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 16) {
+                    float *__restrict A_out = &A_pack[(ii / 16) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        for (int br = 0; br < 16; ++br) {
+                            cfloat val = (ii + br < current_M) ? CMUL(A[(i + ii + br) * lda + k + kk],alpha,conja,0) : czero;
+                            A_out[kk * 32 + br] = creal(val);
+                            A_out[kk * 32 + 16 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 16) {
+                    int current_N_block = stdmin(16, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 16) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 16) {
+                        int current_M_block = stdmin(16, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 16) * panel_stride_A];
+                        cfloat *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 16 && current_N_block == 16) {
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) cfloat C_buffer[256];
+                                for (int idx = 0; idx < 256; ++idx) {
+                                    C_buffer[idx] =  czero;
+                                }
+                            if (beta_mode != 0) {
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 16] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            
+                            }
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_buffer, 16, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 16];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void cgemm_sme_NT(int M, int N, int K, const cfloat alpha, const cfloat *A, int lda, const cfloat *b, int ldb, const cfloat beta, cfloat *C, int ldc, bool conja, bool conjb)
+{
+    if (alpha == czero || K == 0) {
+        if (beta == czero) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = czero;
+                }
+            }
+        } else {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] *= beta;
+                }
+            }
+	}
+        return;
+    }
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC * 2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 15) & ~15;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta==czero) ? 0 : (beta != cone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 16) {
+                float *__restrict B_out = &B_pack[(jj / 16) * panel_stride_B];
+                for (int kk = 0; kk < current_K; kk++) {
+                    for (int bc = 0; bc < 16; ++bc) {
+                        cfloat val = (jj + bc < current_N) ? b[(k + kk) * ldb + j + jj + bc] : czero ;
+                        B_out[kk * 32 + bc] = creal(val);
+                        B_out[kk * 32 + 16 + bc] = (conjb) ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 15) & ~15;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 16) {
+                    float *__restrict A_out = &A_pack[(ii / 16) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; kk++) {
+                        for (int br = 0; br < 16; ++br) {
+                            cfloat val = (ii + br < current_M) ? CMUL(A[(k + kk) * lda + i + ii + br],alpha,conja,0) :  czero;
+                            A_out[kk * 32 + br] = creal(val);
+                            A_out[kk * 32 + 16 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 16) {
+                    int current_N_block = stdmin(16, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 16) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 16) {
+                        int current_M_block = stdmin(16, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 16) * panel_stride_A];
+                        cfloat *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 16 && current_N_block == 16) {
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) cfloat C_buffer[256];
+                            if (beta_mode != 0) {
+                                for (int idx = 0; idx < 256; ++idx) {
+                                    C_buffer[idx] =  czero;
+                                }
+                            }
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 16] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                        
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_buffer, 16, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 16];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void cgemm_sme_TT(int M, int N, int K, const cfloat alpha, const cfloat *A, int lda, const cfloat *b, int ldb, const cfloat beta, cfloat *C, int ldc, bool conja, bool conjb)
+{
+    if (alpha == czero || K == 0) {
+        if (beta == czero) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = czero;
+                }
+            }
+        } else {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] *= beta;
+                }
+            }
+	}
+        return;
+    }
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC * 2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 15) & ~15;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == czero) ? 0 : (beta != cone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 16) {
+                float *__restrict B_out = &B_pack[(jj / 16) * panel_stride_B];
+                for (int kk = 0; kk < current_K; kk++) {
+                    for (int bc = 0; bc < 16; ++bc) {
+                        cfloat val = (jj + bc < current_N) ? b[(k + kk) * ldb + j + jj + bc] : czero;
+                        B_out[kk * 32 + bc] = creal(val);
+                        B_out[kk * 32 + 16 + bc] = (conjb) ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 15) & ~15;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 16) {
+                    float *__restrict A_out = &A_pack[(ii / 16) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; kk++) {
+                        for (int br = 0; br < 16; ++br) {
+                            cfloat val = (ii + br < current_M) ? CMUL( A[(i + ii + br) * lda + k + kk],alpha,conja,0) : czero;
+                            A_out[kk * 32 + br] = creal(val);
+                            A_out[kk * 32 + 16 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 16) {
+                    int current_N_block = stdmin(16, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 16) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 16) {
+                        int current_M_block = stdmin(16, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 16) * panel_stride_A];
+                        cfloat *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 16 && current_N_block == 16) {
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) cfloat C_buffer[256];
+                            if (beta_mode != 0) {
+                                for (int idx = 0; idx < 256; ++idx) {
+                                    C_buffer[idx] =  czero;
+                                }
+			    }
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 16] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            
+                            cgemm_sme_compute_16x16_tile(current_K, A_ptr, B_ptr, C_buffer, 16, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 16];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+void CNAME(const char *transa, const char *transb, const BLASLONG m, const BLASLONG n, const BLASLONG k, const float alpha_r, const float alpha_i, const float *a, const BLASLONG lda, const float *b, const BLASLONG ldb, const float beta_r, const float beta_i, float *c, const BLASLONG ldc)
+{
+cfloat alpha={alpha_r,alpha_i};
+cfloat beta={beta_r,beta_i};
+    bool trans_a = (*transa == 'T' || *transa == 't' || *transa == 'C' || *transa == 'c');
+    bool trans_b = (*transb == 'T' || *transb == 't' || *transb == 'C' || *transb == 'c');
+    if (!trans_a && !trans_b) {
+    bool conja=(*transa == 'R' || *transa == 'r');
+    bool conjb=(*transb == 'R' || *transb == 'r');
+        cgemm_sme_NN(m, n, k, alpha, (const cfloat*) a, lda, (const cfloat*) b, ldb, beta, (cfloat*)c, ldc, conja, conjb);
+    }
+    else if (trans_a && !trans_b) {
+    bool conja=(*transa == 'C' || *transa == 'c');
+    bool conjb=(*transb == 'R' || *transb == 'r');
+        cgemm_sme_TN(m, n, k, alpha, (const cfloat*) a, lda, (const cfloat*) b, ldb, beta, (cfloat*)c, ldc, conja, conjb);
+    }
+    else if (!trans_a && trans_b) {
+    bool conja=(*transa == 'R' || *transa == 'r');
+    bool conjb=(*transb == 'C' || *transb == 'c');
+        cgemm_sme_NT(m, n, k, alpha, (const cfloat*) a, lda, (const cfloat*) b, ldb, beta, (cfloat*)c, ldc, conja, conjb);
+    }
+    else {
+    bool conja=(*transa == 'C' || *transa == 'c');
+    bool conjb=(*transb == 'C' || *transb == 'c');
+        cgemm_sme_TT(m, n, k, alpha, (const cfloat*) a, lda, (const cfloat*) b, ldb, beta, (cfloat*)c, ldc, conja, conjb);
+    }
+}
@@ -0,0 +1,633 @@
+//#include <algorithm>
+#include <arm_neon.h>
+//#include <cstddef>
+#include <stddef.h>
+#include <stdalign.h>
+#include <stdbool.h>
+#include "common.h"
+#ifndef stdmin
+#define stdmin(a,b)   (a>b? b:a)
+#endif
+
+#define MC 512
+#define KC 1024
+#define NC 2048
+
+static inline void sgemm_sme_compute_32x32_tile(blasint current_K, const float *A_ptr, const float *B_ptr, float *C_ptr, size_t ldc, blasint beta_mode, const float *beta_ptr)
+{
+    size_t ldc_bytes = ldc * sizeof(float);
+
+    asm volatile("smstart\n\t"
+		 "ptrue p0.s\n\t"
+                 "cmp %w[beta_mode], #0\n\t"
+                 "b.eq 10f\n\t"
+                 "cmp %w[beta_mode], #1\n\t"
+                 "b.eq 11f\n\t"
+                 "b 12f\n\t"
+
+                 "10:\n\t"
+                 "zero {za}\n\t"
+                 "b 19f\n\t"
+
+                 "11:\n\t"
+                 "ld1rw z31.s, p0/z, [%[beta_ptr]]\n\t"
+                 "mov w12, #0\n\t"
+                 "mov x13, %[c]\n\t"
+                 "110:\n\t"
+                 "ld1w z0.s, p0/z, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "ld1w z1.s, p0/z, [x14]\n\t"
+                 "fmul z0.s, p0/m, z0.s, z31.s\n\t"
+                 "fmul z1.s, p0/m, z1.s, z31.s\n\t"
+                 "mova za0v.s[w12, 0], p0/m, z0.s\n\t"
+                 "mova za2v.s[w12, 0], p0/m, z1.s\n\t" // FIXED: za2 (Bottom-Left)
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w12, w12, #1\n\t"
+                 "cmp w12, #16\n\t"
+                 "b.ne 110b\n\t"
+
+                 "mov w15, #0\n\t"
+                 "111:\n\t"
+                 "ld1w z0.s, p0/z, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "ld1w z1.s, p0/z, [x14]\n\t"
+                 "fmul z0.s, p0/m, z0.s, z31.s\n\t"
+                 "fmul z1.s, p0/m, z1.s, z31.s\n\t"
+                 "mova za1v.s[w15, 0], p0/m, z0.s\n\t" // FIXED: za1 (Top-Right)
+                 "mova za3v.s[w15, 0], p0/m, z1.s\n\t"
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w15, w15, #1\n\t"
+                 "cmp w15, #16\n\t"
+                 "b.ne 111b\n\t"
+                 "b 19f\n\t"
+
+                 "12:\n\t"
+                 "mov w12, #0\n\t"
+                 "mov x13, %[c]\n\t"
+                 "100:\n\t"
+                 "ld1w {za0v.s[w12, 0]}, p0/z, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "ld1w {za2v.s[w12, 0]}, p0/z, [x14]\n\t" // FIXED: za2
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w12, w12, #1\n\t"
+                 "cmp w12, #16\n\t"
+                 "b.ne 100b\n\t"
+
+                 "mov w15, #0\n\t"
+                 "101:\n\t"
+                 "ld1w {za1v.s[w15, 0]}, p0/z, [x13]\n\t" // FIXED: za1
+                 "add x14, x13, #64\n\t"
+                 "ld1w {za3v.s[w15, 0]}, p0/z, [x14]\n\t"
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w15, w15, #1\n\t"
+                 "cmp w15, #16\n\t"
+                 "b.ne 101b\n\t"
+
+                 "19:\n\t"
+                 "mov w10, %w[k]\n\t"
+                 "cbz w10, 3f\n\t"
+                 "11:\n\t"
+                 "ld1w z0.s, p0/z, [%[a], #0, mul vl]\n\t"
+                 "ld1w z1.s, p0/z, [%[a], #1, mul vl]\n\t"
+                 "ld1w z2.s, p0/z, [%[b], #0, mul vl]\n\t"
+                 "ld1w z3.s, p0/z, [%[b], #1, mul vl]\n\t"
+                 "fmopa za0.s, p0/m, p0/m, z0.s, z2.s\n\t"
+                 "fmopa za1.s, p0/m, p0/m, z0.s, z3.s\n\t"
+                 "fmopa za2.s, p0/m, p0/m, z1.s, z2.s\n\t"
+                 "fmopa za3.s, p0/m, p0/m, z1.s, z3.s\n\t"
+                 "add %[a], %[a], #128\n\t"
+                 "add %[b], %[b], #128\n\t"
+                 "subs w10, w10, #1\n\t"
+                 "b.ne 11b\n\t"
+
+                 "3:\n\t"
+                 "mov w12, #0\n\t"
+                 "mov x13, %[c]\n\t"
+                 "200:\n\t"
+                 "st1w {za0v.s[w12, 0]}, p0, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "st1w {za2v.s[w12, 0]}, p0, [x14]\n\t" // FIXED: za2
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w12, w12, #1\n\t"
+                 "cmp w12, #16\n\t"
+                 "b.ne 200b\n\t"
+
+                 "mov w15, #0\n\t"
+                 "201:\n\t"
+                 "st1w {za1v.s[w15, 0]}, p0, [x13]\n\t" // FIXED: za1
+                 "add x14, x13, #64\n\t"
+                 "st1w {za3v.s[w15, 0]}, p0, [x14]\n\t"
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w15, w15, #1\n\t"
+                 "cmp w15, #16\n\t"
+                 "b.ne 201b\n\t"
+		 "smstop\n\t"
+		"msr fpsr, xzr\n\t"
+    : [a] "+&r"(A_ptr), [b] "+&r"(B_ptr)
+    : [k] "r"(current_K), [c] "r"(C_ptr), [ldc_bytes] "r"(ldc_bytes), [beta_mode] "r"(beta_mode), [beta_ptr] "r"(beta_ptr)
+    : "x10", "x12", "x13", "x14", "x15", "memory", "cc", 
+                         "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7",
+                         "v8", "v9", "v10", "v11", "v12", "v13", "v14", "v15",
+                         "v16", "v17", "v18", "v19", "v20", "v21", "v22", "v23",
+                         "v24", "v25", "v26", "v27", "v28", "v29", "v30", "v31",
+                         "z0", "z1", "z2", "z3", "z4", "z5", "z6", "z7",
+                         "z8", "z9", "z10", "z11", "z12", "z13", "z14", "z15",
+                         "z16", "z17", "z18", "z19", "z20", "z21", "z22", "z23",
+                         "z24", "z25", "z26", "z27", "z28", "z29", "z30", "z31",
+      			 "p0", "p1", "p2", "p3", "p4", "p5", "p6", "p7",
+                         "p8", "p9", "p10", "p11", "p12", "p13", "p14", "p15","za");
+
+
+}
+
+static void sgemm_sme_NN(blasint M, blasint N, blasint K, float alpha, const float *A, blasint lda, const float *B, blasint ldb, float beta, float *C, blasint ldc)
+{
+    if (alpha == 0.0f || K == 0) {
+        if (beta != 1.0f) {
+            for (blasint j = 0; j < N; ++j) {
+                for (blasint i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == 0.0f) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC];
+
+#pragma omp parallel for schedule(dynamic, 1)
+    for (blasint j = 0; j < N; j += NC) {
+        blasint current_N = stdmin(NC, N - j);
+        blasint N_pad = (current_N + 31) & ~31;
+        for (blasint k = 0; k < K; k += KC) {
+            blasint current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == 0.0f) ? 0 : (beta != 1.0f ? 1 : 2)) : 2;
+
+            blasint N_main = current_N & ~31;
+            for (blasint jj = 0; jj < N_main; jj += 32) {
+                float *__restrict B_out = &B_pack[(jj / 32) * panel_stride_B];
+                for (blasint kk = 0; kk < current_K; ++kk) {
+                    const float *__restrict B_in = &B[(j + jj) * (size_t)ldb + k + kk];
+                    for (int bc = 0; bc < 32; ++bc) {
+                        B_out[kk * 32 + bc] = B_in[bc * ldb];
+                    }
+                }
+            }
+            if (N_main < current_N) {
+                float *__restrict B_out = &B_pack[(N_main / 32) * panel_stride_B];
+                for (int bc = 0; bc < 32; ++bc) {
+                    if (N_main + bc < current_N) {
+                        const float *__restrict B_col = &B[(j + N_main + bc) * (size_t)ldb + k];
+                        for (blasint kk = 0; kk < current_K; ++kk) {
+                            B_out[kk * 32 + bc] = B_col[kk];
+                        }
+                    }
+                    else {
+                        for (blasint kk = 0; kk < current_K; ++kk) {
+                            B_out[kk * 32 + bc] = 0.0f;
+                        }
+                    }
+                }
+            }
+
+            for (blasint i = 0; i < M; i += MC) {
+                blasint current_M = stdmin(MC, M - i);
+                blasint M_pad = (current_M + 31) & ~31;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+                float32x4_t valpha = vdupq_n_f32(alpha);
+
+                blasint M_main = current_M & ~31;
+                for (blasint ii = 0; ii < M_main; ii += 32) {
+                    float *__restrict A_out = &A_pack[(ii / 32) * panel_stride_A];
+                    for (blasint kk = 0; kk < current_K; ++kk) {
+                        const float *__restrict A_col = &A[(k + kk) * (size_t)lda + i + ii];
+                        float *__restrict A_out_ptr = &A_out[kk * 32];
+                        for (int v = 0; v < 8; ++v) {
+                            vst1q_f32(&A_out_ptr[v * 4], vmulq_f32(vld1q_f32(&A_col[v * 4]), valpha));
+                        }
+                    }
+                }
+                if (M_main < current_M) {
+                    float *__restrict A_out = &A_pack[(M_main / 32) * panel_stride_A];
+                    for (blasint kk = 0; kk < current_K; ++kk) {
+                        const float *__restrict A_col = &A[(k + kk) * (size_t)lda + i + M_main];
+                        for (blasint br = 0; br < current_M - M_main; ++br) {
+                            A_out[kk * 32 + br] = A_col[br] * alpha;
+                        }
+                        for (int br = current_M - M_main; br < 32; ++br) {
+                            A_out[kk * 32 + br] = 0.0f;
+                        }
+                    }
+                }
+
+                for (blasint jj = 0; jj < N_pad; jj += 32) {
+                    blasint current_N_block = stdmin(32, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 32) * panel_stride_B];
+                    for (blasint ii = 0; ii < M_pad; ii += 32) {
+                        int current_M_block = stdmin(32, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 32) * panel_stride_A];
+                        float *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+
+                        if (current_M_block == 32 && current_N_block == 32) {
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) float C_buffer[1024];
+                            if (beta_mode != 0) {
+                                for (int idx = 0; idx < 1024; ++idx) {
+                                    C_buffer[idx] = 0.0f;
+                                }
+                                for (blasint bc = 0; bc < current_N_block; ++bc) {
+                                    for (blasint br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 32] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_buffer, 32, beta_mode, &beta);
+                            for (blasint bc = 0; bc < current_N_block; ++bc) {
+                                for (blasint br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 32];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void sgemm_sme_TN(int M, int N, int K, float alpha, const float *A, int lda, const float *B, int ldb, float beta, float *C, int ldc)
+{
+    if (alpha == 0.0f || K == 0) {
+        if (beta != 1.0f) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == 0.0f) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC];
+
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 31) & ~31;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == 0.0f) ? 0 : (beta != 1.0f ? 1 : 2)) : 2;
+
+            int N_main = current_N & ~31;
+            for (int jj = 0; jj < N_main; jj += 32) {
+                float *__restrict B_out = &B_pack[(jj / 32) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    const float *__restrict B_in = &B[(j + jj) * (size_t)ldb + k + kk];
+                    for (int bc = 0; bc < 32; ++bc) {
+                        B_out[kk * 32 + bc] = B_in[bc * ldb];
+                    }
+                }
+            }
+            if (N_main < current_N) {
+                float *__restrict B_out = &B_pack[(N_main / 32) * panel_stride_B];
+                for (int bc = 0; bc < 32; ++bc) {
+                    if (N_main + bc < current_N) {
+                        const float *__restrict B_col = &B[(j + N_main + bc) * (size_t)ldb + k];
+                        for (int kk = 0; kk < current_K; ++kk) {
+                            B_out[kk * 32 + bc] = B_col[kk];
+                        }
+                    }
+                    else {
+                        for (int kk = 0; kk < current_K; ++kk) {
+                            B_out[kk * 32 + bc] = 0.0f;
+                        }
+                    }
+                }
+            }
+
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 31) & ~31;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+
+                int M_main = current_M & ~31;
+                for (int ii = 0; ii < M_main; ii += 32) {
+                    float *__restrict A_out = &A_pack[(ii / 32) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        const float *__restrict A_in = &A[(i + ii) * (size_t)lda + k + kk];
+                        for (int br = 0; br < 32; ++br) {
+                            A_out[kk * 32 + br] = A_in[br * lda] * alpha;
+                        }
+                    }
+                }
+                if (M_main < current_M) {
+                    float *__restrict A_out = &A_pack[(M_main / 32) * panel_stride_A];
+                    for (int br = 0; br < 32; ++br) {
+                        if (M_main + br < current_M) {
+                            const float *__restrict A_row = &A[(i + M_main + br) * (size_t)lda + k];
+                            for (int kk = 0; kk < current_K; ++kk) {
+                                A_out[kk * 32 + br] = A_row[kk] * alpha;
+                            }
+                        }
+                        else {
+                            for (int kk = 0; kk < current_K; ++kk) {
+                                A_out[kk * 32 + br] = 0.0f;
+                            }
+                        }
+                    }
+                }
+
+                for (int jj = 0; jj < N_pad; jj += 32) {
+                    int current_N_block = stdmin(32, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 32) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 32) {
+                        int current_M_block = stdmin(32, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 32) * panel_stride_A];
+                        float *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+
+                        if (current_M_block == 32 && current_N_block == 32) {
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) float C_buffer[1024];
+                            if (beta_mode != 0) {
+                                for (int idx = 0; idx < 1024; ++idx) {
+                                    C_buffer[idx] = 0.0f;
+                                }
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 32] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_buffer, 32, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 32];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void sgemm_sme_NT(int M, int N, int K, float alpha, const float *A, int lda, const float *B, int ldb, float beta, float *C, int ldc)
+{
+    if (alpha == 0.0f || K == 0) {
+        if (beta != 1.0f) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == 0.0f) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC];
+
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 31) & ~31;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == 0.0f) ? 0 : (beta != 1.0f ? 1 : 2)) : 2;
+
+            int N_main = current_N & ~31;
+            for (int jj = 0; jj < N_main; jj += 32) {
+                float *__restrict B_out = &B_pack[(jj / 32) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    const float *__restrict B_col = &B[(k + kk) * (size_t)ldb + j + jj];
+                    float *__restrict B_out_ptr = &B_out[kk * 32];
+                    for (int v = 0; v < 8; ++v) {
+                        vst1q_f32(&B_out_ptr[v * 4], vld1q_f32(&B_col[v * 4]));
+                    }
+                }
+            }
+            if (N_main < current_N) {
+                float *__restrict B_out = &B_pack[(N_main / 32) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    const float *__restrict B_row = &B[(k + kk) * (size_t)ldb + j + N_main];
+                    for (int bc = 0; bc < current_N - N_main; ++bc) {
+                        B_out[kk * 32 + bc] = B_row[bc];
+                    }
+                    for (int bc = current_N - N_main; bc < 32; ++bc) {
+                        B_out[kk * 32 + bc] = 0.0f;
+                    }
+                }
+            }
+
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 31) & ~31;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+                float32x4_t valpha = vdupq_n_f32(alpha);
+
+                int M_main = current_M & ~31;
+                for (int ii = 0; ii < M_main; ii += 32) {
+                    float *__restrict A_out = &A_pack[(ii / 32) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        const float *__restrict A_col = &A[(k + kk) * (size_t)lda + i + ii];
+                        float *__restrict A_out_ptr = &A_out[kk * 32];
+                        for (int v = 0; v < 8; ++v) {
+                            vst1q_f32(&A_out_ptr[v * 4], vmulq_f32(vld1q_f32(&A_col[v * 4]), valpha));
+                        }
+                    }
+                }
+                if (M_main < current_M) {
+                    float *__restrict A_out = &A_pack[(M_main / 32) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        const float *__restrict A_col = &A[(k + kk) * (size_t)lda + i + M_main];
+                        for (int br = 0; br < current_M - M_main; ++br) {
+                            A_out[kk * 32 + br] = A_col[br] * alpha;
+                        }
+                        for (int br = current_M - M_main; br < 32; ++br) {
+                            A_out[kk * 32 + br] = 0.0f;
+                        }
+                    }
+                }
+
+                for (int jj = 0; jj < N_pad; jj += 32) {
+                    int current_N_block = stdmin(32, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 32) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 32) {
+                        int current_M_block = stdmin(32, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 32) * panel_stride_A];
+                        float *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+
+                        if (current_M_block == 32 && current_N_block == 32) {
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) float C_buffer[1024];
+                            if (beta_mode != 0) {
+                                for (int idx = 0; idx < 1024; ++idx) {
+                                    C_buffer[idx] = 0.0f;
+                                }
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 32] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_buffer, 32, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 32];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void sgemm_sme_TT(int M, int N, int K, float alpha, const float *A, int lda, const float *B, int ldb, float beta, float *C, int ldc)
+{
+    if (alpha == 0.0f || K == 0) {
+        if (beta != 1.0f) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == 0.0f) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+
+    alignas(256) /*thread_local*/ static float A_pack[MC * KC];
+    alignas(256) /*thread_local*/ static float B_pack[KC * NC];
+
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 31) & ~31;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 32 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == 0.0f) ? 0 : (beta != 1.0f ? 1 : 2)) : 2;
+
+            int N_main = current_N & ~31;
+            for (int jj = 0; jj < N_main; jj += 32) {
+                float *__restrict B_out = &B_pack[(jj / 32) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    const float *__restrict B_col = &B[(k + kk) * (size_t)ldb + j + jj];
+                    float *__restrict B_out_ptr = &B_out[kk * 32];
+                    for (int v = 0; v < 8; ++v) {
+                        vst1q_f32(&B_out_ptr[v * 4], vld1q_f32(&B_col[v * 4]));
+                    }
+                }
+            }
+            if (N_main < current_N) {
+                float *__restrict B_out = &B_pack[(N_main / 32) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    const float *__restrict B_row = &B[(k + kk) * (size_t)ldb + j + N_main];
+                    for (int bc = 0; bc < current_N - N_main; ++bc) {
+                        B_out[kk * 32 + bc] = B_row[bc];
+                    }
+                    for (int bc = current_N - N_main; bc < 32; ++bc) {
+                        B_out[kk * 32 + bc] = 0.0f;
+                    }
+                }
+            }
+
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 31) & ~31;
+                size_t panel_stride_A = 32 * (size_t)current_K;
+
+                int M_main = current_M & ~31;
+                for (int ii = 0; ii < M_main; ii += 32) {
+                    float *__restrict A_out = &A_pack[(ii / 32) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        const float *__restrict A_in = &A[(i + ii) * (size_t)lda + k + kk];
+                        for (int br = 0; br < 32; ++br) {
+                            A_out[kk * 32 + br] = A_in[br * lda] * alpha;
+                        }
+                    }
+                }
+                if (M_main < current_M) {
+                    float *__restrict A_out = &A_pack[(M_main / 32) * panel_stride_A];
+                    for (int br = 0; br < 32; ++br) {
+                        if (M_main + br < current_M) {
+                            const float *__restrict A_row = &A[(i + M_main + br) * (size_t)lda + k];
+                            for (int kk = 0; kk < current_K; ++kk) {
+                                A_out[kk * 32 + br] = A_row[kk] * alpha;
+                            }
+                        }
+                        else {
+                            for (int kk = 0; kk < current_K; ++kk) {
+                                A_out[kk * 32 + br] = 0.0f;
+                            }
+                        }
+                    }
+                }
+
+                for (int jj = 0; jj < N_pad; jj += 32) {
+                    int current_N_block = stdmin(32, current_N - jj);
+                    float *B_ptr = &B_pack[(jj / 32) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 32) {
+                        int current_M_block = stdmin(32, current_M - ii);
+                        float *A_ptr = &A_pack[(ii / 32) * panel_stride_A];
+                        float *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+
+                        if (current_M_block == 32 && current_N_block == 32) {
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) float C_buffer[1024];
+                            if (beta_mode != 0) {
+                                for (int idx = 0; idx < 1024; ++idx) {
+                                    C_buffer[idx] = 0.0f;
+                                }
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 32] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            sgemm_sme_compute_32x32_tile(current_K, A_ptr, B_ptr, C_buffer, 32, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 32];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+void CNAME(char *transa,  char *transb, BLASLONG m,  BLASLONG n,  BLASLONG k,  float *alpha,  float *a,  BLASLONG lda,  float *b,  BLASLONG ldb,  float *beta, float *c,  BLASLONG ldc)
+{
+    bool trans_a = (*transa == 'T' || *transa == 't' || *transa == 'C' || *transa == 'c');
+    bool trans_b = (*transb == 'T' || *transb == 't' || *transb == 'C' || *transb == 'c');
+    if (!trans_a && !trans_b) {
+        sgemm_sme_NN(m, n, k, *alpha, a, lda, b, ldb, *beta, c, ldc);
+    }
+    else if (trans_a && !trans_b) {
+        sgemm_sme_TN(m, n, k, *alpha, a, lda, b, ldb, *beta, c, ldc);
+    }
+    else if (!trans_a && trans_b) {
+        sgemm_sme_NT(m, n, k, *alpha, a, lda, b, ldb, *beta, c, ldc);
+    }
+    else {
+        sgemm_sme_TT(m, n, k, *alpha, a, lda, b, ldb, *beta, c, ldc);
+    }
+}
@@ -0,0 +1,540 @@
+#include <arm_neon.h>
+#include <stddef.h>
+#include <stdalign.h>
+#include <stdbool.h>
+#include "common.h"
+#ifndef stdmin
+#define stdmin(a,b)   (a>b? b:a)
+#endif
+
+typedef double _Complex zdouble;
+
+static zdouble CDMUL(zdouble a, zdouble b,bool conja, bool conjb) {
+double ra=creal(a);
+double rb=creal(b);
+double ia=conja ? -cimag(a) : cimag(a);
+double ib=conjb ? -cimag(b) : cimag(b);
+double r1=ra*rb;
+double r2=ia*ib;
+double r=r1-r2;
+double i=(ra+ia)*(rb+ib)-r1-r2;
+zdouble res={r,i};
+return res;
+}
+
+
+static zdouble cdzero={0.,0.};
+static zdouble cdone={1.,0.};
+
+
+#define MC 128
+#define KC 1024
+#define NC 512
+
+static inline void zgemm_sme_compute_8x8_tile(int current_K, const double *A_ptr, const double *B_ptr, zdouble *C_ptr, size_t ldc, int beta_mode, const zdouble *beta_ptr)
+{
+    size_t ldc_bytes = ldc * sizeof(zdouble);
+
+    asm volatile("smstart\n\t"
+		 "ptrue p0.d\n\t"
+                 "zero {za}\n\t"
+                 "cmp %w[beta_mode], #0\n\t"
+                 "b.eq 19f\n\t"
+
+                 // FIX: Load Beta Real and Imaginary components and replicate them
+                 "cmp %w[beta_mode], #1\n\t"
+                 "b.ne 10f\n\t"
+                 "ld1rd z30.d, p0/z, [%[beta_ptr]]\n\t" // Load beta.real()
+                 "add x15, %[beta_ptr], #8\n\t"
+                 "ld1rd z31.d, p0/z, [x15]\n\t" // Load beta.imag()
+                 "10:\n\t"
+
+                 "mov w12, #0\n\t"
+                 "mov x13, %[c]\n\t"
+                 "100:\n\t"
+                 "ld1d z0.d, p0/z, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "ld1d z1.d, p0/z, [x14]\n\t"
+                 "uzp1 z2.d, z0.d, z1.d\n\t" // z2 = C_re
+                 "uzp2 z3.d, z0.d, z1.d\n\t" // z3 = C_im
+
+                 "cmp %w[beta_mode], #1\n\t"
+                 "b.ne 101f\n\t"
+
+                 // FIX: Fully implemented Complex Beta Multiplication using movprfx
+                 "movprfx z4, z2\n\t"
+                 "fmul z4.d, p0/m, z4.d, z30.d\n\t" // z4 = Cre * Bre
+                 "movprfx z5, z3\n\t"
+                 "fmul z5.d, p0/m, z5.d, z31.d\n\t" // z5 = Cim * Bim
+                 "movprfx z6, z4\n\t"
+                 "fsub z6.d, p0/m, z6.d, z5.d\n\t" // z6 = Cre' (Cre*Bre - Cim*Bim)
+
+                 "movprfx z8, z2\n\t"
+                 "fmul z8.d, p0/m, z8.d, z31.d\n\t" // z8 = Cre * Bim
+                 "movprfx z9, z3\n\t"
+                 "fmul z9.d, p0/m, z9.d, z30.d\n\t" // z9 = Cim * Bre
+                 "movprfx z7, z8\n\t"
+                 "fadd z7.d, p0/m, z7.d, z9.d\n\t" // z7 = Cim' (Cre*Bim + Cim*Bre)
+
+                 "mova za0v.d[w12, 0], p0/m, z6.d\n\t"
+                 "mova za2v.d[w12, 0], p0/m, z7.d\n\t"
+                 "b 102f\n\t"
+
+                 "101:\n\t" // Fallback: beta == 1.0
+                 "mova za0v.d[w12, 0], p0/m, z2.d\n\t"
+                 "mova za2v.d[w12, 0], p0/m, z3.d\n\t"
+
+                 "102:\n\t"
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w12, w12, #1\n\t"
+                 "cmp w12, #8\n\t"
+                 "b.ne 100b\n\t"
+
+                 "19:\n\t"
+                 "mov w10, %w[k]\n\t"
+                 "cbz w10, 3f\n\t"
+                 "11:\n\t"
+                 "ld1d z0.d, p0/z, [%[a], #0, mul vl]\n\t"
+                 "ld1d z1.d, p0/z, [%[a], #1, mul vl]\n\t"
+                 "ld1d z2.d, p0/z, [%[b], #0, mul vl]\n\t"
+                 "ld1d z3.d, p0/z, [%[b], #1, mul vl]\n\t"
+                 "fmopa za0.d, p0/m, p0/m, z0.d, z2.d\n\t"
+                 "fmopa za1.d, p0/m, p0/m, z1.d, z3.d\n\t"
+                 "fmopa za2.d, p0/m, p0/m, z0.d, z3.d\n\t"
+                 "fmopa za3.d, p0/m, p0/m, z1.d, z2.d\n\t"
+                 "add %[a], %[a], #128\n\t"
+                 "add %[b], %[b], #128\n\t"
+                 "subs w10, w10, #1\n\t"
+                 "b.ne 11b\n\t"
+
+                 "3:\n\t"
+                 "mov w12, #0\n\t"
+                 "mov x13, %[c]\n\t"
+                 "200:\n\t"
+                 "mova z0.d, p0/m, za0v.d[w12, 0]\n\t"
+                 "mova z1.d, p0/m, za1v.d[w12, 0]\n\t"
+                 "mova z2.d, p0/m, za2v.d[w12, 0]\n\t"
+                 "mova z3.d, p0/m, za3v.d[w12, 0]\n\t"
+
+                 // FIX: Non-destructive SVE arithmetic
+                 "movprfx z4, z0\n\t"
+                 "fsub z4.d, p0/m, z4.d, z1.d\n\t"
+                 "movprfx z5, z2\n\t"
+                 "fadd z5.d, p0/m, z5.d, z3.d\n\t"
+
+                 "zip1 z6.d, z4.d, z5.d\n\t"
+                 "zip2 z7.d, z4.d, z5.d\n\t"
+                 "st1d z6.d, p0, [x13]\n\t"
+                 "add x14, x13, #64\n\t"
+                 "st1d z7.d, p0, [x14]\n\t"
+                 "add x13, x13, %[ldc_bytes]\n\t"
+                 "add w12, w12, #1\n\t"
+                 "cmp w12, #8\n\t"
+                 "b.ne 200b\n\t"
+ 		 "smstop\n\t"
+                 "msr fpsr, xzr\n\t"
+    : [a] "+&r"(A_ptr), [b] "+&r"(B_ptr)
+    : [k] "r"(current_K), [c] "r"(C_ptr), [ldc_bytes] "r"(ldc_bytes), [beta_mode] "r"(beta_mode), [beta_ptr] "r"(beta_ptr)
+    : "x10", "x12", "x13", "x14", "x15", "memory", "cc",
+      "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7",
+      "v8", "v9", "v10", "v11", "v12", "v13", "v14", "v15",
+      "v16", "v17", "v18", "v19", "v20", "v21", "v22", "v23",
+      "v24", "v25", "v26", "v27", "v28", "v29", "v30", "v31",
+      "z0", "z1", "z2", "z3", "z4", "z5", "z6", "z7",
+      "z8", "z9", "z10", "z11", "z12", "z13", "z14", "z15",
+      "z16", "z17", "z18", "z19", "z20", "z21", "z22", "z23",
+      "z24", "z25", "z26", "z27", "z28", "z29", "z30", "z31", "za",
+      "p0", "p1", "p2", "p3", "p4", "p5", "p6", "p7", "p8",
+      "p9", "p10", "p11", "p12", "p13", "p14", "p15");
+
+}
+
+static void zgemm_sme_NN(int M, int N, int K, const zdouble alpha, const zdouble *A, int lda, const zdouble *b, int ldb, const zdouble beta, zdouble *C, int ldc, bool conja, bool conjb)
+{
+
+    if (alpha == cdzero || K == 0) {
+        if (beta != cdone) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == cdzero) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+
+    alignas(256) /*thread_local*/ static double A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static double B_pack[KC * NC * 2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 7) & ~7;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 16 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == cdzero) ? 0 : (beta != cdone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 8) {
+                double *__restrict B_out = &B_pack[(jj / 8) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    for (int bc = 0; bc < 8; ++bc) {
+                        zdouble val = (jj + bc < current_N) ? b[(j + jj + bc) * ldb + k + kk] : cdzero;
+                        B_out[kk * 16 + bc] = creal(val);
+                        B_out[kk * 16 + 8 + bc] = conjb ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 7) & ~7;
+                size_t panel_stride_A = 16 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 8) {
+                    double *__restrict A_out = &A_pack[(ii / 8) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        for (int br = 0; br < 8; ++br) {
+                            zdouble val = (ii + br < current_M) ? CDMUL(A[(k + kk) * lda + i + ii + br] ,alpha,conja,0) : cdzero;
+                            A_out[kk * 16 + br] = creal(val);
+                            A_out[kk * 16 + 8 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 8) {
+                    int current_N_block = stdmin(8, current_N - jj);
+                    double *B_ptr = &B_pack[(jj / 8) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 8) {
+                        int current_M_block = stdmin(8, current_M - ii);
+                        double *A_ptr = &A_pack[(ii / 8) * panel_stride_A];
+                        zdouble *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 8 && current_N_block == 8) {
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) zdouble C_buffer[64];
+                            
+                                for (int idx = 0; idx < 64; ++idx) {
+                                    C_buffer[idx] = cdzero;
+                                }
+                            if (beta_mode != 0) {
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 8] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_buffer, 8, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 8];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void zgemm_sme_TN(int M, int N, int K, const zdouble alpha, const zdouble *A, int lda, const zdouble *b, int ldb, const zdouble beta, zdouble *C, int ldc, bool conja, bool conjb)
+{
+    if (alpha == cdzero || K == 0) {
+        if (beta != cdone) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == cdzero) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+
+    alignas(256) /*thread_local*/ static double A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static double B_pack[KC * NC * 2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 7) & ~7;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 16 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == cdzero) ? 0 : (beta != cdone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 8) {
+                double *__restrict B_out = &B_pack[(jj / 8) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    for (int bc = 0; bc < 8; ++bc) {
+                        zdouble val = (jj + bc < current_N) ? b[(j + jj + bc) * ldb + k + kk] : cdzero;
+                        B_out[kk * 16 + bc] = creal(val);
+                        B_out[kk * 16 + 8 + bc] = conjb ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 7) & ~7;
+                size_t panel_stride_A = 16 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 8) {
+                    double *__restrict A_out = &A_pack[(ii / 8) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        for (int br = 0; br < 8; ++br) {
+                            zdouble val = (ii + br < current_M) ? CDMUL(A[(i + ii + br) * lda + k + kk],alpha,conja,0) : cdzero;
+                            A_out[kk * 16 + br] = creal(val);
+                            A_out[kk * 16 + 8 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 8) {
+                    int current_N_block = stdmin(8, current_N - jj);
+                    double *B_ptr = &B_pack[(jj / 8) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 8) {
+                        int current_M_block = stdmin(8, current_M - ii);
+                        double *A_ptr = &A_pack[(ii / 8) * panel_stride_A];
+                        zdouble *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 8 && current_N_block == 8) {
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) zdouble C_buffer[64];
+                            
+                                for (int idx = 0; idx < 64; ++idx) {
+                                    C_buffer[idx] = cdzero;
+                                }
+                            if (beta_mode != 0) {
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 8] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_buffer, 8, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 8];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void zgemm_sme_NT(blasint M, blasint N, blasint K, const zdouble alpha, const zdouble *A, blasint lda, const zdouble *b, blasint ldb, const zdouble beta, zdouble *C, blasint ldc, bool conja, bool conjb)
+{
+#if 0
+    if (alpha == cdzero || K == 0) {
+        if (beta != cdone) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == cdzero) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+#endif
+    if (alpha == cdzero || K == 0) {
+        if (beta == cdzero) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = cdzero;
+                }
+            }
+        } else {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] *= beta;
+                }
+            }
+        }
+        return;
+    }
+
+
+    alignas(256) /*thread_local*/ static double A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static double B_pack[KC * NC * 2 *2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 7) & ~7;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 16 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == cdzero) ? 0 : (beta != cdone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 8) {
+                double *__restrict B_out = &B_pack[(jj / 8) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    for (int bc = 0; bc < 8; ++bc) {
+                        zdouble val = (jj + bc < current_N) ? b[(k + kk) * ldb + j + jj + bc] : cdzero;
+                        B_out[kk * 16 + bc] = creal(val);
+                        B_out[kk * 16 + 8 + bc] = conjb ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 7) & ~7;
+                size_t panel_stride_A = 16 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 8) {
+                    double *__restrict A_out = &A_pack[(ii / 8) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        for (int br = 0; br < 8; ++br) {
+                            zdouble val = (ii + br < current_M) ? CDMUL(A[(k + kk) * lda + i + ii + br] ,alpha,conja,0) : cdzero;
+                            A_out[kk * 16 + br] = creal(val);
+                            A_out[kk * 16 + 8 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 8) {
+                    int current_N_block = stdmin(8, current_N - jj);
+                    double *B_ptr = &B_pack[(jj / 8) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 8) {
+                        int current_M_block = stdmin(8, current_M - ii);
+                        double *A_ptr = &A_pack[(ii / 8) * panel_stride_A];
+                        zdouble *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 8 && current_N_block == 8) {
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) zdouble C_buffer[64];
+                            
+                                for (int idx = 0; idx < 64; ++idx) {
+                                    C_buffer[idx] = cdzero;
+                                }
+                            if (beta_mode != 0) {
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 8] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_buffer, 8, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 8];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+static void zgemm_sme_TT(int M, int N, int K, const zdouble alpha, const zdouble *A, int lda, const zdouble *b, int ldb, const zdouble beta, zdouble *C, int ldc, bool conja, bool conjb)
+{
+    if (alpha == cdzero || K == 0) {
+        if (beta != cdone) {
+            for (int j = 0; j < N; ++j) {
+                for (int i = 0; i < M; ++i) {
+                    C[i + j * ldc] = (beta == cdzero) ? 0.0f : C[i + j * ldc] * beta;
+                }
+            }
+        }
+        return;
+    }
+
+    alignas(256) /*thread_local*/ static double A_pack[MC * KC * 2];
+    alignas(256) /*thread_local*/ static double B_pack[KC * NC * 2];
+#pragma omp parallel for schedule(dynamic, 1)
+    for (int j = 0; j < N; j += NC) {
+        int current_N = stdmin(NC, N - j);
+        int N_pad = (current_N + 7) & ~7;
+        for (int k = 0; k < K; k += KC) {
+            int current_K = stdmin(KC, K - k);
+            size_t panel_stride_B = 16 * (size_t)current_K;
+            int beta_mode = (k == 0) ? ((beta == cdzero) ? 0 : (beta != cdone ? 1 : 2)) : 2;
+
+            for (int jj = 0; jj < current_N; jj += 8) {
+                double *__restrict B_out = &B_pack[(jj / 8) * panel_stride_B];
+                for (int kk = 0; kk < current_K; ++kk) {
+                    for (int bc = 0; bc < 8; ++bc) {
+                        zdouble val = (jj + bc < current_N) ? b[(k + kk) * ldb + j + jj + bc] : cdzero;
+                        B_out[kk * 16 + bc] = creal(val);
+                        B_out[kk * 16 + 8 + bc] = conjb ? -cimag(val) : cimag(val);
+                    }
+                }
+            }
+            for (int i = 0; i < M; i += MC) {
+                int current_M = stdmin(MC, M - i);
+                int M_pad = (current_M + 7) & ~7;
+                size_t panel_stride_A = 16 * (size_t)current_K;
+                for (int ii = 0; ii < current_M; ii += 8) {
+                    double *__restrict A_out = &A_pack[(ii / 8) * panel_stride_A];
+                    for (int kk = 0; kk < current_K; ++kk) {
+                        for (int br = 0; br < 8; ++br) {
+                            zdouble val = (ii + br < current_M) ? CDMUL(A[(i + ii + br) * lda + k + kk],alpha,conja,0) : cdzero;
+                            A_out[kk * 16 + br] = creal(val);
+                            A_out[kk * 16 + 8 + br] = cimag(val);
+                        }
+                    }
+                }
+                for (int jj = 0; jj < N_pad; jj += 8) {
+                    int current_N_block = stdmin(8, current_N - jj);
+                    double *B_ptr = &B_pack[(jj / 8) * panel_stride_B];
+                    for (int ii = 0; ii < M_pad; ii += 8) {
+                        int current_M_block = stdmin(8, current_M - ii);
+                        double *A_ptr = &A_pack[(ii / 8) * panel_stride_A];
+                        zdouble *C_ptr = &C[(i + ii) + (j + jj) * ldc];
+                        if (current_M_block == 8 && current_N_block == 8) {
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_ptr, ldc, beta_mode, &beta);
+                        }
+                        else {
+                            alignas(256) zdouble C_buffer[64];
+                            
+                                for (int idx = 0; idx < 64; ++idx) {
+                                    C_buffer[idx] = cdzero;
+                                }
+                            if (beta_mode != 0) {
+                                for (int bc = 0; bc < current_N_block; ++bc) {
+                                    for (int br = 0; br < current_M_block; ++br) {
+                                        C_buffer[br + bc * 8] = C_ptr[br + bc * ldc];
+                                    }
+                                }
+                            }
+                            zgemm_sme_compute_8x8_tile(current_K, A_ptr, B_ptr, C_buffer, 8, beta_mode, &beta);
+                            for (int bc = 0; bc < current_N_block; ++bc) {
+                                for (int br = 0; br < current_M_block; ++br) {
+                                    C_ptr[br + bc * ldc] = C_buffer[br + bc * 8];
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+void CNAME(const char *transa, const char *transb, const BLASLONG m, const BLASLONG n, const BLASLONG k, const double alpha_r, const double alpha_i, const double *a, const BLASLONG lda, const double *b, const BLASLONG ldb, const double beta_r, const double beta_i, double *c, const BLASLONG ldc)
+{
+zdouble alpha={alpha_r,alpha_i};
+zdouble beta={beta_r,beta_i};
+
+    bool trans_a = (*transa == 'T' || *transa == 't' || *transa == 'C' || *transa == 'c');
+    bool trans_b = (*transb == 'T' || *transb == 't' || *transb == 'C' || *transb == 'c');
+    if (!trans_a && !trans_b) {
+	bool conja = (*transa == 'R' || *transa == 'r');
+	bool conjb = (*transb == 'R' || *transb == 'r');
+        zgemm_sme_NN(m, n, k, alpha, (const zdouble*)a, lda, (const zdouble*)b, ldb, beta, (zdouble*)c, ldc, conja, conjb);
+    }
+    else if (trans_a && !trans_b) {
+	bool conja = (*transa == 'C' || *transa == 'c');
+	bool conjb = (*transb == 'R' || *transb == 'r');
+        zgemm_sme_TN(m, n, k, alpha, (const zdouble*)a, lda, (const zdouble*)b, ldb, beta, (zdouble*)c, ldc, conja, conjb);
+    }
+    else if (!trans_a && trans_b) {
+	bool conja = (*transa == 'R' || *transa == 'r');
+	bool conjb = (*transb == 'C' || *transb == 'c');
+        zgemm_sme_NT(m, n, k, alpha, (const zdouble*)a, lda, (const zdouble*)b, ldb, beta, (zdouble*)c, ldc, conja, conjb);
+    }
+    else {
+	bool conja = (*transa == 'C' || *transa == 'c');
+	bool conjb = (*transb == 'C' || *transb == 'c');
+        zgemm_sme_TT(m, n, k, alpha, (const zdouble*)a, lda, (const zdouble*)b, ldb, beta, (zdouble*)c, ldc, conja, conjb);
+    }
+}
+
@@ -0,0 +1,346 @@
+/*********************************************************************/
+/* Copyright 2009, 2010 The University of Texas at Austin.           */
+/* All rights reserved.                                              */
+/*                                                                   */
+/* Redistribution and use in source and binary forms, with or        */
+/* without modification, are permitted provided that the following   */
+/* conditions are met:                                               */
+/*                                                                   */
+/*   1. Redistributions of source code must retain the above         */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer.                                                  */
+/*                                                                   */
+/*   2. Redistributions in binary form must reproduce the above      */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer in the documentation and/or other materials       */
+/*      provided with the distribution.                              */
+/*                                                                   */
+/*    THIS  SOFTWARE IS PROVIDED  BY THE  UNIVERSITY OF  TEXAS AT    */
+/*    AUSTIN  ``AS IS''  AND ANY  EXPRESS OR  IMPLIED WARRANTIES,    */
+/*    INCLUDING, BUT  NOT LIMITED  TO, THE IMPLIED  WARRANTIES OF    */
+/*    MERCHANTABILITY  AND FITNESS FOR  A PARTICULAR  PURPOSE ARE    */
+/*    DISCLAIMED.  IN  NO EVENT SHALL THE UNIVERSITY  OF TEXAS AT    */
+/*    AUSTIN OR CONTRIBUTORS BE  LIABLE FOR ANY DIRECT, INDIRECT,    */
+/*    INCIDENTAL,  SPECIAL, EXEMPLARY,  OR  CONSEQUENTIAL DAMAGES    */
+/*    (INCLUDING, BUT  NOT LIMITED TO,  PROCUREMENT OF SUBSTITUTE    */
+/*    GOODS  OR  SERVICES; LOSS  OF  USE,  DATA,  OR PROFITS;  OR    */
+/*    BUSINESS INTERRUPTION) HOWEVER CAUSED  AND ON ANY THEORY OF    */
+/*    LIABILITY, WHETHER  IN CONTRACT, STRICT  LIABILITY, OR TORT    */
+/*    (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY WAY OUT    */
+/*    OF  THE  USE OF  THIS  SOFTWARE,  EVEN  IF ADVISED  OF  THE    */
+/*    POSSIBILITY OF SUCH DAMAGE.                                    */
+/*                                                                   */
+/* The views and conclusions contained in the software and           */
+/* documentation are those of the authors and should not be          */
+/* interpreted as representing official policies, either expressed   */
+/* or implied, of The University of Texas at Austin.                 */
+/*********************************************************************/
+
+#include "common.h"
+
+/* UNROLL_M=6-aware TRSM kernel.
+ * Stock generic/trsm_kernel_*.c assumes power-of-2 UNROLL_M via
+ * GEMM_UNROLL_M_SHIFT (incorrectly set to 2 for M=6) and bit masks
+ * (m & (M-1)). That overruns packed buffers and corrupts the heap —
+ * reproduces as free()/double-free in HPL panel factor (dtrsm).
+ * Pattern mirrors loongarch64/trsm_kernel_*_UNROLLN6.c for N=6.
+ */
+
+static FLOAT dm1 = -1.;
+
+#ifdef CONJ
+#define GEMM_KERNEL   GEMM_KERNEL_L
+#else
+#define GEMM_KERNEL   GEMM_KERNEL_N
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 1
+#define GEMM_UNROLL_M_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 2
+#define GEMM_UNROLL_M_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 4
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 6
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 8
+#define GEMM_UNROLL_M_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 16
+#define GEMM_UNROLL_M_SHIFT 4
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 1
+#define GEMM_UNROLL_N_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 2
+#define GEMM_UNROLL_N_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 4
+#define GEMM_UNROLL_N_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 8
+#define GEMM_UNROLL_N_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 16
+#define GEMM_UNROLL_N_SHIFT 4
+#endif
+
+#ifndef COMPLEX
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa,  bb;
+
+  int i, j, k;
+
+  a += (m - 1) * m;
+  b += (m - 1) * n;
+
+  for (i = m - 1; i >= 0; i--) {
+
+    aa = *(a + i);
+
+    for (j = 0; j < n; j ++) {
+      bb = *(c + i + j * ldc);
+      bb *= aa;
+      *b             = bb;
+      *(c + i + j * ldc) = bb;
+      b ++;
+
+      for (k = 0; k < i; k ++){
+	*(c + k + j * ldc) -= bb * *(a + k);
+      }
+
+    }
+    a -= m;
+    b -= 2 * n;
+  }
+
+}
+
+#else
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa1, aa2;
+  FLOAT bb1, bb2;
+  FLOAT cc1, cc2;
+
+  int i, j, k;
+
+  ldc *= 2;
+  a += (m - 1) * m * 2;
+  b += (m - 1) * n * 2;
+
+  for (i = m - 1; i >= 0; i--) {
+
+    aa1 = *(a + i * 2 + 0);
+    aa2 = *(a + i * 2 + 1);
+
+    for (j = 0; j < n; j ++) {
+      bb1 = *(c + i * 2 + 0 + j * ldc);
+      bb2 = *(c + i * 2 + 1 + j * ldc);
+
+#ifndef CONJ
+      cc1 = aa1 * bb1 - aa2 * bb2;
+      cc2 = aa1 * bb2 + aa2 * bb1;
+#else
+      cc1 = aa1 * bb1 + aa2 * bb2;
+      cc2 = aa1 * bb2 - aa2 * bb1;
+#endif
+
+
+      *(b + 0) = cc1;
+      *(b + 1) = cc2;
+      *(c + i * 2 + 0 + j * ldc) = cc1;
+      *(c + i * 2 + 1 + j * ldc) = cc2;
+      b += 2;
+
+      for (k = 0; k < i; k ++){
+#ifndef CONJ
+	*(c + k * 2 + 0 + j * ldc) -= cc1 * *(a + k * 2 + 0) - cc2 * *(a + k * 2 + 1);
+	*(c + k * 2 + 1 + j * ldc) -= cc1 * *(a + k * 2 + 1) + cc2 * *(a + k * 2 + 0);
+#else
+	*(c + k * 2 + 0 + j * ldc) -=   cc1 * *(a + k * 2 + 0) + cc2 * *(a + k * 2 + 1);
+	*(c + k * 2 + 1 + j * ldc) -= - cc1 * *(a + k * 2 + 1) + cc2 * *(a + k * 2 + 0);
+#endif
+      }
+
+    }
+    a -= m * 2;
+    b -= 4 * n;
+  }
+
+}
+
+#endif
+
+
+int CNAME(BLASLONG m, BLASLONG n, BLASLONG k,  FLOAT dummy1,
+#ifdef COMPLEX
+	   FLOAT dummy2,
+#endif
+	   FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG offset){
+
+  BLASLONG i, j;
+  FLOAT *aa, *cc;
+  BLASLONG  kk;
+  BLASLONG  mmodM = m - (m / GEMM_UNROLL_M) * GEMM_UNROLL_M;
+
+#if 0
+  fprintf(stderr, "TRSM KERNEL LN : m = %3ld  n = %3ld  k = %3ld offset = %3ld\n",
+	  m, n, k, offset);
+#endif
+
+  j = (n >> GEMM_UNROLL_N_SHIFT);
+
+  while (j > 0) {
+
+    kk = m + offset;
+
+    if (mmodM) {
+      for (i = 1; i < GEMM_UNROLL_M; i *= 2){
+	if (mmodM & i) {
+	  aa = a + ((m & ~(i - 1)) - i) * k * COMPSIZE;
+	  cc = c + ((m & ~(i - 1)) - i)     * COMPSIZE;
+
+	  if (k - kk > 0) {
+	    GEMM_KERNEL(i, GEMM_UNROLL_N, k - kk, dm1,
+#ifdef COMPLEX
+			ZERO,
+#endif
+			aa + i             * kk * COMPSIZE,
+			b  + GEMM_UNROLL_N * kk * COMPSIZE,
+			cc,
+			ldc);
+	  }
+
+	  solve(i, GEMM_UNROLL_N,
+		aa + (kk - i) * i             * COMPSIZE,
+		b  + (kk - i) * GEMM_UNROLL_N * COMPSIZE,
+		cc, ldc);
+
+	  kk -= i;
+	}
+      }
+    }
+
+    i = (m / GEMM_UNROLL_M);
+    if (i > 0) {
+      aa = a + ((m / GEMM_UNROLL_M) * GEMM_UNROLL_M - GEMM_UNROLL_M) * k * COMPSIZE;
+      cc = c + ((m / GEMM_UNROLL_M) * GEMM_UNROLL_M - GEMM_UNROLL_M)     * COMPSIZE;
+
+      do {
+	if (k - kk > 0) {
+	  GEMM_KERNEL(GEMM_UNROLL_M, GEMM_UNROLL_N, k - kk, dm1,
+#ifdef COMPLEX
+		      ZERO,
+#endif
+		      aa + GEMM_UNROLL_M * kk * COMPSIZE,
+		      b +  GEMM_UNROLL_N * kk * COMPSIZE,
+		      cc,
+		      ldc);
+	}
+
+	solve(GEMM_UNROLL_M, GEMM_UNROLL_N,
+	      aa + (kk - GEMM_UNROLL_M) * GEMM_UNROLL_M * COMPSIZE,
+	      b  + (kk - GEMM_UNROLL_M) * GEMM_UNROLL_N * COMPSIZE,
+	      cc, ldc);
+
+	aa -= GEMM_UNROLL_M * k * COMPSIZE;
+	cc -= GEMM_UNROLL_M     * COMPSIZE;
+	kk -= GEMM_UNROLL_M;
+	i --;
+      } while (i > 0);
+    }
+
+    b += GEMM_UNROLL_N * k * COMPSIZE;
+    c += GEMM_UNROLL_N * ldc * COMPSIZE;
+    j --;
+  }
+
+  if (n & (GEMM_UNROLL_N - 1)) {
+
+    j = (GEMM_UNROLL_N >> 1);
+    while (j > 0) {
+      if (n & j) {
+
+	kk = m + offset;
+
+	if (mmodM) {
+	  for (i = 1; i < GEMM_UNROLL_M; i *= 2){
+	    if (mmodM & i) {
+	      aa = a + ((m & ~(i - 1)) - i) * k * COMPSIZE;
+	      cc = c + ((m & ~(i - 1)) - i)     * COMPSIZE;
+
+	      if (k - kk > 0) {
+		GEMM_KERNEL(i, j, k - kk, dm1,
+#ifdef COMPLEX
+			    ZERO,
+#endif
+			    aa + i * kk * COMPSIZE,
+			    b  + j * kk * COMPSIZE,
+			    cc, ldc);
+	      }
+
+	      solve(i, j,
+		    aa + (kk - i) * i * COMPSIZE,
+		    b  + (kk - i) * j * COMPSIZE,
+		    cc, ldc);
+
+	      kk -= i;
+	    }
+	  }
+	}
+
+	i = (m / GEMM_UNROLL_M);
+	if (i > 0) {
+	  aa = a + ((m / GEMM_UNROLL_M) * GEMM_UNROLL_M - GEMM_UNROLL_M) * k * COMPSIZE;
+	  cc = c + ((m / GEMM_UNROLL_M) * GEMM_UNROLL_M - GEMM_UNROLL_M)     * COMPSIZE;
+
+	  do {
+	    if (k - kk > 0) {
+	      GEMM_KERNEL(GEMM_UNROLL_M, j, k - kk, dm1,
+#ifdef COMPLEX
+			  ZERO,
+#endif
+			  aa + GEMM_UNROLL_M * kk * COMPSIZE,
+			  b +  j             * kk * COMPSIZE,
+			  cc,
+			  ldc);
+	    }
+
+	    solve(GEMM_UNROLL_M, j,
+		  aa + (kk - GEMM_UNROLL_M) * GEMM_UNROLL_M * COMPSIZE,
+		  b  + (kk - GEMM_UNROLL_M) * j             * COMPSIZE,
+		  cc, ldc);
+
+	    aa -= GEMM_UNROLL_M * k * COMPSIZE;
+	    cc -= GEMM_UNROLL_M     * COMPSIZE;
+	    kk -= GEMM_UNROLL_M;
+	    i --;
+	  } while (i > 0);
+	}
+
+	b += j * k   * COMPSIZE;
+	c += j * ldc * COMPSIZE;
+      }
+      j >>= 1;
+    }
+  }
+
+  return 0;
+}
@@ -0,0 +1,330 @@
+/*********************************************************************/
+/* Copyright 2009, 2010 The University of Texas at Austin.           */
+/* All rights reserved.                                              */
+/*                                                                   */
+/* Redistribution and use in source and binary forms, with or        */
+/* without modification, are permitted provided that the following   */
+/* conditions are met:                                               */
+/*                                                                   */
+/*   1. Redistributions of source code must retain the above         */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer.                                                  */
+/*                                                                   */
+/*   2. Redistributions in binary form must reproduce the above      */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer in the documentation and/or other materials       */
+/*      provided with the distribution.                              */
+/*                                                                   */
+/*    THIS  SOFTWARE IS PROVIDED  BY THE  UNIVERSITY OF  TEXAS AT    */
+/*    AUSTIN  ``AS IS''  AND ANY  EXPRESS OR  IMPLIED WARRANTIES,    */
+/*    INCLUDING, BUT  NOT LIMITED  TO, THE IMPLIED  WARRANTIES OF    */
+/*    MERCHANTABILITY  AND FITNESS FOR  A PARTICULAR  PURPOSE ARE    */
+/*    DISCLAIMED.  IN  NO EVENT SHALL THE UNIVERSITY  OF TEXAS AT    */
+/*    AUSTIN OR CONTRIBUTORS BE  LIABLE FOR ANY DIRECT, INDIRECT,    */
+/*    INCIDENTAL,  SPECIAL, EXEMPLARY,  OR  CONSEQUENTIAL DAMAGES    */
+/*    (INCLUDING, BUT  NOT LIMITED TO,  PROCUREMENT OF SUBSTITUTE    */
+/*    GOODS  OR  SERVICES; LOSS  OF  USE,  DATA,  OR PROFITS;  OR    */
+/*    BUSINESS INTERRUPTION) HOWEVER CAUSED  AND ON ANY THEORY OF    */
+/*    LIABILITY, WHETHER  IN CONTRACT, STRICT  LIABILITY, OR TORT    */
+/*    (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY WAY OUT    */
+/*    OF  THE  USE OF  THIS  SOFTWARE,  EVEN  IF ADVISED  OF  THE    */
+/*    POSSIBILITY OF SUCH DAMAGE.                                    */
+/*                                                                   */
+/* The views and conclusions contained in the software and           */
+/* documentation are those of the authors and should not be          */
+/* interpreted as representing official policies, either expressed   */
+/* or implied, of The University of Texas at Austin.                 */
+/*********************************************************************/
+
+#include "common.h"
+
+/* UNROLL_M=6-aware TRSM kernel.
+ * Stock generic/trsm_kernel_*.c assumes power-of-2 UNROLL_M via
+ * GEMM_UNROLL_M_SHIFT (incorrectly set to 2 for M=6) and bit masks
+ * (m & (M-1)). That overruns packed buffers and corrupts the heap —
+ * reproduces as free()/double-free in HPL panel factor (dtrsm).
+ * Pattern mirrors loongarch64/trsm_kernel_*_UNROLLN6.c for N=6.
+ */
+
+static FLOAT dm1 = -1.;
+
+#ifdef CONJ
+#define GEMM_KERNEL   GEMM_KERNEL_L
+#else
+#define GEMM_KERNEL   GEMM_KERNEL_N
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 1
+#define GEMM_UNROLL_M_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 2
+#define GEMM_UNROLL_M_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 4
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 6
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 8
+#define GEMM_UNROLL_M_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 16
+#define GEMM_UNROLL_M_SHIFT 4
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 1
+#define GEMM_UNROLL_N_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 2
+#define GEMM_UNROLL_N_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 4
+#define GEMM_UNROLL_N_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 8
+#define GEMM_UNROLL_N_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 16
+#define GEMM_UNROLL_N_SHIFT 4
+#endif
+
+#ifndef COMPLEX
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa, bb;
+
+  int i, j, k;
+
+  for (i = 0; i < m; i++) {
+
+    aa = *(a + i);
+
+    for (j = 0; j < n; j ++) {
+      bb = *(c + i + j * ldc);
+      bb *= aa;
+      *b             = bb;
+      *(c + i + j * ldc) = bb;
+      b ++;
+
+      for (k = i + 1; k < m; k ++){
+	*(c + k + j * ldc) -= bb * *(a + k);
+      }
+
+    }
+    a += m;
+  }
+}
+
+#else
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa1, aa2;
+  FLOAT bb1, bb2;
+  FLOAT cc1, cc2;
+
+  int i, j, k;
+
+  ldc *= 2;
+
+  for (i = 0; i < m; i++) {
+
+    aa1 = *(a + i * 2 + 0);
+    aa2 = *(a + i * 2 + 1);
+
+    for (j = 0; j < n; j ++) {
+      bb1 = *(c + i * 2 + 0 + j * ldc);
+      bb2 = *(c + i * 2 + 1 + j * ldc);
+
+#ifndef CONJ
+      cc1 = aa1 * bb1 - aa2 * bb2;
+      cc2 = aa1 * bb2 + aa2 * bb1;
+#else
+      cc1 = aa1 * bb1 + aa2 * bb2;
+      cc2 = aa1 * bb2 - aa2 * bb1;
+#endif
+
+      *(b + 0) = cc1;
+      *(b + 1) = cc2;
+      *(c + i * 2 + 0 + j * ldc) = cc1;
+      *(c + i * 2 + 1 + j * ldc) = cc2;
+      b += 2;
+
+      for (k = i + 1; k < m; k ++){
+#ifndef CONJ
+	*(c + k * 2 + 0 + j * ldc) -= cc1 * *(a + k * 2 + 0) - cc2 * *(a + k * 2 + 1);
+	*(c + k * 2 + 1 + j * ldc) -= cc1 * *(a + k * 2 + 1) + cc2 * *(a + k * 2 + 0);
+#else
+	*(c + k * 2 + 0 + j * ldc) -= cc1 * *(a + k * 2 + 0) + cc2 * *(a + k * 2 + 1);
+	*(c + k * 2 + 1 + j * ldc) -= -cc1 * *(a + k * 2 + 1) + cc2 * *(a + k * 2 + 0);
+#endif
+      }
+
+    }
+    a += m * 2;
+  }
+}
+
+#endif
+
+
+int CNAME(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT dummy1,
+#ifdef COMPLEX
+	   FLOAT dummy2,
+#endif
+	   FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG offset){
+
+  FLOAT *aa, *cc;
+  BLASLONG  kk;
+  BLASLONG  mmodM = m - (m / GEMM_UNROLL_M) * GEMM_UNROLL_M;
+  BLASLONG i, j, jj;
+
+#if 0
+  fprintf(stderr, "TRSM KERNEL LT : m = %3ld  n = %3ld  k = %3ld offset = %3ld\n",
+	  m, n, k, offset);
+#endif
+
+  jj = 0;
+
+  j = (n >> GEMM_UNROLL_N_SHIFT);
+
+  while (j > 0) {
+
+    kk = offset;
+    aa = a;
+    cc = c;
+
+    i = (m / GEMM_UNROLL_M);
+
+    while (i > 0) {
+
+	if (kk > 0) {
+	  GEMM_KERNEL(GEMM_UNROLL_M, GEMM_UNROLL_N, kk, dm1,
+#ifdef COMPLEX
+		      ZERO,
+#endif
+		      aa, b, cc, ldc);
+	}
+
+	solve(GEMM_UNROLL_M, GEMM_UNROLL_N,
+	      aa + kk * GEMM_UNROLL_M * COMPSIZE,
+	      b  + kk * GEMM_UNROLL_N * COMPSIZE,
+	      cc, ldc);
+
+      aa += GEMM_UNROLL_M * k * COMPSIZE;
+      cc += GEMM_UNROLL_M     * COMPSIZE;
+      kk += GEMM_UNROLL_M;
+      i --;
+    }
+
+    if (mmodM) {
+      i = 4; /* largest power-of-2 edge < UNROLL_M=6 */
+      while (i > 0) {
+	if (mmodM & i) {
+	    if (kk > 0) {
+	      GEMM_KERNEL(i, GEMM_UNROLL_N, kk, dm1,
+#ifdef COMPLEX
+			  ZERO,
+#endif
+			  aa, b, cc, ldc);
+	    }
+	  solve(i, GEMM_UNROLL_N,
+		aa + kk * i             * COMPSIZE,
+		b  + kk * GEMM_UNROLL_N * COMPSIZE,
+		cc, ldc);
+
+	  aa += i * k * COMPSIZE;
+	  cc += i     * COMPSIZE;
+	  kk += i;
+	}
+	i >>= 1;
+      }
+    }
+
+    b += GEMM_UNROLL_N * k   * COMPSIZE;
+    c += GEMM_UNROLL_N * ldc * COMPSIZE;
+    j --;
+    jj += GEMM_UNROLL_M;
+  }
+
+  if (n & (GEMM_UNROLL_N - 1)) {
+
+    j = (GEMM_UNROLL_N >> 1);
+    while (j > 0) {
+      if (n & j) {
+
+	kk = offset;
+	aa = a;
+	cc = c;
+
+	i = (m / GEMM_UNROLL_M);
+
+	while (i > 0) {
+	  if (kk > 0) {
+	    GEMM_KERNEL(GEMM_UNROLL_M, j, kk, dm1,
+#ifdef COMPLEX
+			ZERO,
+#endif
+			aa,
+			b,
+			cc,
+			ldc);
+	  }
+
+	  solve(GEMM_UNROLL_M, j,
+		aa + kk * GEMM_UNROLL_M * COMPSIZE,
+		b  + kk * j             * COMPSIZE, cc, ldc);
+
+	  aa += GEMM_UNROLL_M * k * COMPSIZE;
+	  cc += GEMM_UNROLL_M     * COMPSIZE;
+	  kk += GEMM_UNROLL_M;
+	  i --;
+	}
+
+	if (mmodM) {
+	  i = 4; /* largest power-of-2 edge < UNROLL_M=6 */
+	  while (i > 0) {
+	    if (mmodM & i) {
+	      if (kk > 0) {
+		GEMM_KERNEL(i, j, kk, dm1,
+#ifdef COMPLEX
+			    ZERO,
+#endif
+			    aa,
+			    b,
+			    cc,
+			    ldc);
+	      }
+
+	      solve(i, j,
+		    aa + kk * i * COMPSIZE,
+		    b  + kk * j * COMPSIZE, cc, ldc);
+
+	      aa += i * k * COMPSIZE;
+	      cc += i     * COMPSIZE;
+	      kk += i;
+	      }
+	    i >>= 1;
+	  }
+	}
+
+	b += j * k   * COMPSIZE;
+	c += j * ldc * COMPSIZE;
+      }
+      j >>= 1;
+    }
+  }
+
+  return 0;
+}
@@ -0,0 +1,328 @@
+/*********************************************************************/
+/* Copyright 2009, 2010 The University of Texas at Austin.           */
+/* All rights reserved.                                              */
+/*                                                                   */
+/* Redistribution and use in source and binary forms, with or        */
+/* without modification, are permitted provided that the following   */
+/* conditions are met:                                               */
+/*                                                                   */
+/*   1. Redistributions of source code must retain the above         */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer.                                                  */
+/*                                                                   */
+/*   2. Redistributions in binary form must reproduce the above      */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer in the documentation and/or other materials       */
+/*      provided with the distribution.                              */
+/*                                                                   */
+/*    THIS  SOFTWARE IS PROVIDED  BY THE  UNIVERSITY OF  TEXAS AT    */
+/*    AUSTIN  ``AS IS''  AND ANY  EXPRESS OR  IMPLIED WARRANTIES,    */
+/*    INCLUDING, BUT  NOT LIMITED  TO, THE IMPLIED  WARRANTIES OF    */
+/*    MERCHANTABILITY  AND FITNESS FOR  A PARTICULAR  PURPOSE ARE    */
+/*    DISCLAIMED.  IN  NO EVENT SHALL THE UNIVERSITY  OF TEXAS AT    */
+/*    AUSTIN OR CONTRIBUTORS BE  LIABLE FOR ANY DIRECT, INDIRECT,    */
+/*    INCIDENTAL,  SPECIAL, EXEMPLARY,  OR  CONSEQUENTIAL DAMAGES    */
+/*    (INCLUDING, BUT  NOT LIMITED TO,  PROCUREMENT OF SUBSTITUTE    */
+/*    GOODS  OR  SERVICES; LOSS  OF  USE,  DATA,  OR PROFITS;  OR    */
+/*    BUSINESS INTERRUPTION) HOWEVER CAUSED  AND ON ANY THEORY OF    */
+/*    LIABILITY, WHETHER  IN CONTRACT, STRICT  LIABILITY, OR TORT    */
+/*    (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY WAY OUT    */
+/*    OF  THE  USE OF  THIS  SOFTWARE,  EVEN  IF ADVISED  OF  THE    */
+/*    POSSIBILITY OF SUCH DAMAGE.                                    */
+/*                                                                   */
+/* The views and conclusions contained in the software and           */
+/* documentation are those of the authors and should not be          */
+/* interpreted as representing official policies, either expressed   */
+/* or implied, of The University of Texas at Austin.                 */
+/*********************************************************************/
+
+#include "common.h"
+
+/* UNROLL_M=6-aware TRSM kernel.
+ * Stock generic/trsm_kernel_*.c assumes power-of-2 UNROLL_M via
+ * GEMM_UNROLL_M_SHIFT (incorrectly set to 2 for M=6) and bit masks
+ * (m & (M-1)). That overruns packed buffers and corrupts the heap —
+ * reproduces as free()/double-free in HPL panel factor (dtrsm).
+ * Pattern mirrors loongarch64/trsm_kernel_*_UNROLLN6.c for N=6.
+ */
+
+static FLOAT dm1 = -1.;
+
+#ifdef CONJ
+#define GEMM_KERNEL   GEMM_KERNEL_R
+#else
+#define GEMM_KERNEL   GEMM_KERNEL_N
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 1
+#define GEMM_UNROLL_M_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 2
+#define GEMM_UNROLL_M_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 4
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 6
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 8
+#define GEMM_UNROLL_M_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 16
+#define GEMM_UNROLL_M_SHIFT 4
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 1
+#define GEMM_UNROLL_N_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 2
+#define GEMM_UNROLL_N_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 4
+#define GEMM_UNROLL_N_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 8
+#define GEMM_UNROLL_N_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 16
+#define GEMM_UNROLL_N_SHIFT 4
+#endif
+
+#ifndef COMPLEX
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa, bb;
+
+  int i, j, k;
+
+  for (i = 0; i < n; i++) {
+
+    bb = *(b + i);
+
+    for (j = 0; j < m; j ++) {
+      aa = *(c + j + i * ldc);
+      aa *= bb;
+      *a  = aa;
+      *(c + j + i * ldc) = aa;
+      a ++;
+
+      for (k = i + 1; k < n; k ++){
+	*(c + j + k * ldc) -= aa * *(b + k);
+      }
+
+    }
+    b += n;
+  }
+}
+
+#else
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa1, aa2;
+  FLOAT bb1, bb2;
+  FLOAT cc1, cc2;
+
+  int i, j, k;
+
+  ldc *= 2;
+
+  for (i = 0; i < n; i++) {
+
+    bb1 = *(b + i * 2 + 0);
+    bb2 = *(b + i * 2 + 1);
+
+    for (j = 0; j < m; j ++) {
+      aa1 = *(c + j * 2 + 0 + i * ldc);
+      aa2 = *(c + j * 2 + 1 + i * ldc);
+
+#ifndef CONJ
+      cc1 = aa1 * bb1 - aa2 * bb2;
+      cc2 = aa1 * bb2 + aa2 * bb1;
+#else
+      cc1 =  aa1 * bb1 + aa2 * bb2;
+      cc2 = -aa1 * bb2 + aa2 * bb1;
+#endif
+
+      *(a + 0) = cc1;
+      *(a + 1) = cc2;
+      *(c + j * 2 + 0 + i * ldc) = cc1;
+      *(c + j * 2 + 1 + i * ldc) = cc2;
+      a += 2;
+
+      for (k = i + 1; k < n; k ++){
+#ifndef CONJ
+	*(c + j * 2 + 0 + k * ldc) -= cc1 * *(b + k * 2 + 0) - cc2 * *(b + k * 2 + 1);
+	*(c + j * 2 + 1 + k * ldc) -= cc1 * *(b + k * 2 + 1) + cc2 * *(b + k * 2 + 0);
+#else
+	*(c + j * 2 + 0 + k * ldc) -=   cc1 * *(b + k * 2 + 0) + cc2 * *(b + k * 2 + 1);
+	*(c + j * 2 + 1 + k * ldc) -= - cc1 * *(b + k * 2 + 1) + cc2 * *(b + k * 2 + 0);
+#endif
+      }
+
+    }
+    b += n * 2;
+  }
+}
+
+#endif
+
+
+int CNAME(BLASLONG m, BLASLONG n, BLASLONG k, FLOAT dummy1,
+#ifdef COMPLEX
+	   FLOAT dummy2,
+#endif
+	   FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG offset){
+
+  FLOAT *aa, *cc;
+  BLASLONG  kk;
+  BLASLONG  mmodM = m - (m / GEMM_UNROLL_M) * GEMM_UNROLL_M;
+  BLASLONG i, j, jj;
+
+#if 0
+  fprintf(stderr, "TRSM RN KERNEL m = %3ld  n = %3ld  k = %3ld offset = %3ld\n",
+	  m, n, k, offset);
+#endif
+
+  jj = 0;
+  j = (n >> GEMM_UNROLL_N_SHIFT);
+  kk = -offset;
+
+  while (j > 0) {
+
+    aa = a;
+    cc = c;
+
+    i = (m / GEMM_UNROLL_M);
+
+    if (i > 0) {
+      do {
+	if (kk > 0) {
+	  GEMM_KERNEL(GEMM_UNROLL_M, GEMM_UNROLL_N, kk, dm1,
+#ifdef COMPLEX
+		      ZERO,
+#endif
+		      aa, b, cc, ldc);
+	}
+
+	solve(GEMM_UNROLL_M, GEMM_UNROLL_N,
+	      aa + kk * GEMM_UNROLL_M * COMPSIZE,
+	      b  + kk * GEMM_UNROLL_N * COMPSIZE,
+	      cc, ldc);
+
+	aa += GEMM_UNROLL_M * k * COMPSIZE;
+	cc += GEMM_UNROLL_M     * COMPSIZE;
+	i --;
+      } while (i > 0);
+    }
+
+
+    if (mmodM) {
+      i = 4; /* largest power-of-2 edge < UNROLL_M=6 */
+      while (i > 0) {
+	if (mmodM & i) {
+	    if (kk > 0) {
+	      GEMM_KERNEL(i, GEMM_UNROLL_N, kk, dm1,
+#ifdef COMPLEX
+			  ZERO,
+#endif
+			  aa, b, cc, ldc);
+	    }
+	  solve(i, GEMM_UNROLL_N,
+		aa + kk * i             * COMPSIZE,
+		b  + kk * GEMM_UNROLL_N * COMPSIZE,
+		cc, ldc);
+
+	  aa += i * k * COMPSIZE;
+	  cc += i     * COMPSIZE;
+	}
+	i >>= 1;
+      }
+    }
+
+    kk += GEMM_UNROLL_N;
+    b += GEMM_UNROLL_N * k   * COMPSIZE;
+    c += GEMM_UNROLL_N * ldc * COMPSIZE;
+    j --;
+    jj += GEMM_UNROLL_M;
+  }
+
+  if (n & (GEMM_UNROLL_N - 1)) {
+
+    j = (GEMM_UNROLL_N >> 1);
+    while (j > 0) {
+      if (n & j) {
+
+	aa = a;
+	cc = c;
+
+	i = (m / GEMM_UNROLL_M);
+
+	while (i > 0) {
+	  if (kk > 0) {
+	    GEMM_KERNEL(GEMM_UNROLL_M, j, kk, dm1,
+#ifdef COMPLEX
+			ZERO,
+#endif
+			aa,
+			b,
+			cc,
+			ldc);
+	  }
+
+	  solve(GEMM_UNROLL_M, j,
+		aa + kk * GEMM_UNROLL_M * COMPSIZE,
+		b  + kk * j             * COMPSIZE, cc, ldc);
+
+	  aa += GEMM_UNROLL_M * k * COMPSIZE;
+	  cc += GEMM_UNROLL_M     * COMPSIZE;
+	  i --;
+	}
+
+	if (mmodM) {
+	  i = 4; /* largest power-of-2 edge < UNROLL_M=6 */
+	  while (i > 0) {
+	    if (mmodM & i) {
+	      if (kk > 0) {
+		GEMM_KERNEL(i, j, kk, dm1,
+#ifdef COMPLEX
+			    ZERO,
+#endif
+			    aa,
+			    b,
+			    cc,
+			    ldc);
+	      }
+
+	      solve(i, j,
+		    aa + kk * i * COMPSIZE,
+		    b  + kk * j * COMPSIZE, cc, ldc);
+
+	      aa += i * k * COMPSIZE;
+	      cc += i     * COMPSIZE;
+	      }
+	    i >>= 1;
+	  }
+	}
+
+	b += j * k   * COMPSIZE;
+	c += j * ldc * COMPSIZE;
+	kk += j;
+      }
+      j >>= 1;
+    }
+  }
+
+  return 0;
+}
@@ -0,0 +1,355 @@
+/*********************************************************************/
+/* Copyright 2009, 2010 The University of Texas at Austin.           */
+/* All rights reserved.                                              */
+/*                                                                   */
+/* Redistribution and use in source and binary forms, with or        */
+/* without modification, are permitted provided that the following   */
+/* conditions are met:                                               */
+/*                                                                   */
+/*   1. Redistributions of source code must retain the above         */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer.                                                  */
+/*                                                                   */
+/*   2. Redistributions in binary form must reproduce the above      */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer in the documentation and/or other materials       */
+/*      provided with the distribution.                              */
+/*                                                                   */
+/*    THIS  SOFTWARE IS PROVIDED  BY THE  UNIVERSITY OF  TEXAS AT    */
+/*    AUSTIN  ``AS IS''  AND ANY  EXPRESS OR  IMPLIED WARRANTIES,    */
+/*    INCLUDING, BUT  NOT LIMITED  TO, THE IMPLIED  WARRANTIES OF    */
+/*    MERCHANTABILITY  AND FITNESS FOR  A PARTICULAR  PURPOSE ARE    */
+/*    DISCLAIMED.  IN  NO EVENT SHALL THE UNIVERSITY  OF TEXAS AT    */
+/*    AUSTIN OR CONTRIBUTORS BE  LIABLE FOR ANY DIRECT, INDIRECT,    */
+/*    INCIDENTAL,  SPECIAL, EXEMPLARY,  OR  CONSEQUENTIAL DAMAGES    */
+/*    (INCLUDING, BUT  NOT LIMITED TO,  PROCUREMENT OF SUBSTITUTE    */
+/*    GOODS  OR  SERVICES; LOSS  OF  USE,  DATA,  OR PROFITS;  OR    */
+/*    BUSINESS INTERRUPTION) HOWEVER CAUSED  AND ON ANY THEORY OF    */
+/*    LIABILITY, WHETHER  IN CONTRACT, STRICT  LIABILITY, OR TORT    */
+/*    (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY WAY OUT    */
+/*    OF  THE  USE OF  THIS  SOFTWARE,  EVEN  IF ADVISED  OF  THE    */
+/*    POSSIBILITY OF SUCH DAMAGE.                                    */
+/*                                                                   */
+/* The views and conclusions contained in the software and           */
+/* documentation are those of the authors and should not be          */
+/* interpreted as representing official policies, either expressed   */
+/* or implied, of The University of Texas at Austin.                 */
+/*********************************************************************/
+
+#include "common.h"
+
+/* UNROLL_M=6-aware TRSM kernel.
+ * Stock generic/trsm_kernel_*.c assumes power-of-2 UNROLL_M via
+ * GEMM_UNROLL_M_SHIFT (incorrectly set to 2 for M=6) and bit masks
+ * (m & (M-1)). That overruns packed buffers and corrupts the heap —
+ * reproduces as free()/double-free in HPL panel factor (dtrsm).
+ * Pattern mirrors loongarch64/trsm_kernel_*_UNROLLN6.c for N=6.
+ */
+
+static FLOAT dm1 = -1.;
+
+#ifdef CONJ
+#define GEMM_KERNEL   GEMM_KERNEL_R
+#else
+#define GEMM_KERNEL   GEMM_KERNEL_N
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 1
+#define GEMM_UNROLL_M_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 2
+#define GEMM_UNROLL_M_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 4
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 6
+#define GEMM_UNROLL_M_SHIFT 2
+#endif
+
+
+#if GEMM_DEFAULT_UNROLL_M == 8
+#define GEMM_UNROLL_M_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_M == 16
+#define GEMM_UNROLL_M_SHIFT 4
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 1
+#define GEMM_UNROLL_N_SHIFT 0
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 2
+#define GEMM_UNROLL_N_SHIFT 1
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 4
+#define GEMM_UNROLL_N_SHIFT 2
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 8
+#define GEMM_UNROLL_N_SHIFT 3
+#endif
+
+#if GEMM_DEFAULT_UNROLL_N == 16
+#define GEMM_UNROLL_N_SHIFT 4
+#endif
+
+
+#ifndef COMPLEX
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa,  bb;
+
+  int i, j, k;
+
+  a += (n - 1) * m;
+  b += (n - 1) * n;
+
+  for (i = n - 1; i >= 0; i--) {
+
+    bb = *(b + i);
+
+    for (j = 0; j < m; j ++) {
+      aa = *(c + j + i * ldc);
+      aa *= bb;
+      *a   = aa;
+      *(c + j + i * ldc) = aa;
+      a ++;
+
+      for (k = 0; k < i; k ++){
+	*(c + j + k * ldc) -= aa * *(b + k);
+      }
+
+    }
+    b -= n;
+    a -= 2 * m;
+  }
+
+}
+
+#else
+
+static inline void solve(BLASLONG m, BLASLONG n, FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc) {
+
+  FLOAT aa1, aa2;
+  FLOAT bb1, bb2;
+  FLOAT cc1, cc2;
+
+  int i, j, k;
+
+  ldc *= 2;
+
+  a += (n - 1) * m * 2;
+  b += (n - 1) * n * 2;
+
+  for (i = n - 1; i >= 0; i--) {
+
+    bb1 = *(b + i * 2 + 0);
+    bb2 = *(b + i * 2 + 1);
+
+    for (j = 0; j < m; j ++) {
+
+      aa1 = *(c + j * 2 + 0 + i * ldc);
+      aa2 = *(c + j * 2 + 1 + i * ldc);
+
+#ifndef CONJ
+      cc1 = aa1 * bb1 - aa2 * bb2;
+      cc2 = aa1 * bb2 + aa2 * bb1;
+#else
+      cc1 =  aa1 * bb1  + aa2 * bb2;
+      cc2 = - aa1 * bb2 + aa2 * bb1;
+#endif
+
+      *(a + 0) = cc1;
+      *(a + 1) = cc2;
+
+      *(c + j * 2 + 0 + i * ldc) = cc1;
+      *(c + j * 2 + 1 + i * ldc) = cc2;
+      a += 2;
+
+      for (k = 0; k < i; k ++){
+#ifndef CONJ
+	*(c + j * 2 + 0 + k * ldc) -= cc1 * *(b + k * 2 + 0) - cc2 * *(b + k * 2 + 1);
+	*(c + j * 2 + 1 + k * ldc) -= cc1 * *(b + k * 2 + 1) + cc2 * *(b + k * 2 + 0);
+#else
+	*(c + j * 2 + 0 + k * ldc) -=   cc1 * *(b + k * 2 + 0) + cc2 * *(b + k * 2 + 1);
+	*(c + j * 2 + 1 + k * ldc) -=  -cc1 * *(b + k * 2 + 1) + cc2 * *(b + k * 2 + 0);
+#endif
+      }
+
+    }
+    b -= n * 2;
+    a -= 4 * m;
+  }
+
+}
+
+#endif
+
+int CNAME(BLASLONG m, BLASLONG n, BLASLONG k,  FLOAT dummy1,
+#ifdef COMPLEX
+	   FLOAT dummy2,
+#endif
+	   FLOAT *a, FLOAT *b, FLOAT *c, BLASLONG ldc, BLASLONG offset){
+
+  BLASLONG i, j;
+  FLOAT *aa, *cc;
+  BLASLONG  kk;
+  BLASLONG  mmodM = m - (m / GEMM_UNROLL_M) * GEMM_UNROLL_M;
+
+#if 0
+  fprintf(stderr, "TRSM RT KERNEL m = %3ld  n = %3ld  k = %3ld offset = %3ld\n",
+	  m, n, k, offset);
+#endif
+
+  kk = n - offset;
+  c += n * ldc * COMPSIZE;
+  b += n * k   * COMPSIZE;
+
+  if (n & (GEMM_UNROLL_N - 1)) {
+
+    j = 1;
+    while (j < GEMM_UNROLL_N) {
+      if (n & j) {
+
+	aa  = a;
+	b -= j * k  * COMPSIZE;
+	c -= j * ldc* COMPSIZE;
+	cc  = c;
+
+	i = (m / GEMM_UNROLL_M);
+	if (i > 0) {
+
+	  do {
+	    if (k - kk > 0) {
+	      GEMM_KERNEL(GEMM_UNROLL_M, j, k - kk, dm1,
+#ifdef COMPLEX
+			  ZERO,
+#endif
+			  aa + GEMM_UNROLL_M * kk * COMPSIZE,
+			  b  +  j            * kk * COMPSIZE,
+			  cc,
+			  ldc);
+	    }
+
+	    solve(GEMM_UNROLL_M, j,
+		  aa + (kk - j) * GEMM_UNROLL_M * COMPSIZE,
+		  b  + (kk - j) * j             * COMPSIZE,
+		  cc, ldc);
+
+	    aa += GEMM_UNROLL_M * k * COMPSIZE;
+	    cc += GEMM_UNROLL_M     * COMPSIZE;
+	    i --;
+	  } while (i > 0);
+	}
+
+	if (mmodM) {
+	  i = 4; /* largest power-of-2 edge < UNROLL_M=6 */
+	  do {
+	    if (mmodM & i) {
+
+	      if (k - kk > 0) {
+		GEMM_KERNEL(i, j, k - kk, dm1,
+#ifdef COMPLEX
+			    ZERO,
+#endif
+			    aa + i * kk * COMPSIZE,
+			    b  + j * kk * COMPSIZE,
+			    cc, ldc);
+	      }
+
+	      solve(i, j,
+		    aa + (kk - j) * i * COMPSIZE,
+		    b  + (kk - j) * j * COMPSIZE,
+		    cc, ldc);
+
+	      aa += i * k * COMPSIZE;
+	      cc += i     * COMPSIZE;
+
+	    }
+	    i >>= 1;
+	  } while (i > 0);
+	}
+	kk -= j;
+      }
+      j <<= 1;
+    }
+  }
+
+  j = (n >> GEMM_UNROLL_N_SHIFT);
+
+  if (j > 0) {
+
+    do {
+      aa  = a;
+      b -= GEMM_UNROLL_N * k   * COMPSIZE;
+      c -= GEMM_UNROLL_N * ldc * COMPSIZE;
+      cc  = c;
+
+      i = (m / GEMM_UNROLL_M);
+      if (i > 0) {
+	do {
+	  if (k - kk > 0) {
+	    GEMM_KERNEL(GEMM_UNROLL_M, GEMM_UNROLL_N, k - kk, dm1,
+#ifdef COMPLEX
+			ZERO,
+#endif
+			aa + GEMM_UNROLL_M * kk * COMPSIZE,
+			b  + GEMM_UNROLL_N * kk * COMPSIZE,
+			cc,
+			ldc);
+	  }
+
+	  solve(GEMM_UNROLL_M, GEMM_UNROLL_N,
+		aa + (kk - GEMM_UNROLL_N) * GEMM_UNROLL_M * COMPSIZE,
+		b  + (kk - GEMM_UNROLL_N) * GEMM_UNROLL_N * COMPSIZE,
+		cc, ldc);
+
+	  aa += GEMM_UNROLL_M * k * COMPSIZE;
+	  cc += GEMM_UNROLL_M     * COMPSIZE;
+	  i --;
+	} while (i > 0);
+      }
+
+      if (mmodM) {
+	i = 4; /* largest power-of-2 edge < UNROLL_M=6 */
+	do {
+	  if (mmodM & i) {
+	    if (k - kk > 0) {
+	      GEMM_KERNEL(i, GEMM_UNROLL_N, k - kk, dm1,
+#ifdef COMPLEX
+			  ZERO,
+#endif
+			  aa + i             * kk * COMPSIZE,
+			  b  + GEMM_UNROLL_N * kk * COMPSIZE,
+			  cc,
+			  ldc);
+	    }
+
+	    solve(i, GEMM_UNROLL_N,
+		  aa + (kk - GEMM_UNROLL_N) * i             * COMPSIZE,
+		  b  + (kk - GEMM_UNROLL_N) * GEMM_UNROLL_N * COMPSIZE,
+		  cc, ldc);
+
+	    aa += i * k * COMPSIZE;
+	    cc += i     * COMPSIZE;
+	  }
+	  i >>= 1;
+	} while (i > 0);
+      }
+
+      kk -= GEMM_UNROLL_N;
+      j --;
+    } while (j > 0);
+  }
+
+  return 0;
+}
+
+
@@ -0,0 +1,445 @@
+/*********************************************************************/
+/* Copyright 2009, 2010 The University of Texas at Austin.           */
+/* All rights reserved.                                              */
+/*                                                                   */
+/* Redistribution and use in source and binary forms, with or        */
+/* without modification, are permitted provided that the following   */
+/* conditions are met:                                               */
+/*                                                                   */
+/*   1. Redistributions of source code must retain the above         */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer.                                                  */
+/*                                                                   */
+/*   2. Redistributions in binary form must reproduce the above      */
+/*      copyright notice, this list of conditions and the following  */
+/*      disclaimer in the documentation and/or other materials       */
+/*      provided with the distribution.                              */
+/*                                                                   */
+/*    THIS  SOFTWARE IS PROVIDED  BY THE  UNIVERSITY OF  TEXAS AT    */
+/*    AUSTIN  ``AS IS''  AND ANY  EXPRESS OR  IMPLIED WARRANTIES,    */
+/*    INCLUDING, BUT  NOT LIMITED  TO, THE IMPLIED  WARRANTIES OF    */
+/*    MERCHANTABILITY  AND FITNESS FOR  A PARTICULAR  PURPOSE ARE    */
+/*    DISCLAIMED.  IN  NO EVENT SHALL THE UNIVERSITY  OF TEXAS AT    */
+/*    AUSTIN OR CONTRIBUTORS BE  LIABLE FOR ANY DIRECT, INDIRECT,    */
+/*    INCIDENTAL,  SPECIAL, EXEMPLARY,  OR  CONSEQUENTIAL DAMAGES    */
+/*    (INCLUDING, BUT  NOT LIMITED TO,  PROCUREMENT OF SUBSTITUTE    */
+/*    GOODS  OR  SERVICES; LOSS  OF  USE,  DATA,  OR PROFITS;  OR    */
+/*    BUSINESS INTERRUPTION) HOWEVER CAUSED  AND ON ANY THEORY OF    */
+/*    LIABILITY, WHETHER  IN CONTRACT, STRICT  LIABILITY, OR TORT    */
+/*    (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY WAY OUT    */
+/*    OF  THE  USE OF  THIS  SOFTWARE,  EVEN  IF ADVISED  OF  THE    */
+/*    POSSIBILITY OF SUCH DAMAGE.                                    */
+/*                                                                   */
+/* The views and conclusions contained in the software and           */
+/* documentation are those of the authors and should not be          */
+/* interpreted as representing official policies, either expressed   */
+/* or implied, of The University of Texas at Austin.                 */
+/*********************************************************************/
+
+#include <stdio.h>
+#include "common.h"
+
+/* Contiguous MR=6 negating transpose pack (needed when DGEMM_UNROLL_M=6).
+ * Layout matches the Goto-style panel packing of neg_tcopy_4/8: full
+ * 6-row panels first, then m remainders of 4 / 2 / 1. */
+
+int CNAME(BLASLONG m, BLASLONG n, FLOAT *a, BLASLONG lda, FLOAT *b){
+
+  BLASLONG i, j;
+
+  FLOAT *a_offset, *a_offset1, *a_offset2, *a_offset3, *a_offset4;
+  FLOAT *a_offset5, *a_offset6;
+  FLOAT *b_offset, *b_offset1, *b_offset2, *b_offset3;
+  FLOAT  ctemp1,  ctemp2,  ctemp3,  ctemp4;
+  FLOAT  ctemp5,  ctemp6,  ctemp7,  ctemp8;
+  FLOAT  ctemp9, ctemp10, ctemp11, ctemp12;
+  FLOAT ctemp13, ctemp14, ctemp15, ctemp16;
+  FLOAT ctemp17, ctemp18, ctemp19, ctemp20;
+  FLOAT ctemp21, ctemp22, ctemp23, ctemp24;
+
+  a_offset   = a;
+  b_offset   = b;
+
+  b_offset2  = b + m  * (n & ~3);
+  b_offset3  = b + m  * (n & ~1);
+
+  j = (m / 6);
+  if (j > 0){
+    do{
+      a_offset1  = a_offset;
+      a_offset2  = a_offset1 + lda;
+      a_offset3  = a_offset2 + lda;
+      a_offset4  = a_offset3 + lda;
+      a_offset5  = a_offset4 + lda;
+      a_offset6  = a_offset5 + lda;
+      a_offset  += 6 * lda;
+
+      b_offset1  = b_offset;
+      b_offset  += 24;
+
+      i = (n >> 2);
+      if (i > 0){
+	do{
+	  ctemp1  = *(a_offset1 + 0);
+	  ctemp2  = *(a_offset1 + 1);
+	  ctemp3  = *(a_offset1 + 2);
+	  ctemp4  = *(a_offset1 + 3);
+
+	  ctemp5  = *(a_offset2 + 0);
+	  ctemp6  = *(a_offset2 + 1);
+	  ctemp7  = *(a_offset2 + 2);
+	  ctemp8  = *(a_offset2 + 3);
+
+	  ctemp9  = *(a_offset3 + 0);
+	  ctemp10 = *(a_offset3 + 1);
+	  ctemp11 = *(a_offset3 + 2);
+	  ctemp12 = *(a_offset3 + 3);
+
+	  ctemp13 = *(a_offset4 + 0);
+	  ctemp14 = *(a_offset4 + 1);
+	  ctemp15 = *(a_offset4 + 2);
+	  ctemp16 = *(a_offset4 + 3);
+
+	  ctemp17 = *(a_offset5 + 0);
+	  ctemp18 = *(a_offset5 + 1);
+	  ctemp19 = *(a_offset5 + 2);
+	  ctemp20 = *(a_offset5 + 3);
+
+	  ctemp21 = *(a_offset6 + 0);
+	  ctemp22 = *(a_offset6 + 1);
+	  ctemp23 = *(a_offset6 + 2);
+	  ctemp24 = *(a_offset6 + 3);
+
+	  a_offset1 += 4;
+	  a_offset2 += 4;
+	  a_offset3 += 4;
+	  a_offset4 += 4;
+	  a_offset5 += 4;
+	  a_offset6 += 4;
+
+	  *(b_offset1 +  0) = -ctemp1;
+	  *(b_offset1 +  1) = -ctemp2;
+	  *(b_offset1 +  2) = -ctemp3;
+	  *(b_offset1 +  3) = -ctemp4;
+
+	  *(b_offset1 +  4) = -ctemp5;
+	  *(b_offset1 +  5) = -ctemp6;
+	  *(b_offset1 +  6) = -ctemp7;
+	  *(b_offset1 +  7) = -ctemp8;
+
+	  *(b_offset1 +  8) = -ctemp9;
+	  *(b_offset1 +  9) = -ctemp10;
+	  *(b_offset1 + 10) = -ctemp11;
+	  *(b_offset1 + 11) = -ctemp12;
+
+	  *(b_offset1 + 12) = -ctemp13;
+	  *(b_offset1 + 13) = -ctemp14;
+	  *(b_offset1 + 14) = -ctemp15;
+	  *(b_offset1 + 15) = -ctemp16;
+
+	  *(b_offset1 + 16) = -ctemp17;
+	  *(b_offset1 + 17) = -ctemp18;
+	  *(b_offset1 + 18) = -ctemp19;
+	  *(b_offset1 + 19) = -ctemp20;
+
+	  *(b_offset1 + 20) = -ctemp21;
+	  *(b_offset1 + 21) = -ctemp22;
+	  *(b_offset1 + 22) = -ctemp23;
+	  *(b_offset1 + 23) = -ctemp24;
+
+	  b_offset1 += m * 4;
+	  i --;
+	}while(i > 0);
+      }
+
+      if (n & 2) {
+	ctemp1  = *(a_offset1 + 0);
+	ctemp2  = *(a_offset1 + 1);
+
+	ctemp3  = *(a_offset2 + 0);
+	ctemp4  = *(a_offset2 + 1);
+
+	ctemp5  = *(a_offset3 + 0);
+	ctemp6  = *(a_offset3 + 1);
+
+	ctemp7  = *(a_offset4 + 0);
+	ctemp8  = *(a_offset4 + 1);
+
+	ctemp9  = *(a_offset5 + 0);
+	ctemp10 = *(a_offset5 + 1);
+
+	ctemp11 = *(a_offset6 + 0);
+	ctemp12 = *(a_offset6 + 1);
+
+	a_offset1 += 2;
+	a_offset2 += 2;
+	a_offset3 += 2;
+	a_offset4 += 2;
+	a_offset5 += 2;
+	a_offset6 += 2;
+
+	*(b_offset2 +  0) = -ctemp1;
+	*(b_offset2 +  1) = -ctemp2;
+	*(b_offset2 +  2) = -ctemp3;
+	*(b_offset2 +  3) = -ctemp4;
+
+	*(b_offset2 +  4) = -ctemp5;
+	*(b_offset2 +  5) = -ctemp6;
+	*(b_offset2 +  6) = -ctemp7;
+	*(b_offset2 +  7) = -ctemp8;
+
+	*(b_offset2 +  8) = -ctemp9;
+	*(b_offset2 +  9) = -ctemp10;
+	*(b_offset2 + 10) = -ctemp11;
+	*(b_offset2 + 11) = -ctemp12;
+
+	b_offset2 += 12;
+      }
+
+      if (n & 1) {
+	ctemp1  = *(a_offset1 + 0);
+	ctemp2  = *(a_offset2 + 0);
+	ctemp3  = *(a_offset3 + 0);
+	ctemp4  = *(a_offset4 + 0);
+	ctemp5  = *(a_offset5 + 0);
+	ctemp6  = *(a_offset6 + 0);
+
+	*(b_offset3 +  0) = -ctemp1;
+	*(b_offset3 +  1) = -ctemp2;
+	*(b_offset3 +  2) = -ctemp3;
+	*(b_offset3 +  3) = -ctemp4;
+	*(b_offset3 +  4) = -ctemp5;
+	*(b_offset3 +  5) = -ctemp6;
+
+	b_offset3 += 6;
+      }
+
+      j--;
+    }while(j > 0);
+  }
+
+  /* Remainders after m/6 panels are 0..5 — cascade 4, then 2, then 1.
+   * Do not use (m & 2)/(m & 1) on the original m: those bits are not the
+   * leftover after a non-power-of-2 main panel. */
+  j = m % 6;
+
+  if (j >= 4){
+    a_offset1  = a_offset;
+    a_offset2  = a_offset1 + lda;
+    a_offset3  = a_offset2 + lda;
+    a_offset4  = a_offset3 + lda;
+    a_offset  += 4 * lda;
+
+    b_offset1  = b_offset;
+    b_offset  += 16;
+
+    i = (n >> 2);
+    if (i > 0){
+      do{
+	ctemp1  = *(a_offset1 + 0);
+	ctemp2  = *(a_offset1 + 1);
+	ctemp3  = *(a_offset1 + 2);
+	ctemp4  = *(a_offset1 + 3);
+
+	ctemp5  = *(a_offset2 + 0);
+	ctemp6  = *(a_offset2 + 1);
+	ctemp7  = *(a_offset2 + 2);
+	ctemp8  = *(a_offset2 + 3);
+
+	ctemp9  = *(a_offset3 + 0);
+	ctemp10 = *(a_offset3 + 1);
+	ctemp11 = *(a_offset3 + 2);
+	ctemp12 = *(a_offset3 + 3);
+
+	ctemp13 = *(a_offset4 + 0);
+	ctemp14 = *(a_offset4 + 1);
+	ctemp15 = *(a_offset4 + 2);
+	ctemp16 = *(a_offset4 + 3);
+
+	a_offset1 += 4;
+	a_offset2 += 4;
+	a_offset3 += 4;
+	a_offset4 += 4;
+
+	*(b_offset1 +  0) = -ctemp1;
+	*(b_offset1 +  1) = -ctemp2;
+	*(b_offset1 +  2) = -ctemp3;
+	*(b_offset1 +  3) = -ctemp4;
+
+	*(b_offset1 +  4) = -ctemp5;
+	*(b_offset1 +  5) = -ctemp6;
+	*(b_offset1 +  6) = -ctemp7;
+	*(b_offset1 +  7) = -ctemp8;
+
+	*(b_offset1 +  8) = -ctemp9;
+	*(b_offset1 +  9) = -ctemp10;
+	*(b_offset1 + 10) = -ctemp11;
+	*(b_offset1 + 11) = -ctemp12;
+
+	*(b_offset1 + 12) = -ctemp13;
+	*(b_offset1 + 13) = -ctemp14;
+	*(b_offset1 + 14) = -ctemp15;
+	*(b_offset1 + 15) = -ctemp16;
+
+	b_offset1 += m * 4;
+	i --;
+      }while(i > 0);
+    }
+
+    if (n & 2) {
+      ctemp1  = *(a_offset1 + 0);
+      ctemp2  = *(a_offset1 + 1);
+
+      ctemp3  = *(a_offset2 + 0);
+      ctemp4  = *(a_offset2 + 1);
+
+      ctemp5  = *(a_offset3 + 0);
+      ctemp6  = *(a_offset3 + 1);
+
+      ctemp7  = *(a_offset4 + 0);
+      ctemp8  = *(a_offset4 + 1);
+
+      a_offset1 += 2;
+      a_offset2 += 2;
+      a_offset3 += 2;
+      a_offset4 += 2;
+
+      *(b_offset2 +  0) = -ctemp1;
+      *(b_offset2 +  1) = -ctemp2;
+      *(b_offset2 +  2) = -ctemp3;
+      *(b_offset2 +  3) = -ctemp4;
+
+      *(b_offset2 +  4) = -ctemp5;
+      *(b_offset2 +  5) = -ctemp6;
+      *(b_offset2 +  6) = -ctemp7;
+      *(b_offset2 +  7) = -ctemp8;
+
+      b_offset2 += 8;
+    }
+
+    if (n & 1) {
+      ctemp1  = *(a_offset1 + 0);
+      ctemp2  = *(a_offset2 + 0);
+      ctemp3  = *(a_offset3 + 0);
+      ctemp4  = *(a_offset4 + 0);
+
+      *(b_offset3 +  0) = -ctemp1;
+      *(b_offset3 +  1) = -ctemp2;
+      *(b_offset3 +  2) = -ctemp3;
+      *(b_offset3 +  3) = -ctemp4;
+
+      b_offset3 += 4;
+    }
+    j -= 4;
+  }
+
+  if (j >= 2){
+    a_offset1  = a_offset;
+    a_offset2  = a_offset1 + lda;
+    a_offset  += 2 * lda;
+
+    b_offset1  = b_offset;
+    b_offset  += 8;
+
+    i = (n >> 2);
+    if (i > 0){
+      do{
+	ctemp1  = *(a_offset1 + 0);
+	ctemp2  = *(a_offset1 + 1);
+	ctemp3  = *(a_offset1 + 2);
+	ctemp4  = *(a_offset1 + 3);
+
+	ctemp5  = *(a_offset2 + 0);
+	ctemp6  = *(a_offset2 + 1);
+	ctemp7  = *(a_offset2 + 2);
+	ctemp8  = *(a_offset2 + 3);
+
+	a_offset1 += 4;
+	a_offset2 += 4;
+
+	*(b_offset1 +  0) = -ctemp1;
+	*(b_offset1 +  1) = -ctemp2;
+	*(b_offset1 +  2) = -ctemp3;
+	*(b_offset1 +  3) = -ctemp4;
+
+	*(b_offset1 +  4) = -ctemp5;
+	*(b_offset1 +  5) = -ctemp6;
+	*(b_offset1 +  6) = -ctemp7;
+	*(b_offset1 +  7) = -ctemp8;
+
+	b_offset1 += m * 4;
+	i --;
+      }while(i > 0);
+    }
+
+    if (n & 2) {
+      ctemp1  = *(a_offset1 + 0);
+      ctemp2  = *(a_offset1 + 1);
+
+      ctemp3  = *(a_offset2 + 0);
+      ctemp4  = *(a_offset2 + 1);
+
+      a_offset1 += 2;
+      a_offset2 += 2;
+
+      *(b_offset2 +  0) = -ctemp1;
+      *(b_offset2 +  1) = -ctemp2;
+      *(b_offset2 +  2) = -ctemp3;
+      *(b_offset2 +  3) = -ctemp4;
+
+      b_offset2 += 4;
+    }
+
+    if (n & 1) {
+      ctemp1  = *(a_offset1 + 0);
+      ctemp2  = *(a_offset2 + 0);
+
+      *(b_offset3 +  0) = -ctemp1;
+      *(b_offset3 +  1) = -ctemp2;
+      b_offset3 += 2;
+    }
+    j -= 2;
+  }
+
+  if (j >= 1){
+    a_offset1  = a_offset;
+    b_offset1  = b_offset;
+
+    i = (n >> 2);
+    if (i > 0){
+      do{
+	ctemp1  = *(a_offset1 + 0);
+	ctemp2  = *(a_offset1 + 1);
+	ctemp3  = *(a_offset1 + 2);
+	ctemp4  = *(a_offset1 + 3);
+
+	a_offset1 += 4;
+
+	*(b_offset1 +  0) = -ctemp1;
+	*(b_offset1 +  1) = -ctemp2;
+	*(b_offset1 +  2) = -ctemp3;
+	*(b_offset1 +  3) = -ctemp4;
+
+	b_offset1 += 4 * m;
+
+	i --;
+      }while(i > 0);
+    }
+
+    if (n & 2) {
+      ctemp1  = *(a_offset1 + 0);
+      ctemp2  = *(a_offset1 + 1);
+      a_offset1 += 2;
+
+      *(b_offset2 +  0) = -ctemp1;
+      *(b_offset2 +  1) = -ctemp2;
+    }
+
+    if (n & 1) {
+      ctemp1  = *(a_offset1 + 0);
+      *(b_offset3 +  0) = -ctemp1;
+    }
+  }
+
+  return 0;
+}
@@ -143,4 +143,25 @@ ZTRSMKERNEL_LN  = ../generic/trsm_kernel_LN.c
 ZTRSMKERNEL_LT  = ../generic/trsm_kernel_LT.c
 ZTRSMKERNEL_RN  = ../generic/trsm_kernel_RN.c
 ZTRSMKERNEL_RT  = ../generic/trsm_kernel_RT.c
+
+SGEMM_SMALL_M_PERMIT = sgemm_small_matrix_permit.c
+SGEMM_SMALL_K_NN = sgemm_small_kernel_nn_lsx.S
+SGEMM_SMALL_K_B0_NN = sgemm_small_kernel_nn_lsx.S
+SGEMM_SMALL_K_NT = sgemm_small_kernel_nt_lsx.S
+SGEMM_SMALL_K_B0_NT = sgemm_small_kernel_nt_lsx.S
+SGEMM_SMALL_K_TN = sgemm_small_kernel_tn_lsx.S
+SGEMM_SMALL_K_B0_TN = sgemm_small_kernel_tn_lsx.S
+SGEMM_SMALL_K_TT = sgemm_small_kernel_tt_lsx.S
+SGEMM_SMALL_K_B0_TT = sgemm_small_kernel_tt_lsx.S
+
+DGEMM_SMALL_M_PERMIT = dgemm_small_matrix_permit.c
+DGEMM_SMALL_K_NN = dgemm_small_kernel_nn_lsx.S
+DGEMM_SMALL_K_B0_NN = dgemm_small_kernel_nn_lsx.S
+DGEMM_SMALL_K_NT = dgemm_small_kernel_nt_lsx.S
+DGEMM_SMALL_K_B0_NT = dgemm_small_kernel_nt_lsx.S
+DGEMM_SMALL_K_TN = dgemm_small_kernel_tn_lsx.S
+DGEMM_SMALL_K_B0_TN = dgemm_small_kernel_tn_lsx.S
+DGEMM_SMALL_K_TT = dgemm_small_kernel_tt_lsx.S
+DGEMM_SMALL_K_B0_TT = dgemm_small_kernel_tt_lsx.S
+
 endif
--- a/Show More
+++ b/Show More