Revert "Fix undefined CC in f_check (again)"

Fix undefined CC in f_check (again)

.drone.yml

@@ -190,3 +190,27 @@ steps:
     - make -C ctest $COMMON_FLAGS
     - make -C utest $COMMON_FLAGS
     - make -C cpp_thread_test dgemm_tester
+---
+kind: pipeline
+name: arm64_gcc10
+
+platform:
+  os: linux
+  arch: arm64
+
+steps:
+- name: Build and Test
+  image: ubuntu:20.04
+  environment:
+    CC: gcc-10
+    FC: gfortran-10
+    COMMON_FLAGS: 'TARGET=ARMV8 DYNAMIC_ARCH=1'
+  commands:
+    - echo "MAKE_FLAGS:= $COMMON_FLAGS"
+    - apt-get update -y
+    - apt-get install -y make $CC gfortran-10 perl python g++
+    - $CC --version
+    - make QUIET_MAKE=1 $COMMON_FLAGS
+    - make -C utest $COMMON_FLAGS
+    - make -C test $COMMON_FLAGS
+    

.github/workflows/nightly-Homebrew-build.yml

@@ -44,6 +44,11 @@ jobs:
         if: github.event_name != 'pull_request'
         run: brew update || true
 
+      - name: unlink installed gcc to allow updating
+        run: |
+          brew unlink gcc@8
+          brew unlink gcc@9
+          
       - name: Install prerequisites
         run: brew install --fetch-HEAD --HEAD --only-dependencies --keep-tmp openblas
 

.gitignore

@@ -91,3 +91,4 @@ benchmark/*.goto
 benchmark/smallscaling
 CMakeCache.txt
 CMakeFiles/*
+.vscode

Makefile.arm64

@@ -1,4 +1,4 @@
-
+ifneq ($(C_COMPILER), PGI)
 ifeq ($(CORE), ARMV8)
 CCOMMON_OPT += -march=armv8-a
 FCOMMON_OPT += -march=armv8-a
@@ -77,4 +77,4 @@ CCOMMON_OPT += -march=armv8.2-a -mtune=tsv110
 FCOMMON_OPT += -march=armv8.2-a -mtune=tsv110
 endif
 endif
-
+endif

Makefile.system

@@ -21,6 +21,8 @@ ifeq ($(ARCH), amd64)
 override ARCH=x86_64
 else ifeq ($(ARCH), powerpc64)
 override ARCH=power
+else ifeq ($(ARCH), powerpc64le)
+override ARCH=power
 else ifeq ($(ARCH), powerpc)
 override ARCH=power
 else ifeq ($(ARCH), i386)
@@ -623,6 +625,11 @@ DYNAMIC_CORE += THUNDERX2T99
 DYNAMIC_CORE += TSV110
 DYNAMIC_CORE += EMAG8180
 DYNAMIC_CORE += THUNDERX3T110
+ifdef DYNAMIC_LIST
+override DYNAMIC_CORE = ARMV8 $(DYNAMIC_LIST)
+XCCOMMON_OPT = -DDYNAMIC_LIST -DDYN_ARMV8
+XCCOMMON_OPT += $(foreach dcore,$(DYNAMIC_LIST),-DDYN_$(dcore))
+endif
 endif
 
 ifeq ($(ARCH), mips64)
@@ -1277,6 +1284,10 @@ CCOMMON_OPT	 += -DUSE_PAPI
 EXTRALIB	 += -lpapi -lperfctr
 endif
 
+ifdef BUFFERSIZE
+CCOMMON_OPT	 += -DBUFFERSIZE=$(BUFFERSIZE)
+endif
+
 ifdef DYNAMIC_THREADS
 CCOMMON_OPT	 += -DDYNAMIC_THREADS
 endif

README.md

@@ -212,7 +212,8 @@ Please note that it is not possible to combine support for different architectur
 - **Android**: Supported by the community. Please read <https://github.com/xianyi/OpenBLAS/wiki/How-to-build-OpenBLAS-for-Android>.
 - **AIX**: Supported on PPC up to POWER8
 - **Haiku**: Supported by the community. We don't actively test the library on this OS.
-- **SunOS**: Supported by the community. We don't actively test the library on this OS:
+- **SunOS**: Supported by the community. We don't actively test the library on this OS.
+- **Cortex-M**: Supported by the community. Please read <https://github.com/xianyi/OpenBLAS/wiki/How-to-use-OpenBLAS-on-Cortex-M>.
 
 ## Usage
 

appveyor.yml

@@ -30,10 +30,11 @@ environment:
     CONDA_INSTALL_LOCN: C:\\Miniconda36-x64
   matrix:
     - COMPILER: clang-cl
-      WITH_FORTRAN: yes
+      WITH_FORTRAN: ON
+      USE_OPENMP: ON
     - COMPILER: clang-cl
       DYNAMIC_ARCH: ON
-      WITH_FORTRAN: no
+      WITH_FORTRAN: OFF
     - COMPILER: cl
     - COMPILER: MinGW64-gcc-7.2.0-mingw
       DYNAMIC_ARCH: OFF
@@ -47,12 +48,7 @@ environment:
 install:
   - if [%COMPILER%]==[clang-cl] call %CONDA_INSTALL_LOCN%\Scripts\activate.bat
   - if [%COMPILER%]==[clang-cl] conda config --add channels conda-forge --force
-  - if [%COMPILER%]==[clang-cl] conda install --yes --quiet clangdev cmake
-
-  - if [%WITH_FORTRAN%]==[no] conda install --yes --quiet ninja
-  - if [%WITH_FORTRAN%]==[yes] conda install --yes --quiet -c isuruf kitware-ninja
-  - if [%WITH_FORTRAN%]==[yes] conda install --yes --quiet flang
-
+  - if [%COMPILER%]==[clang-cl] conda install --yes --quiet clangdev cmake ninja flang=11.0.1
   - if [%COMPILER%]==[clang-cl] call "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvarsall.bat" x64
   - if [%COMPILER%]==[clang-cl] set "LIB=%CONDA_INSTALL_LOCN%\Library\lib;%LIB%"
   - if [%COMPILER%]==[clang-cl] set "CPATH=%CONDA_INSTALL_LOCN%\Library\include;%CPATH%"
@@ -68,8 +64,9 @@ before_build:
   - if [%COMPILER%]==[MinGW64-gcc-7.2.0-mingw] cmake -G "MinGW Makefiles" -DNOFORTRAN=1 ..
   - if [%COMPILER%]==[MinGW-gcc-6.3.0-32] cmake -G "MSYS Makefiles" -DNOFORTRAN=1 ..
   - if [%COMPILER%]==[MinGW-gcc-5.3.0] cmake -G "MSYS Makefiles" -DNOFORTRAN=1 ..
-  - if [%WITH_FORTRAN%]==[no] cmake -G "Ninja" -DCMAKE_CXX_COMPILER=clang-cl -DCMAKE_C_COMPILER=clang-cl -DMSVC_STATIC_CRT=ON ..
-  - if [%WITH_FORTRAN%]==[yes] cmake -G "Ninja" -DCMAKE_CXX_COMPILER=clang-cl -DCMAKE_C_COMPILER=clang-cl -DCMAKE_Fortran_COMPILER=flang -DBUILD_WITHOUT_LAPACK=no -DNOFORTRAN=0 ..
+  - if [%WITH_FORTRAN%]==[OFF] cmake -G "Ninja" -DCMAKE_CXX_COMPILER=clang-cl -DCMAKE_C_COMPILER=clang-cl -DMSVC_STATIC_CRT=ON ..
+  - if [%WITH_FORTRAN%]==[ON] cmake -G "Ninja" -DCMAKE_CXX_COMPILER=clang-cl -DCMAKE_C_COMPILER=clang-cl -DCMAKE_Fortran_COMPILER=flang -DBUILD_WITHOUT_LAPACK=no -DNOFORTRAN=0 ..
+  - if [%USE_OPENMP%]==[ON] cmake -DUSE_OPENMP=ON ..
   - if [%DYNAMIC_ARCH%]==[ON] cmake -DDYNAMIC_ARCH=ON -DDYNAMIC_LIST='CORE2;NEHALEM;SANDYBRIDGE;BULLDOZER;HASWELL' ..
 
 build_script:

c_check

@@ -1,4 +1,4 @@
-#!/usr/bin/perl
+#!/usr/bin/env perl
 
 #use File::Basename;
 # use File::Temp qw(tempfile);

cblas.h

@@ -125,9 +125,14 @@ void cblas_zswap(OPENBLAS_CONST blasint n, void *x, OPENBLAS_CONST blasint incx,
 
 void cblas_srot(OPENBLAS_CONST blasint N, float *X, OPENBLAS_CONST blasint incX, float *Y, OPENBLAS_CONST blasint incY, OPENBLAS_CONST float c, OPENBLAS_CONST float s);
 void cblas_drot(OPENBLAS_CONST blasint N, double *X, OPENBLAS_CONST blasint incX, double *Y, OPENBLAS_CONST blasint incY, OPENBLAS_CONST double c, OPENBLAS_CONST double  s);
+void cblas_csrot(OPENBLAS_CONST blasint n, OPENBLAS_CONST void *x, OPENBLAS_CONST blasint incx, void *y, OPENBLAS_CONST blasint incY, OPENBLAS_CONST float c, OPENBLAS_CONST float s);
+void cblas_zdrot(OPENBLAS_CONST blasint n, OPENBLAS_CONST void *x, OPENBLAS_CONST blasint incx, void *y, OPENBLAS_CONST blasint incY, OPENBLAS_CONST double c, OPENBLAS_CONST double s);
 
 void cblas_srotg(float *a, float *b, float *c, float *s);
 void cblas_drotg(double *a, double *b, double *c, double *s);
+void cblas_crotg(void *a, void *b, float *c, void *s);
+void cblas_zrotg(void *a, void *b, double *c, void *s);
+
 
 void cblas_srotm(OPENBLAS_CONST blasint N, float *X, OPENBLAS_CONST blasint incX, float *Y, OPENBLAS_CONST blasint incY, OPENBLAS_CONST float *P);
 void cblas_drotm(OPENBLAS_CONST blasint N, double *X, OPENBLAS_CONST blasint incX, double *Y, OPENBLAS_CONST blasint incY, OPENBLAS_CONST double *P);

cmake/arch.cmake

@@ -45,6 +45,9 @@ endif ()
 if (DYNAMIC_ARCH)
   if (ARM64)
     set(DYNAMIC_CORE ARMV8 CORTEXA53 CORTEXA57 CORTEXA72 CORTEXA73 FALKOR THUNDERX THUNDERX2T99 TSV110 EMAG8180 NEOVERSEN1 THUNDERX3T110)
+    if (DYNAMIC_LIST)
+	    set(DYNAMIC_CORE ARMV8 ${DYNAMIC_LIST})
+    endif ()
   endif ()
   
   if (POWER)

cmake/lapacke.cmake

@@ -2499,6 +2499,5 @@ foreach (Utils_FILE ${Utils_SRC})
 endforeach ()
 
 set(lapacke_include_dir "${NETLIB_LAPACK_DIR}/LAPACKE/include")
-configure_file("${lapacke_include_dir}/lapacke_mangling_with_flags.h.in" "${lapacke_include_dir}/lapacke_mangling.h" COPYONLY)
 include_directories(${lapacke_include_dir})
 set_source_files_properties(${LAPACKE_SOURCES} PROPERTIES COMPILE_FLAGS "${LAPACK_CFLAGS}")

cmake/system.cmake

@@ -148,16 +148,20 @@ endif ()
 include("${PROJECT_SOURCE_DIR}/cmake/prebuild.cmake")
 if (DEFINED TARGET)
   if (${TARGET} STREQUAL COOPERLAKE AND NOT NO_AVX512)
-#    if (${CMAKE_C_COMPILER_ID} STREQUAL "GNU")
+    if (${CMAKE_C_COMPILER_ID} STREQUAL "GNU")
       execute_process(COMMAND ${CMAKE_C_COMPILER} -dumpversion OUTPUT_VARIABLE GCC_VERSION)
-        if (${GCC_VERSION} VERSION_GREATER 10.1 OR ${GCC_VERSION} VERSION_EQUAL 10.1)
+        if (${CMAKE_C_COMPILER_VERSION} VERSION_GREATER 10.09)
           set (KERNEL_DEFINITIONS "${KERNEL_DEFINITIONS} -march=cooperlake")
         else()
           set (KERNEL_DEFINITIONS "${KERNEL_DEFINITIONS} -march=skylake-avx512")
         endif()
-#    elseif (${CMAKE_C_COMPILER_ID} STREQUAL "CLANG")
-#      set (KERNEL_DEFINITIONS "${KERNEL_DEFINITIONS} -mavx2")
-#    endif()    
+    elseif (${CMAKE_C_COMPILER_ID} STREQUAL "Clang" OR ${CMAKE_C_COMPILER_ID} STREQUAL "AppleClang")
+         if (${CMAKE_C_COMPILER_VERSION} VERSION_GREATER 8.99)
+          set (KERNEL_DEFINITIONS "${KERNEL_DEFINITIONS} -march=cooperlake")
+        else()
+          set (KERNEL_DEFINITIONS "${KERNEL_DEFINITIONS} -march=skylake-avx512")
+        endif()
+    endif()    
   endif()
   if (${TARGET} STREQUAL SKYLAKEX AND NOT NO_AVX512)
     set (KERNEL_DEFINITIONS "${KERNEL_DEFINITIONS} -march=skylake-avx512")
@@ -233,6 +237,11 @@ if (BINARY64)
   endif ()
 endif ()
 
+if(EMBEDDED)
+  set(CCOMMON_OPT "${CCOMMON_OPT} -DOS_EMBEDDED")
+  set(CCOMMON_OPT "${CCOMMON_OPT} -mthumb -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16")
+endif()
+
 if (NEED_PIC)
   if (${CMAKE_C_COMPILER} STREQUAL "IBM")
     set(CCOMMON_OPT "${CCOMMON_OPT} -qpic=large")

cmake/utils.cmake

@@ -74,6 +74,9 @@ macro(ParseMakefileVars MAKEFILE_IN)
           string(REGEX MATCH "ifneq \\(\\$\\(([_A-Z]+)\\),[ \t]*([0-9_A-Z]+)\\)" line_match "${makefile_line}")
           if (NOT "${line_match}" STREQUAL "")
 #            message(STATUS "IFNEQ: ${line_match} first: ${CMAKE_MATCH_1} second: ${CMAKE_MATCH_2}")
+            if ( ${CMAKE_MATCH_1} STREQUAL C_COMPILER)
+                 set (CMAKE_MATCH_1 CMAKE_C_COMPILER)
+            endif ()
             if (NOT ( ${${CMAKE_MATCH_1}} STREQUAL ${CMAKE_MATCH_2}))
 #              message (STATUS "condition is true")
               set (IfElse 1)

common.h

@@ -122,7 +122,7 @@ extern "C" {
 #define ATOM GOTO_ATOM
 #undef  GOTO_ATOM
 #endif
-#else
+#elif !defined(OS_EMBEDDED)
 #include <sys/mman.h>
 #ifndef NO_SYSV_IPC
 #include <sys/shm.h>
@@ -134,6 +134,9 @@ extern "C" {
 #if defined(SMP) || defined(USE_LOCKING)
 #include <pthread.h>
 #endif
+#else
+#include <time.h>
+#include <math.h>
 #endif
 
 #if defined(OS_SUNOS)
@@ -488,10 +491,12 @@ static inline unsigned long long rpcc(void){
   struct timespec ts;
   clock_gettime(CLOCK_MONOTONIC, &ts);
   return (unsigned long long)ts.tv_sec * 1000000000ull + ts.tv_nsec;
-#else
+#elif !defined(OS_EMBEDDED)
   struct timeval tv;
   gettimeofday(&tv,NULL);
   return (unsigned long long)tv.tv_sec * 1000000000ull + tv.tv_usec * 1000;
+#else
+  return 0;
 #endif
 }
 #define RPCC_DEFINED
@@ -521,6 +526,10 @@ static void __inline blas_lock(volatile BLASULONG *address){
 #include "common_linux.h"
 #endif
 
+#ifdef OS_EMBEDDED
+#define DTB_DEFAULT_ENTRIES 64
+#endif
+
 #define MMAP_ACCESS (PROT_READ | PROT_WRITE)
 
 #ifdef __NetBSD__

common_arm64.h

@@ -39,7 +39,7 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #define INLINE inline
 
-#ifdef F_INTERFACE_FLANG
+#if defined( F_INTERFACE_FLANG) || defined(F_INTERFACE_PGI)
 #define RETURN_BY_STACK
 #else
 #define RETURN_BY_COMPLEX

cpuid_x86.c

@@ -1418,6 +1418,15 @@ int get_cpuname(void){
       case 9:
       case 8:      
         switch (model) {
+        case 12: // Tiger Lake
+          if(support_avx512())
+            return CPUTYPE_SKYLAKEX;
+          if(support_avx2())
+            return CPUTYPE_HASWELL;
+          if(support_avx())
+            return CPUTYPE_SANDYBRIDGE;
+          else
+          return CPUTYPE_NEHALEM;
 	case 14: // Kaby Lake and refreshes
           if(support_avx2())
             return CPUTYPE_HASWELL;
@@ -1436,6 +1445,15 @@ int get_cpuname(void){
         return CPUTYPE_SANDYBRIDGE;
           else
         return CPUTYPE_NEHALEM;
+    case 7: // Rocket Lake           
+	  if(support_avx512())
+            return CPUTYPE_SKYLAKEX;
+          if(support_avx2())
+            return CPUTYPE_HASWELL;
+          if(support_avx())
+	    return CPUTYPE_SANDYBRIDGE;
+	  else
+	  return CPUTYPE_NEHALEM;
 	}
 	break;    
       }
@@ -2014,6 +2032,19 @@ int get_coretype(void){
   #endif
             else
               return CORE_NEHALEM;
+	  case 7:// Rocket Lake
+#ifndef NO_AVX512
+	  if(support_avx512())
+            return CORE_SKYLAKEX;
+#endif
+#ifndef NO_AVX2
+	  if(support_avx2())
+            return CORE_HASWELL;
+#endif
+	  if(support_avx())
+	    return CORE_SANDYBRIDGE;
+	  else
+	  return CORE_NEHALEM;
         }
       case 5:
         switch (model) {
@@ -2102,6 +2133,16 @@ int get_coretype(void){
         break;    	
       case 9:
       case 8:
+       if (model == 12) { // Tiger Lake
+          if(support_avx512())
+            return CPUTYPE_SKYLAKEX;
+          if(support_avx2())
+            return CPUTYPE_HASWELL;
+          if(support_avx())
+            return CPUTYPE_SANDYBRIDGE;
+          else
+          return CPUTYPE_NEHALEM;
+        }
         if (model == 14) { // Kaby Lake 
 	  if(support_avx())
 #ifndef NO_AVX2

driver/others/blas_server.c

@@ -1024,38 +1024,39 @@ int BLASFUNC(blas_thread_shutdown)(void){
 
   int i;
 
-  if (!blas_server_avail) return 0;
-
   LOCK_COMMAND(&server_lock);
 
-  for (i = 0; i < blas_num_threads - 1; i++) {
+  if (blas_server_avail) {
+
+    for (i = 0; i < blas_num_threads - 1; i++) {
 
 
-    pthread_mutex_lock (&thread_status[i].lock);
+      pthread_mutex_lock (&thread_status[i].lock);
 
-    atomic_store_queue(&thread_status[i].queue, (blas_queue_t *)-1);
-    thread_status[i].status = THREAD_STATUS_WAKEUP;
-    pthread_cond_signal (&thread_status[i].wakeup);
+      atomic_store_queue(&thread_status[i].queue, (blas_queue_t *)-1);
+      thread_status[i].status = THREAD_STATUS_WAKEUP;
+      pthread_cond_signal (&thread_status[i].wakeup);
 
-    pthread_mutex_unlock(&thread_status[i].lock);
+      pthread_mutex_unlock(&thread_status[i].lock);
 
-  }
+    }
 
-  for(i = 0; i < blas_num_threads - 1; i++){
-    pthread_join(blas_threads[i], NULL);
-  }
+    for(i = 0; i < blas_num_threads - 1; i++){
+      pthread_join(blas_threads[i], NULL);
+    }
 
-  for(i = 0; i < blas_num_threads - 1; i++){
-    pthread_mutex_destroy(&thread_status[i].lock);
-    pthread_cond_destroy (&thread_status[i].wakeup);
-  }
+    for(i = 0; i < blas_num_threads - 1; i++){
+      pthread_mutex_destroy(&thread_status[i].lock);
+      pthread_cond_destroy (&thread_status[i].wakeup);
+    }
 
 #ifdef NEED_STACKATTR
-  pthread_attr_destory(&attr);
+    pthread_attr_destroy(&attr);
 #endif
 
-  blas_server_avail = 0;
+    blas_server_avail = 0;
 
+  }
   UNLOCK_COMMAND(&server_lock);
 
   return 0;

driver/others/dynamic.c

@@ -644,6 +644,21 @@ static gotoblas_t *get_coretype(void){
         return NULL;  
       case 9:
       case 8:
+        if (model == 12) { // Tiger Lake
+          if (support_avx512()) 
+            return &gotoblas_SKYLAKEX;
+          if(support_avx2()){
+            openblas_warning(FALLBACK_VERBOSE, HASWELL_FALLBACK);
+            return &gotoblas_HASWELL;
+          }
+          if(support_avx()) {
+            openblas_warning(FALLBACK_VERBOSE, SANDYBRIDGE_FALLBACK);
+            return &gotoblas_SANDYBRIDGE;
+          } else {
+          openblas_warning(FALLBACK_VERBOSE, NEHALEM_FALLBACK);
+          return &gotoblas_NEHALEM;
+          }
+        }
 	if (model == 14 ) { // Kaby Lake, Coffee Lake
 	  if(support_avx2())
 	    return &gotoblas_HASWELL;
@@ -656,7 +671,7 @@ static gotoblas_t *get_coretype(void){
 	  }
 	}
       case 10:
-    if (model == 5 || model == 6) {
+        if (model == 5 || model == 6) {
 	  if(support_avx2())
 	    return &gotoblas_HASWELL;
 	  if(support_avx()) {
@@ -666,7 +681,20 @@ static gotoblas_t *get_coretype(void){
 	    openblas_warning(FALLBACK_VERBOSE, NEHALEM_FALLBACK);
 	    return &gotoblas_NEHALEM; //OS doesn't support AVX. Use old kernels.
 	  }
-    }
+        }
+        if (model == 7) {
+	  if (support_avx512()) 
+	    return &gotoblas_SKYLAKEX;
+	  if(support_avx2())
+	    return &gotoblas_HASWELL;
+	  if(support_avx()) {
+	    openblas_warning(FALLBACK_VERBOSE, SANDYBRIDGE_FALLBACK);
+	    return &gotoblas_SANDYBRIDGE;
+	  } else {
+	    openblas_warning(FALLBACK_VERBOSE, NEHALEM_FALLBACK);
+	    return &gotoblas_NEHALEM; //OS doesn't support AVX. Use old kernels.
+	  }
+        }      
 	return NULL;
       }
       case 0xf:

driver/others/dynamic_arm64.c

@@ -43,6 +43,63 @@
 #endif
 
 extern gotoblas_t  gotoblas_ARMV8;
+#ifdef DYNAMIC_LIST
+#ifdef DYN_CORTEXA53
+extern gotoblas_t  gotoblas_CORTEXA53;
+#else
+#define gotoblas_CORTEXA53 gotoblas_ARMV8
+#endif
+#ifdef DYN_CORTEXA57
+extern gotoblas_t  gotoblas_CORTEXA57;
+#else
+#define gotoblas_CORTEXA57 gotoblas_ARMV8
+#endif
+#ifdef DYN_CORTEXA72
+extern gotoblas_t  gotoblas_CORTEXA72;
+#else
+#define gotoblas_CORTEXA72 gotoblas_ARMV8
+#endif
+#ifdef DYN_CORTEXA73
+extern gotoblas_t  gotoblas_CORTEXA73;
+#else
+#define gotoblas_CORTEXA73 gotoblas_ARMV8
+#endif
+#ifdef DYN_FALKOR
+extern gotoblas_t  gotoblas_FALKOR;
+#else
+#define gotoblas_FALKOR gotoblas_ARMV8
+#endif
+#ifdef DYN_TSV110
+extern gotoblas_t  gotoblas_TSV110;
+#else
+#define gotoblas_TSV110 gotoblas_ARMV8
+#endif
+#ifdef DYN_THUNDERX
+extern gotoblas_t  gotoblas_THUNDERX;
+#else
+#define gotoblas_THUNDERX gotoblas_ARMV8
+#endif
+#ifdef DYN_THUNDERX2T99
+extern gotoblas_t  gotoblas_THUNDERX2T99;
+#else
+#define gotoblas_THUNDERX2T99 gotoblas_ARMV8
+#endif
+#ifdef DYN_THUNDERX3T110
+extern gotoblas_t  gotoblas_THUNDERX3T110;
+#else
+#define gotoblas_THUNDERX3T110 gotoblas_ARMV8
+#endif
+#ifdef DYN_EMAG8180
+extern gotoblas_t  gotoblas_EMAG8180;
+#else
+#define gotoblas_EMAG8180 gotoblas_ARMV8
+#endif
+#ifdef DYN_NEOVERSEN1
+extern gotoblas_t  gotoblas_NEOVERSEN1;
+#else
+#define gotoblas_NEOVERSEN1 gotoblas_ARMV8
+#endif
+#else
 extern gotoblas_t  gotoblas_CORTEXA53;
 extern gotoblas_t  gotoblas_CORTEXA57;
 extern gotoblas_t  gotoblas_CORTEXA72;
@@ -54,6 +111,7 @@ extern gotoblas_t  gotoblas_TSV110;
 extern gotoblas_t  gotoblas_EMAG8180;
 extern gotoblas_t  gotoblas_NEOVERSEN1;
 extern gotoblas_t  gotoblas_THUNDERX3T110;
+#endif
 
 extern void openblas_warning(int verbose, const char * msg);
 
@@ -68,7 +126,7 @@ extern void openblas_warning(int verbose, const char * msg);
 #endif
 
 #define get_cpu_ftr(id, var) ({					\
-		__asm__("mrs %0, "#id : "=r" (var));		\
+		__asm__ ("mrs %0, "#id : "=r" (var));		\
 	})
 
 static char *corename[] = {

driver/others/dynamic_power.c

@@ -202,7 +202,7 @@ static gotoblas_t *get_coretype(void) {
 		return &gotoblas_POWER10;
 #endif
 	/* Fall back to the POWER9 implementation if the toolchain is too old or the MMA feature is not set */
-#if (!defined __GNUC__) || ( __GNUC__ >= 6)
+#if (!defined __GNUC__) || ( __GNUC__ >= 11) || (__GNUC__ == 10 && __GNUC_MINOR__ >= 2)
 	if (__builtin_cpu_is("power10"))
 		return &gotoblas_POWER9;
 #endif	

driver/others/memory.c

@@ -222,11 +222,11 @@ int get_num_procs(void);
 #else
 int get_num_procs(void) {
   static int nums = 0;
+
+#if defined(__GLIBC_PREREQ)
   cpu_set_t cpuset,*cpusetp;
   size_t size;
   int ret;
-
-#if defined(__GLIBC_PREREQ)
 #if !__GLIBC_PREREQ(2, 7)
   int i;
 #if !__GLIBC_PREREQ(2, 6)
@@ -1241,7 +1241,7 @@ UNLOCK_COMMAND(&alloc_lock);
 
       func = &memoryalloc[0];
 
-      while ((func != NULL) && (map_address == (void *) -1)) {
+      while ((*func != NULL) && (map_address == (void *) -1)) {
 
         map_address = (*func)((void *)base_address);
 
@@ -1619,10 +1619,12 @@ static int on_process_term(void)
 #else
 #pragma data_seg(".CRT$XLB")
 #endif
-static void (APIENTRY *dll_callback)(HINSTANCE h, DWORD ul_reason_for_call, PVOID pv) = DllMain;
+
 #ifdef _WIN64
+static const PIMAGE_TLS_CALLBACK dll_callback(HINSTANCE h, DWORD ul_reason_for_call, PVOID pv) = DllMain;
 #pragma const_seg()
 #else
+static void (APIENTRY *dll_callback)(HINSTANCE h, DWORD ul_reason_for_call, PVOID pv) = DllMain;
 #pragma data_seg()
 #endif
 
@@ -1631,10 +1633,12 @@ static void (APIENTRY *dll_callback)(HINSTANCE h, DWORD ul_reason_for_call, PVOI
 #else
 #pragma data_seg(".CRT$XTU")
 #endif
-static int(*p_process_term)(void) = on_process_term;
+
 #ifdef _WIN64
+static const int(*p_process_term)(void) = on_process_term;
 #pragma const_seg()
 #else
+static int(*p_process_term)(void) = on_process_term;
 #pragma data_seg()
 #endif
 #endif
@@ -1668,16 +1672,23 @@ void gotoblas_dummy_for_PGI(void) {
 #ifndef MEM_LARGE_PAGES
 #define MEM_LARGE_PAGES  0x20000000
 #endif
-#else
+#elif !defined(OS_EMBEDDED)
 #define ALLOC_MMAP
 #define ALLOC_MALLOC
+#else
+#define ALLOC_MALLOC
+
+inline int puts(const char *str) { return 0; }
+inline int printf(const char *format, ...) { return 0; }
+inline char *getenv(const char *name) { return ""; }
+inline int atoi(const char *str) { return 0; }
 #endif
 
 #include <stdlib.h>
 #include <stdio.h>
 #include <fcntl.h>
 
-#if !defined(OS_WINDOWS) || defined(OS_CYGWIN_NT)
+#if (!defined(OS_WINDOWS) || defined(OS_CYGWIN_NT)) && !defined(OS_EMBEDDED)
 #include <sys/mman.h>
 #ifndef NO_SYSV_IPC
 #include <sys/shm.h>

exports/gensymbol

@@ -1,4 +1,4 @@
-#!/usr/bin/perl
+#!/usr/bin/env perl
 
 # Changelog
 # 2017/09/03 staticfloat

f_check

@@ -1,4 +1,4 @@
-#!/usr/bin/perl
+#!/usr/bin/env perl
 
 $hostos   = `uname -s | sed -e s/\-.*//`;    chop($hostos);
 
@@ -75,6 +75,7 @@ if ($compiler eq "") {
 
 	} elsif ($data =~ /GNU/ || $data =~ /GCC/ ) {
 
+            $data =~ s/\(+.*?\)+//g;
 	    $data =~ /(\d+)\.(\d+).(\d+)/;
 	    $major = $1;
 	    $minor = $2;

interface/Makefile

@@ -316,7 +316,7 @@ CCBLAS1OBJS   = \
 	cblas_cscal.$(SUFFIX) cblas_csscal.$(SUFFIX) \
 	cblas_cswap.$(SUFFIX) cblas_scnrm2.$(SUFFIX) \
 	cblas_caxpby.$(SUFFIX) \
-	cblas_icmin.$(SUFFIX) cblas_icmax.$(SUFFIX) cblas_scsum.$(SUFFIX)
+	cblas_icmin.$(SUFFIX) cblas_icmax.$(SUFFIX) cblas_scsum.$(SUFFIX) cblas_csrot.$(SUFFIX) cblas_crotg.$(SUFFIX)
 
 CCBLAS2OBJS   = \
 	cblas_cgemv.$(SUFFIX) cblas_cgerc.$(SUFFIX) cblas_cgeru.$(SUFFIX) \
@@ -346,7 +346,7 @@ CZBLAS1OBJS   = \
 	cblas_zscal.$(SUFFIX) cblas_zdscal.$(SUFFIX) \
 	cblas_zswap.$(SUFFIX) cblas_dznrm2.$(SUFFIX) \
 	cblas_zaxpby.$(SUFFIX) \
-	cblas_izmin.$(SUFFIX) cblas_izmax.$(SUFFIX) cblas_dzsum.$(SUFFIX)
+	cblas_izmin.$(SUFFIX) cblas_izmax.$(SUFFIX) cblas_dzsum.$(SUFFIX) cblas_zdrot.$(SUFFIX) cblas_zrotg.$(SUFFIX)
 
 
 CZBLAS2OBJS   = \
@@ -1634,6 +1634,12 @@ cblas_srotg.$(SUFFIX) cblas_srotg.$(PSUFFIX): rotg.c
 cblas_drotg.$(SUFFIX) cblas_drotg.$(PSUFFIX): rotg.c
 	$(CC) $(CFLAGS) -DCBLAS -c $< -o $(@F)
 
+cblas_crotg.$(SUFFIX) cblas_crotg.$(PSUFFIX): zrotg.c
+	$(CC) -c $(CFLAGS) -DCBLAS $< -o $(@F)
+
+cblas_zrotg.$(SUFFIX) cblas_zrotg.$(PSUFFIX): zrotg.c
+	$(CC) -c $(CFLAGS) -DCBLAS $< -o $(@F)
+
 cblas_srotm.$(SUFFIX) cblas_srotm.$(PSUFFIX): rotm.c
 	$(CC) $(CFLAGS) -DCBLAS -c $< -o $(@F)
 
@@ -1664,6 +1670,12 @@ cblas_csscal.$(SUFFIX) cblas_csscal.$(PSUFFIX) : zscal.c
 cblas_zdscal.$(SUFFIX) cblas_zdscal.$(PSUFFIX) : zscal.c
 	$(CC) $(CFLAGS) -DCBLAS -c -DSSCAL $< -o $(@F)
 
+cblas_csrot.$(SUFFIX) cblas_csrot.$(PSUFFIX) : zrot.c
+	$(CC) $(CFLAGS) -DCBLAS -c $< -o $(@F)
+
+cblas_zdrot.$(SUFFIX) cblas_zdrot.$(PSUFFIX) : zrot.c
+	$(CC) $(CFLAGS) -DCBLAS -c $< -o $(@F)
+
 ifeq ($(BUILD_BFLOAT16),1)
 cblas_sbgemv.$(SUFFIX) cblas_sbgemv.$(PSUFFIX) : sbgemv.c
 	$(CC) -DCBLAS -c $(CFLAGS) $< -o $(@F)

interface/create

@@ -1,4 +1,4 @@
-#!/usr/bin/perl
+#!/usr/bin/env perl
 
 $count = 0;
 

interface/gemm.c

@@ -246,6 +246,7 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANS
 
 #ifdef SMP
   double MNK;
+#if defined(USE_SIMPLE_THREADED_LEVEL3) || !defined(NO_AFFINITY)
 #ifndef COMPLEX
 #ifdef XDOUBLE
   int mode  =  BLAS_XDOUBLE | BLAS_REAL;
@@ -264,6 +265,7 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANS
 #endif
 #endif
 #endif
+#endif
 
 #if defined(SMP) && !defined(NO_AFFINITY) && !defined(USE_SIMPLE_THREADED_LEVEL3)
   int nodes;
@@ -417,8 +419,10 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANS
   sb = (XFLOAT *)(((BLASLONG)sa + ((GEMM_P * GEMM_Q * COMPSIZE * SIZE + GEMM_ALIGN) & ~GEMM_ALIGN)) + GEMM_OFFSET_B);
 
 #ifdef SMP
+#if defined(USE_SIMPLE_THREADED_LEVEL3) || !defined(NO_AFFINITY)
   mode |= (transa << BLAS_TRANSA_SHIFT);
   mode |= (transb << BLAS_TRANSB_SHIFT);
+#endif
 
   MNK = (double) args.m * (double) args.n * (double) args.k;
   if ( MNK <= (SMP_THRESHOLD_MIN  * (double) GEMM_MULTITHREAD_THRESHOLD)  )

interface/rotmg.c

@@ -107,7 +107,6 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
 		dq1 =  dp1 * *dx1;
 		if(ABS(dq1) > ABS(dq2))
 		{
-			dflag = ZERO;
 			dh11  =  ONE;
 			dh22  =  ONE;
 			dh21 = -  dy1 / *dx1;

interface/zhbmv.c

@@ -61,7 +61,7 @@ static  int (*hbmv[])(BLASLONG, BLASLONG, FLOAT, FLOAT, FLOAT *, BLASLONG, FLOAT
 #endif
 };
 
-#ifdef SMP
+#ifdef SMPBUG
 static  int (*hbmv_thread[])(BLASLONG, BLASLONG, FLOAT *, FLOAT *, BLASLONG, FLOAT *, BLASLONG, FLOAT *, BLASLONG, FLOAT *, int) = {
 #ifdef XDOUBLE
   xhbmv_thread_U, xhbmv_thread_L, xhbmv_thread_V, xhbmv_thread_M,
@@ -92,7 +92,7 @@ void NAME(char *UPLO, blasint *N, blasint *K, FLOAT  *ALPHA, FLOAT *a, blasint *
   blasint info;
   int uplo;
   FLOAT *buffer;
-#ifdef SMP
+#ifdef SMPBUG
   int nthreads;
 #endif
 
@@ -144,7 +144,7 @@ void CNAME(enum CBLAS_ORDER order,
   FLOAT *buffer;
   int uplo;
   blasint info;
-#ifdef SMP
+#ifdef SMPBUG
   int nthreads;
 #endif
 
@@ -203,7 +203,7 @@ void CNAME(enum CBLAS_ORDER order,
 
   buffer = (FLOAT *)blas_memory_alloc(1);
 
-#ifdef SMP
+#ifdef SMPBUG
   nthreads = num_cpu_avail(2);
 
   if (nthreads == 1) {
@@ -211,7 +211,7 @@ void CNAME(enum CBLAS_ORDER order,
 
     (hbmv[uplo])(n, k, alpha_r, alpha_i, a, lda, x, incx, y, incy, buffer);
 
-#ifdef SMP
+#ifdef SMPBUG
   } else {
 
     (hbmv_thread[uplo])(n, k, ALPHA, a, lda, x, incx, y, incy, buffer, nthreads);

kernel/CMakeLists.txt

@@ -187,10 +187,11 @@ function (build_core TARGET_CORE KDIR TSUFFIX KERNEL_DEFINITIONS)
     endif ()
     # Makefile.L3
     set(USE_TRMM false)
-    if (ARM OR ARM64 OR (TARGET_CORE MATCHES LONGSOON3B) OR (TARGET_CORE MATCHES GENERIC) OR (TARGET_CORE MATCHES HASWELL) OR (TARGET_CORE MATCHES ZEN) OR (TARGET_CORE MATCHES SKYLAKEX) OR (TARGET_CORE MATCHES COOPERLAKE))
+    string(TOUPPER ${TARGET_CORE} UC_TARGET_CORE)
+    if (ARM OR ARM64 OR (UC_TARGET_CORE MATCHES LONGSOON3B) OR (UC_TARGET_CORE MATCHES GENERIC) OR (UC_TARGET_CORE MATCHES HASWELL) OR (UC_TARGET_CORE MATCHES ZEN) OR (UC_TARGET_CORE MATCHES SKYLAKEX) OR (UC_TARGET_CORE MATCHES COOPERLAKE))
       set(USE_TRMM true)
     endif ()
-    if (ZARCH OR (TARGET_CORE MATCHES POWER8) OR (TARGET_CORE MATCHES POWER9) OR (TARGET_CORE MATCHES POWER10))
+    if (ZARCH OR (UC_TARGET_CORE MATCHES POWER8) OR (UC_TARGET_CORE MATCHES POWER9) OR (UC_TARGET_CORE MATCHES POWER10))
       set(USE_TRMM true)
     endif ()
 

kernel/arm/omatcopy_rt.c

@@ -1,5 +1,5 @@
 /***************************************************************************
-Copyright (c) 2013, The OpenBLAS Project
+Copyright (c) 2021, 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
@@ -27,36 +27,208 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "common.h"
 
-/*****************************************************
- * 2014/06/09 Saar
- *
- * Order rowMajor
- * Trans
- *
-******************************************************/
-
 int CNAME(BLASLONG rows, BLASLONG cols, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *b, BLASLONG ldb)
 {
-	BLASLONG i,j;
-	FLOAT *aptr,*bptr;
+    BLASLONG i, j;
+    FLOAT *a_offset, *a_offset1, *a_offset2, *a_offset3, *a_offset4;
+    FLOAT *b_offset, *b_offset1, *b_offset2, *b_offset3, *b_offset4;
 
-	if ( rows <= 0     )  return(0);
-	if ( cols <= 0     )  return(0);
+    if (rows <= 0)  return 0;
+    if (cols <= 0)  return 0;
 
-	aptr = a;
+    a_offset = a;
+    b_offset = b;
 
-	for ( i=0; i<rows ; i++ )
-	{
-		bptr = &b[i];
-		for(j=0; j<cols; j++)
-		{
-			bptr[j*ldb] = alpha * aptr[j];
-		}
-		aptr += lda;
-	}
+    i = (rows >> 2);
+    if (i > 0) {
+        do {
+            a_offset1 = a_offset;
+            a_offset2 = a_offset1 + lda;
+            a_offset3 = a_offset2 + lda;
+            a_offset4 = a_offset3 + lda;
+            a_offset += 4 * lda;
 
-	return(0);
+            b_offset1 = b_offset;
+            b_offset2 = b_offset1 + ldb;
+            b_offset3 = b_offset2 + ldb;
+            b_offset4 = b_offset3 + ldb;
+            b_offset += 4;
 
+            j = (cols >> 2);
+            if (j > 0) {
+                do {
+                /* Column 1 of MAT_B */
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha; // Row 1 of MAT_A
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+                *(b_offset3 + 0) = *(a_offset1 + 2)*alpha;
+                *(b_offset4 + 0) = *(a_offset1 + 3)*alpha;
+
+                /* Column 2 of MAT_B */
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha; // Row 2 of MAT_A
+                *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+                *(b_offset3 + 1) = *(a_offset2 + 2)*alpha;
+                *(b_offset4 + 1) = *(a_offset2 + 3)*alpha;
+
+                /* Column 3 of MAT_B */
+                *(b_offset1 + 2) = *(a_offset3 + 0)*alpha; // Row 3 of MAT_A
+                *(b_offset2 + 2) = *(a_offset3 + 1)*alpha;
+                *(b_offset3 + 2) = *(a_offset3 + 2)*alpha;
+                *(b_offset4 + 2) = *(a_offset3 + 3)*alpha;
+
+                /* Column 4 of MAT_B */
+                *(b_offset1 + 3) = *(a_offset4 + 0)*alpha; // Row 4 of MAT_A
+                *(b_offset2 + 3) = *(a_offset4 + 1)*alpha;
+                *(b_offset3 + 3) = *(a_offset4 + 2)*alpha;
+                *(b_offset4 + 3) = *(a_offset4 + 3)*alpha;
+
+                a_offset1 += 4;
+                a_offset2 += 4;
+                a_offset3 += 4;
+                a_offset4 += 4;
+                b_offset1 += ldb * 4;
+                b_offset2 += ldb * 4;
+                b_offset3 += ldb * 4;
+                b_offset4 += ldb * 4;
+                
+                j--;
+                } while (j > 0);
+            } // if(j > 0)
+
+
+            if (cols & 2) {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+                *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+
+                *(b_offset1 + 2) = *(a_offset3 + 0)*alpha;
+                *(b_offset2 + 2) = *(a_offset3 + 1)*alpha;
+
+                *(b_offset1 + 3) = *(a_offset4 + 0)*alpha;
+                *(b_offset2 + 3) = *(a_offset4 + 1)*alpha;
+
+                a_offset1 += 2;
+                a_offset2 += 2;
+                a_offset3 += 2;
+                a_offset4 += 2;
+
+                b_offset1 += ldb*2;
+
+            }
+
+            if (cols & 1) {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+
+                *(b_offset1 + 2) = *(a_offset3 + 0)*alpha;
+
+                *(b_offset1 + 3) = *(a_offset4 + 0)*alpha;
+            }
+
+            i--;
+        } while (i > 0);
+    }
+
+
+    if (rows & 2) {
+        a_offset1 = a_offset;
+        a_offset2 = a_offset1 + lda;
+        a_offset += 2 * lda;
+
+        b_offset1 = b_offset;
+        b_offset2 = b_offset1 + ldb;
+        b_offset3 = b_offset2 + ldb;
+        b_offset4 = b_offset3 + ldb;
+        b_offset += 2;
+
+        j = (cols >> 2);
+        if (j > 0){
+            do {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+                *(b_offset3 + 0) = *(a_offset1 + 2)*alpha;
+                *(b_offset4 + 0) = *(a_offset1 + 3)*alpha;
+
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+                *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+                *(b_offset3 + 1) = *(a_offset2 + 2)*alpha;
+                *(b_offset4 + 1) = *(a_offset2 + 3)*alpha;
+                
+                a_offset1 += 4;
+                a_offset2 += 4;
+                b_offset1 += ldb * 4;
+                b_offset2 += ldb * 4;
+                b_offset3 += ldb * 4;
+                b_offset4 += ldb * 4;
+
+                j--;
+            } while (j > 0);
+        }
+
+
+        if (cols & 2){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+            *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+
+            *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+            *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+
+            a_offset1 += 2;
+            a_offset2 += 2;
+            b_offset1 += ldb*2;
+
+        }
+
+
+        if (cols & 1){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+            *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+        }
+    } // if (rows & 2)
+
+
+    if (rows & 1) {
+        a_offset1 = a_offset;
+        a_offset += lda;
+
+        b_offset1 = b_offset;
+        b_offset2 = b_offset1 + ldb;
+        b_offset3 = b_offset2 + ldb;
+        b_offset4 = b_offset3 + ldb;
+
+        j = (cols >> 2);
+        if (j > 0){
+            do {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+                *(b_offset3 + 0) = *(a_offset1 + 2)*alpha;
+                *(b_offset4 + 0) = *(a_offset1 + 3)*alpha;
+                
+                a_offset1 += 4;
+                b_offset1 += ldb * 4;
+                b_offset2 += ldb * 4;
+                b_offset3 += ldb * 4;
+                b_offset4 += ldb * 4;
+                
+                j--;
+            } while (j > 0);
+        }
+
+        if (cols & 2){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+            *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+
+            a_offset1 += 2;
+            b_offset1 += ldb * 2;
+        }
+
+        if (cols & 1){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+        }
+    }
+
+    return 0;
 }
 
-

kernel/arm/zdot.c

@@ -48,7 +48,7 @@ OPENBLAS_COMPLEX_FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLA
 
 	dot[0]=0.0;
 	dot[1]=0.0;
-#if !defined(__PPC__) && !defined(__SunOS)
+#if !defined(__PPC__) && !defined(__SunOS) && !defined(__PGI)
 	CREAL(result) = 0.0 ;
 	CIMAG(result) = 0.0 ;
 #else
@@ -73,7 +73,7 @@ OPENBLAS_COMPLEX_FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLA
 		i++ ;
 
 	}
-#if !defined(__PPC__)	&& !defined(__SunOS)
+#if !defined(__PPC__)	&& !defined(__SunOS) && !defined(__PGI)
         CREAL(result) = dot[0];
 	CIMAG(result) = dot[1];
 #else

kernel/arm64/KERNEL.ARMV8

@@ -97,9 +97,18 @@ CNRM2KERNEL    = znrm2.S
 ZNRM2KERNEL    = znrm2.S
 
 DDOTKERNEL     = dot.S
+ifneq ($(C_COMPILER), PGI)
 SDOTKERNEL     = ../generic/dot.c
+else
+SDOTKERNEL = dot.S
+endif
+ifneq ($(C_COMPILER), PGI)
 CDOTKERNEL     = zdot.S
 ZDOTKERNEL     = zdot.S
+else
+CDOTKERNEL = ../arm/zdot.c
+ZDOTKERNEL = ../arm/zdot.c
+endif
 DSDOTKERNEL    = dot.S
 
 DGEMM_BETA     = dgemm_beta.S

kernel/arm64/KERNEL.CORTEXA53

@@ -96,11 +96,20 @@ DNRM2KERNEL    = nrm2.S
 CNRM2KERNEL    = znrm2.S
 ZNRM2KERNEL    = znrm2.S
 
-DDOTKERNEL     = dot.S
-SDOTKERNEL     = ../generic/dot.c
-CDOTKERNEL     = zdot.S
-ZDOTKERNEL     = zdot.S
-DSDOTKERNEL    = dot.S
+ifneq ($(C_COMPILER), PGI)
+SDOTKERNEL   = ../generic/dot.c
+else
+SDOTKERNEL   = dot.S
+endif
+DDOTKERNEL   = dot.S
+ifneq ($(C_COMPILER), PGI)
+CDOTKERNEL   = zdot.S
+ZDOTKERNEL   = zdot.S
+else
+CDOTKERNEL = ../arm/zdot.c
+ZDOTKERNEL = ../arm/zdot.c
+endif
+DSDOTKERNEL  = dot.S
 
 DGEMM_BETA     = dgemm_beta.S
 SGEMM_BETA     = sgemm_beta.S

kernel/arm64/KERNEL.CORTEXA57

@@ -70,10 +70,19 @@ DCOPYKERNEL  = copy.S
 CCOPYKERNEL  = copy.S
 ZCOPYKERNEL  = copy.S
 
+ifneq ($(C_COMPILER), PGI)
 SDOTKERNEL   = ../generic/dot.c
+else
+SDOTKERNEL   = dot.S
+endif
 DDOTKERNEL   = dot.S
+ifneq ($(C_COMPILER), PGI)
 CDOTKERNEL   = zdot.S
 ZDOTKERNEL   = zdot.S
+else
+CDOTKERNEL = ../arm/zdot.c
+ZDOTKERNEL = ../arm/zdot.c
+endif
 DSDOTKERNEL  = dot.S
 
 SNRM2KERNEL  = nrm2.S

kernel/arm64/KERNEL.NEOVERSEN1

@@ -91,10 +91,10 @@ IDAMAXKERNEL   = iamax_thunderx2t99.c
 ICAMAXKERNEL   = izamax_thunderx2t99.c
 IZAMAXKERNEL   = izamax_thunderx2t99.c
 
-SNRM2KERNEL    = nrm2.S
-DNRM2KERNEL    = nrm2.S
-CNRM2KERNEL    = znrm2.S
-ZNRM2KERNEL    = znrm2.S
+SNRM2KERNEL    = scnrm2_thunderx2t99.c
+DNRM2KERNEL    = dznrm2_thunderx2t99.c
+CNRM2KERNEL    = scnrm2_thunderx2t99.c
+ZNRM2KERNEL    = dznrm2_thunderx2t99.c
 
 DDOTKERNEL     = dot_thunderx2t99.c
 SDOTKERNEL     = dot_thunderx2t99.c

kernel/arm64/KERNEL.THUNDERX

@@ -47,8 +47,13 @@ ZCOPYKERNEL  = copy.S
 
 SDOTKERNEL   = dot_thunderx.c
 DDOTKERNEL   = ddot_thunderx.c
+ifneq ($(C_COMPILER), PGI)
 CDOTKERNEL   = zdot.S
 ZDOTKERNEL   = zdot.S
+else
+CDOTKERNEL = ../arm/zdot.c
+ZDOTKERNEL = ../arm/zdot.c
+endif
 DSDOTKERNEL  = dot.S
 
 SNRM2KERNEL  = nrm2.S

kernel/arm64/KERNEL.THUNDERX2T99

@@ -153,12 +153,12 @@ IDAMAXKERNEL   = iamax_thunderx2t99.c
 ICAMAXKERNEL   = izamax_thunderx2t99.c
 IZAMAXKERNEL   = izamax_thunderx2t99.c
 
-SNRM2KERNEL    = nrm2.S
-CNRM2KERNEL    = nrm2.S
+SNRM2KERNEL    = scnrm2_thunderx2t99.c
+CNRM2KERNEL    = scnrm2_thunderx2t99.c
 #DNRM2KERNEL    = dznrm2_thunderx2t99_fast.c
 #ZNRM2KERNEL    = dznrm2_thunderx2t99_fast.c
-DNRM2KERNEL    = znrm2.S
-ZNRM2KERNEL    = znrm2.S
+DNRM2KERNEL    = dznrm2_thunderx2t99.c
+ZNRM2KERNEL    = dznrm2_thunderx2t99.c
 
 
 DDOTKERNEL     = dot_thunderx2t99.c

kernel/arm64/KERNEL.THUNDERX3T110

@@ -153,16 +153,13 @@ IDAMAXKERNEL   = iamax_thunderx2t99.c
 ICAMAXKERNEL   = izamax_thunderx2t99.c
 IZAMAXKERNEL   = izamax_thunderx2t99.c
 
-#SNRM2KERNEL    = scnrm2_thunderx2t99.c
-#CNRM2KERNEL    = scnrm2_thunderx2t99.c
-##DNRM2KERNEL    = dznrm2_thunderx2t99_fast.c
-##ZNRM2KERNEL    = dznrm2_thunderx2t99_fast.c
-#DNRM2KERNEL    = dznrm2_thunderx2t99.c
-#ZNRM2KERNEL    = dznrm2_thunderx2t99.c
-SNRM2KERNEL = nrm2.S
-DNRM2KERNEL = nrm2.S
-CNRM2KERNEL = znrm2.S
-ZNRM2KERNEL = znrm2.S
+SNRM2KERNEL    = scnrm2_thunderx2t99.c
+CNRM2KERNEL    = scnrm2_thunderx2t99.c
+#DNRM2KERNEL    = dznrm2_thunderx2t99_fast.c
+#ZNRM2KERNEL    = dznrm2_thunderx2t99_fast.c
+DNRM2KERNEL    = dznrm2_thunderx2t99.c
+ZNRM2KERNEL    = dznrm2_thunderx2t99.c
+
 
 DDOTKERNEL     = dot_thunderx2t99.c
 SDOTKERNEL     = dot_thunderx2t99.c

kernel/arm64/KERNEL.TSV110

@@ -72,8 +72,13 @@ ZCOPYKERNEL  = copy.S
 
 SDOTKERNEL   = dot.S
 DDOTKERNEL   = dot.S
+ifneq ($(C_COMPILER), PGI)
 CDOTKERNEL   = zdot.S
 ZDOTKERNEL   = zdot.S
+else
+CDOTKERNEL = ../arm/zdot.c
+ZDOTKERNEL = ../arm/zdot.c
+endif
 DSDOTKERNEL  = dot.S
 
 SNRM2KERNEL  = nrm2.S

kernel/arm64/dznrm2_thunderx2t99.c

@@ -58,6 +58,7 @@ extern int blas_level1_thread_with_return_value(int mode, BLASLONG m, BLASLONG n
 #define CUR_MAXINV	"d8"
 #define CUR_MAXINV_V	"v8.2d"
 #define CUR_MAX_V	"v8.2d"
+#define REGINF		"d9"
 
 static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 		         double *ssq, double *scale)
@@ -79,8 +80,10 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"	ble	9f //nrm2_kernel_L999			\n"
 
 	"1: //nrm2_kernel_F_BEGIN:				\n"
+	"	mov	x6, #0x7FF0000000000000 //+Infinity	\n"
 	"	fmov	"REGZERO", xzr				\n"
 	"	fmov	"REGONE", #1.0				\n"
+	"	fmov	"REGINF", x6				\n"
 	"	lsl	"INC_X", "INC_X", #"INC_SHIFT"		\n"
 	"	mov	"J", "N"				\n"
 	"	cmp	"J", xzr				\n"
@@ -104,6 +107,8 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"	ldr	d4, ["X"]				\n"
 	"	fabs	d4, d4					\n"
 	"	fmax	"CUR_MAX", "SCALE", d4			\n"
+	"	fcmp	"CUR_MAX", "REGINF"			\n"
+	"	beq	10f					\n"
 	"	fdiv	"SCALE", "SCALE", "CUR_MAX"		\n"
 	"	fmul	"SCALE", "SCALE", "SCALE"		\n"
 	"	fmul	"SSQ", "SSQ", "SCALE"			\n"
@@ -116,6 +121,8 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"	ldr	d3, ["X", #8]				\n"
 	"	fabs	d3, d3					\n"
 	"	fmax	"CUR_MAX", "SCALE", d3			\n"
+	"	fcmp	"CUR_MAX", "REGINF"			\n"
+	"	beq	10f					\n"
 	"	fdiv	"SCALE", "SCALE", "CUR_MAX"		\n"
 	"	fmul	"SCALE", "SCALE", "SCALE"		\n"
 	"	fmul	"SSQ", "SSQ", "SCALE"			\n"
@@ -158,6 +165,8 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"	fmaxp	v24.2d, v24.2d, v26.2d			\n"
 	"	fmaxp	v24.2d, v24.2d, v24.2d			\n"
 	"	fmax	"CUR_MAX", "SCALE", d24			\n"
+	"	fcmp	"CUR_MAX", "REGINF"			\n"
+	"	beq	10f					\n"
 	"	fdiv	"CUR_MAXINV", "REGONE", "CUR_MAX"	\n"
 	"	//dup	"CUR_MAX_V", v7.d[0]			\n"
 	"	fdiv	"SCALE", "SCALE", "CUR_MAX"		\n"
@@ -217,6 +226,8 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"	fmaxp	v24.2d, v24.2d, v26.2d			\n"
 	"	fmaxp	v24.2d, v24.2d, v24.2d			\n"
 	"	fmax	"CUR_MAX", "SCALE", d24			\n"
+	"	fcmp	"CUR_MAX", "REGINF"			\n"
+	"	beq	10f					\n"
 	"	fdiv	"CUR_MAXINV", "REGONE", "CUR_MAX"	\n"
 	"	//dup	"CUR_MAX_V", v7.d[0]			\n"
 	"	fdiv	"SCALE", "SCALE", "CUR_MAX"		\n"
@@ -265,6 +276,8 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"	ldr	d4, ["X"]				\n"
 	"	fabs	d4, d4					\n"
 	"	fmax	"CUR_MAX", "SCALE", d4			\n"
+	"	fcmp	"CUR_MAX", "REGINF"			\n"
+	"	beq	10f					\n"
 	"	fdiv	"SCALE", "SCALE", "CUR_MAX"		\n"
 	"	fmul	"SCALE", "SCALE", "SCALE"		\n"
 	"	fmul	"SSQ", "SSQ", "SCALE"			\n"
@@ -276,6 +289,8 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"	ldr	d3, ["X", #8]				\n"
 	"	fabs	d3, d3					\n"
 	"	fmax	"CUR_MAX", "SCALE", d3			\n"
+	"	fcmp	"CUR_MAX", "REGINF"			\n"
+	"	beq	10f					\n"
 	"	fdiv	"SCALE", "SCALE", "CUR_MAX"		\n"
 	"	fmul	"SCALE", "SCALE", "SCALE"		\n"
 	"	fmul	"SSQ", "SSQ", "SCALE"			\n"
@@ -291,6 +306,11 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	"9: //nrm2_kernel_L999:					\n"
 	"	str	"SSQ", [%[SSQ_]]			\n"
 	"	str	"SCALE", [%[SCALE_]]			\n"
+	"	b	11f					\n"
+	"10:							\n"
+	"	str	"REGINF", [%[SSQ_]]			\n"
+	"	str	"REGINF", [%[SCALE_]]			\n"
+	"11:							\n"
 
 	:
 	: [SSQ_]    "r"  (ssq),			//%0
@@ -300,7 +320,7 @@ static void nrm2_compute(BLASLONG n, FLOAT *x, BLASLONG inc_x,
 	  [INCX_]   "r"  (inc_x)		//%4
 	: "cc",
 	  "memory",
-	  "x0", "x1", "x2", "x3", "x4", "x5",
+	  "x0", "x1", "x2", "x3", "x4", "x5", "x6",
 	  "d0", "d1", "d2", "d3", "d4", "d5", "d6", "d7", "d8"
 	);
 
@@ -359,6 +379,12 @@ FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x)
 			cur_ssq = *ptr;
 			cur_scale = *(ptr + 1);
 
+			if (cur_ssq == INFINITY) {
+				ssq = INFINITY;
+				scale = INFINITY;
+				break;
+			}
+
 			if (cur_scale != 0) {
 				if (cur_scale > scale) {
 					scale = (scale / cur_scale);

kernel/power/KERNEL.POWER10

@@ -154,11 +154,7 @@ ZCOPYKERNEL  = zcopy_power10.c
 SDOTKERNEL   =  sdot_power10.c
 DDOTKERNEL   =  ddot_power10.c
 DSDOTKERNEL  =  sdot_power10.c
-ifneq ($(GCCVERSIONGTEQ9),1)
-CDOTKERNEL   =  cdot_power9.S
-else
 CDOTKERNEL   =  cdot.c
-endif
 ZDOTKERNEL   =  zdot.c
 #
 SNRM2KERNEL  = ../arm/nrm2.c

kernel/power/cdot.c

@@ -28,6 +28,9 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #else
 
 #include "common.h"
+#if defined(POWER10)
+#include "cdot_microk_power10.c"
+#else
 #ifndef HAVE_KERNEL_8
 #include <altivec.h> 
 
@@ -99,6 +102,7 @@ static void cdot_kernel_8(BLASLONG n, FLOAT *x, FLOAT *y, float *dot)
  
 }
 #endif
+#endif
  
 
 OPENBLAS_COMPLEX_FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y) {
@@ -116,7 +120,11 @@ OPENBLAS_COMPLEX_FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLA
 
     if ((inc_x == 1) && (inc_y == 1)) {
 
+#if defined(POWER10)
+        BLASLONG n1 = n & -16;
+#else
         BLASLONG n1 = n & -8;
+#endif
         BLASLONG j=0; 
 
         if (n1){

kernel/power/cdot_microk_power10.c

@@ -0,0 +1,177 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+
+#define HAVE_KERNEL_8 1
+
+static void cdot_kernel_8 (long n, float *x, float *y, float *dot)
+{
+  __vector unsigned char mask = { 11,10,9,8,15,14,13,12,3,2,1,0,7,6,5,4};
+  __asm__
+    (
+       "dcbt		0, %2		\n\t"
+       "dcbt		0, %3		\n\t"
+
+       "xxlxor		32, 32,	32	\n\t"
+       "xxlxor		33, 33,	33	\n\t"
+       "xxlxor		34, 34,	34	\n\t"
+       "xxlxor		35, 35,	35	\n\t"
+       "xxlxor		36, 36,	36	\n\t"
+       "xxlxor		37, 37,	37	\n\t"
+       "xxlxor		38, 38,	38	\n\t"
+       "xxlxor		39, 39,	39	\n\t"
+
+       "lxvp            40, 0(%2)       \n\t"
+       "lxvp            42, 32(%2)      \n\t"
+       "lxvp            44, 64(%2)      \n\t"
+       "lxvp            46, 96(%2)      \n\t"
+       "lxvp            48, 0(%3)       \n\t"
+       "lxvp            50, 32(%3)      \n\t"
+       "lxvp            52, 64(%3)      \n\t"
+       "lxvp            54, 96(%3)      \n\t"
+
+       "xxperm          56, 48, %x7     \n\t"
+       "xxperm          57, 49, %x7     \n\t"
+       "xxperm          58, 50, %x7     \n\t"
+       "xxperm          59, 51, %x7     \n\t"
+
+       "xxperm          60, 52, %x7     \n\t"
+       "xxperm          61, 53, %x7     \n\t"
+       "xxperm          62, 54, %x7     \n\t"
+       "xxperm          63, 55, %x7     \n\t"
+
+       "addi		%2, %2, 128	\n\t"
+       "addi		%3, %3, 128	\n\t"
+
+       "addic.		%1, %1, -16	\n\t"
+       "ble		two%=		\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "xvmaddasp	32, 40, 48	\n\t"	// x0_r * y0_r , x0_i * y0_i
+       "xvmaddasp	34, 41, 49	\n\t"	// x1_r * y1_r , x1_i * y1_i
+       "lxvp            48, 0(%3)       \n\t"
+
+       "xvmaddasp	36, 42, 50	\n\t"	// x2_r * y2_r , x2_i * y2_i
+       "xvmaddasp	38, 43, 51	\n\t"	// x3_r * y3_r , x3_i * y3_i
+       "lxvp            50, 32(%3)      \n\t"
+
+       "xvmaddasp	33, 40, 56 	\n\t"	// x0_r * y0_i , x0_i * y0_r
+       "xvmaddasp	35, 41, 57 	\n\t"	// x1_r * y1_i , x1_i * y1_r
+       "lxvp            40, 0(%2)       \n\t"
+
+       "xvmaddasp	37, 42, 58	\n\t"	// x2_r * y2_i , x2_i * y2_r
+       "xvmaddasp	39, 43, 59	\n\t"	// x3_r * y3_i , x3_i * y3_r
+       "lxvp            42, 32(%2)      \n\t"
+
+       "xxperm          56, 48, %x7     \n\t"
+       "xxperm          57, 49, %x7     \n\t"
+       "xxperm          58, 50, %x7     \n\t"
+       "xxperm          59, 51, %x7     \n\t"
+
+       "xvmaddasp	32, 44, 52	\n\t"	// x0_r * y0_r , x0_i * y0_i
+       "xvmaddasp	34, 45, 53	\n\t"	// x1_r * y1_r , x1_i * y1_i
+       "lxvp            52, 64(%3)      \n\t"
+
+       "xvmaddasp	36, 46, 54 	\n\t"	// x2_r * y2_r , x2_i * y2_i
+       "xvmaddasp	38, 47, 55	\n\t"	// x3_r * y3_r , x3_i * y3_i
+       "lxvp            54, 96(%3)      \n\t"
+
+       "xvmaddasp	33, 44, 60	\n\t"	// x0_r * y0_i , x0_i * y0_r
+       "xvmaddasp	35, 45, 61	\n\t"	// x1_r * y1_i , x1_i * y1_r
+       "lxvp            44, 64(%2)      \n\t"
+       "xvmaddasp	37, 46, 62	\n\t"	// x2_r * y2_i , x2_i * y2_r
+       "xvmaddasp	39, 47, 63	\n\t"	// x3_r * y3_i , x3_i * y3_r
+       "lxvp            46, 96(%2)      \n\t"
+
+       "xxperm          60, 52, %x7     \n\t"
+       "xxperm          61, 53, %x7     \n\t"
+       "xxperm          62, 54, %x7     \n\t"
+       "xxperm          63, 55, %x7     \n\t"
+
+       "addi		%2, %2, 128 	\n\t"
+       "addi		%3, %3, 128	\n\t"
+
+       "addic.		%1, %1, -16	\n\t"
+       "bgt		one%=		\n"
+
+     "two%=:				\n\t"
+
+       "xvmaddasp	32, 40, 48	\n\t"	// x0_r * y0_r , x0_i * y0_i
+       "xvmaddasp	34, 41, 49	\n\t"	// x1_r * y1_r , x1_i * y1_i
+       "xvmaddasp	36, 42, 50	\n\t"	// x2_r * y2_r , x2_i * y2_i
+       "xvmaddasp	38, 43, 51	\n\t"	// x3_r * y3_r , x3_i * y3_i
+
+       "xvmaddasp	33, 40, 56	\n\t"	// x0_r * y0_i , x0_i * y0_r
+       "xvmaddasp	35, 41, 57	\n\t"	// x1_r * y1_i , x1_i * y1_r
+       "xvmaddasp	37, 42, 58	\n\t"	// x2_r * y2_i , x2_i * y2_r
+       "xvmaddasp	39, 43, 59	\n\t"	// x3_r * y3_i , x3_i * y3_r
+
+       "xvmaddasp	32, 44, 52	\n\t"	// x0_r * y0_r , x0_i * y0_i
+       "xvmaddasp	34, 45, 53	\n\t"	// x1_r * y1_r , x1_i * y1_i
+       "xvmaddasp	36, 46, 54	\n\t"	// x2_r * y2_r , x2_i * y2_i
+       "xvmaddasp	38, 47, 55	\n\t"	// x3_r * y3_r , x3_i * y3_i
+
+       "xvmaddasp	33, 44, 60	\n\t"	// x0_r * y0_i , x0_i * y0_r
+       "xvmaddasp	35, 45, 61	\n\t"	// x1_r * y1_i , x1_i * y1_r
+       "xvmaddasp	37, 46, 62	\n\t"	// x2_r * y2_i , x2_i * y2_r
+       "xvmaddasp	39, 47, 63	\n\t"	// x3_r * y3_i , x3_i * y3_r
+
+       "xvaddsp		32, 32, 34	\n\t"
+       "xvaddsp		36, 36, 38	\n\t"
+
+       "xvaddsp		33, 33, 35	\n\t"
+       "xvaddsp		37, 37, 39	\n\t"
+
+       "xvaddsp		35, 32, 36	\n\t"
+       "xvaddsp		34, 33, 37	\n\t"
+       "xxswapd		32, 35		\n\t"
+       "xxswapd		33, 34		\n\t"
+       "xvaddsp		35, 35, 32	\n\t"
+       "xvaddsp		34, 34, 33	\n\t"
+       "xxpermdi	34, 34, 35, 2	\n\t"
+       "stxv		34, 0(%6)       \n\t"
+
+     "#n=%1 x=%4=%2 y=%5=%3 dot=%0=%6"
+     :
+       "=m" (*dot),
+       "+r" (n),	// 1
+       "+b" (x),	// 2
+       "+b" (y)		// 3
+     :
+       "m" (*x),
+       "m" (*y),
+       "b" (dot),	// 6
+       "wa" (mask)
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47",
+       "vs48","vs49","vs50","vs51","vs52","vs53","vs54","vs55",
+       "vs56","vs57","vs58","vs59","vs60","vs61","vs62","vs63"
+     );
+}

kernel/power/cscal_microk_power10.c

@@ -0,0 +1,176 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+
+#define HAVE_KERNEL_8 1
+
+static void zscal_kernel_8 (long n, float *x, float alpha_r, float alpha_i)
+{
+  __vector float t0 = {-alpha_i, alpha_i, -alpha_i, alpha_i};
+  __vector unsigned char mask = { 11,10,9,8,15,14,13,12,3,2,1,0,7,6,5,4};
+  __asm__
+    (
+       "dcbt		0, %2		\n\t"
+       "xscvdpspn	32, %x3    \n\t"
+       "xxspltw		32, 32, 0	\n\t"
+
+       "lxvp		40, 0(%2)	\n\t"
+       "lxvp		42, 32(%2)	\n\t"
+       "lxvp		44, 64(%2)	\n\t"
+       "lxvp		46, 96(%2)	\n\t"
+
+       "addic.		%1, %1, -16	\n\t"
+       "ble		two%=		\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "xvmulsp		48, 40, 32	\n\t"	// x0_r * alpha_r, x0_i * alpha_r
+       "xvmulsp		49, 41, 32	\n\t"
+       "xvmulsp		50, 42, 32	\n\t"
+       "xvmulsp		51, 43, 32	\n\t"
+       "xvmulsp		52, 44, 32	\n\t"
+       "xvmulsp		53, 45, 32	\n\t"
+       "xvmulsp		54, 46, 32	\n\t"
+       "xvmulsp		55, 47, 32	\n\t"
+
+       "xxperm 34, 40, %x5 \n\t"
+       "xxperm 35, 41, %x5 \n\t"
+       "xxperm 36, 42, %x5 \n\t"
+       "xxperm 37, 43, %x5 \n\t"
+       "xxperm 38, 44, %x5 \n\t"
+       "xxperm 39, 45, %x5 \n\t"
+       "xxperm 56, 46, %x5 \n\t"
+       "xxperm 57, 47, %x5 \n\t"
+
+       "xvmulsp		34, 34, %x4	\n\t"	// x0_i * -alpha_i, x0_r * alpha_i
+       "xvmulsp		35, 35, %x4	\n\t"
+
+       "lxvp		40, 128(%2)	\n\t"
+
+       "xvmulsp		36, 36, %x4	\n\t"
+       "xvmulsp		37, 37, %x4	\n\t"
+
+       "lxvp		42, 160(%2)	\n\t"
+
+       "xvmulsp		38, 38, %x4	\n\t"
+       "xvmulsp		39, 39, %x4	\n\t"
+
+       "lxvp		44, 192(%2)	\n\t"
+
+       "xvmulsp		56, 56, %x4	\n\t"
+       "xvmulsp		57, 57, %x4	\n\t"
+
+       "lxvp		46, 224(%2)	\n\t"
+
+       "xvaddsp		48, 48, 34	\n\t"
+       "xvaddsp		49, 49, 35	\n\t"
+       "xvaddsp		50, 50, 36	\n\t"
+       "xvaddsp		51, 51, 37	\n\t"
+
+       "stxvp		48, 0(%2)	\n\t"
+
+       "xvaddsp		52, 52, 38	\n\t"
+       "xvaddsp		53, 53, 39	\n\t"
+
+       "stxvp		50, 32(%2)	\n\t"
+
+       "xvaddsp		54, 54, 56	\n\t"
+       "xvaddsp		55, 55, 57	\n\t"
+
+       "stxvp		52, 64(%2)	\n\t"
+       "stxvp		54, 96(%2)	\n\t"
+
+       "addi		%2, %2, 128	\n\t"
+
+       "addic.		%1, %1, -16	\n\t"
+       "bgt		one%=		\n"
+
+     "two%=:				\n\t"
+
+       "xvmulsp		48, 40, 32	\n\t"	// x0_r * alpha_r, x0_i * alpha_r
+       "xvmulsp		49, 41, 32	\n\t"
+       "xvmulsp		50, 42, 32	\n\t"
+       "xvmulsp		51, 43, 32	\n\t"
+       "xvmulsp		52, 44, 32	\n\t"
+       "xvmulsp		53, 45, 32	\n\t"
+       "xvmulsp		54, 46, 32	\n\t"
+       "xvmulsp		55, 47, 32	\n\t"
+
+       "xxperm 34, 40, %x5 \n\t"
+       "xxperm 35, 41, %x5 \n\t"
+       "xxperm 36, 42, %x5 \n\t"
+       "xxperm 37, 43, %x5 \n\t"
+       "xxperm 38, 44, %x5 \n\t"
+       "xxperm 39, 45, %x5 \n\t"
+       "xxperm 56, 46, %x5 \n\t"
+       "xxperm 57, 47, %x5 \n\t"
+
+
+       "xvmulsp		34, 34, %x4	\n\t"	// x0_i * -alpha_i, x0_r * alpha_i
+       "xvmulsp		35, 35, %x4	\n\t"
+       "xvmulsp		36, 36, %x4	\n\t"
+       "xvmulsp		37, 37, %x4	\n\t"
+       "xvmulsp		38, 38, %x4	\n\t"
+       "xvmulsp		39, 39, %x4	\n\t"
+       "xvmulsp		56, 56, %x4	\n\t"
+       "xvmulsp		57, 57, %x4	\n\t"
+
+       "xvaddsp		48, 48, 34	\n\t"
+       "xvaddsp		49, 49, 35	\n\t"
+       "xvaddsp		50, 50, 36	\n\t"
+       "xvaddsp		51, 51, 37	\n\t"
+
+       "stxvp		48, 0(%2)	\n\t"
+
+       "xvaddsp		52, 52, 38	\n\t"
+       "xvaddsp		53, 53, 39	\n\t"
+
+       "stxvp		50, 32(%2)	\n\t"
+
+       "xvaddsp		54, 54, 56	\n\t"
+       "xvaddsp		55, 55, 57	\n\t"
+
+       "stxvp		52, 64(%2)	\n\t"
+       "stxvp		54, 96(%2)	\n\t"
+
+     "#n=%1 x=%0=%2 alpha=(%3,%4)\n"
+     :
+       "+m" (*x),
+       "+r" (n),	// 1
+       "+b" (x)		// 2
+     :
+       "f" (alpha_r),	// 3 
+       "wa" (t0),	// 4 
+       "wa" (mask)	// 5
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47",
+       "vs48","vs49","vs50","vs51","vs52","vs53","vs54","vs55",
+       "vs56","vs57"
+     );
+}

kernel/power/cswap.c

@@ -36,9 +36,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "common.h"
 
 
-#if defined(POWER8)  || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8)  || defined(POWER9)
 #include "cswap_microk_power8.c"
+#elif defined(POWER10)
+#include "cswap_microk_power10.c"
 #endif
 #endif
 

kernel/power/cswap_microk_power10.c

@@ -0,0 +1,127 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+
+#if defined(DOUBLE)
+#define HAVE_KERNEL_16 1
+static void zswap_kernel_16 (long n, double *x, double *y)
+#else
+#define HAVE_KERNEL_32 1
+static void cswap_kernel_32 (long n, float *x, float *y)
+#endif
+{
+  __asm__
+    (
+       ".align	5		\n"
+     "one%=:				\n\t"
+       "lxvp            32, 0(%4)       \n\t"
+       "lxvp            34, 32(%4)      \n\t"
+       "lxvp            36, 64(%4)      \n\t"
+       "lxvp            38, 96(%4)      \n\t"
+
+       "lxvp            40, 128(%4)     \n\t"
+       "lxvp            42, 160(%4)     \n\t"
+       "lxvp            44, 192(%4)     \n\t"
+       "lxvp            46, 224(%4)     \n\t"
+
+       "lxvp            48, 0(%3)       \n\t"
+       "lxvp            50, 32(%3)      \n\t"
+       "lxvp            52, 64(%3)      \n\t"
+       "lxvp            54, 96(%3)      \n\t"
+
+       "lxvp            56, 128(%3)     \n\t"
+       "lxvp            58, 160(%3)     \n\t"
+       "lxvp            60, 192(%3)     \n\t"
+       "lxvp            62, 224(%3)     \n\t"
+
+
+       "stxv		33, 0(%3)	\n\t"
+       "stxv		32, 16(%3)	\n\t"
+       "stxv		35, 32(%3)	\n\t"
+       "stxv		34, 48(%3)	\n\t"
+       "stxv		37, 64(%3)	\n\t"
+       "stxv		36, 80(%3)	\n\t"
+       "stxv		39, 96(%3)	\n\t"
+       "stxv		38, 112(%3)	\n\t"
+
+       "addi		%3, %3, 128	\n\t"
+
+       "stxv		41, 0(%3)	\n\t"
+       "stxv		40, 16(%3)	\n\t"
+       "stxv		43, 32(%3)	\n\t"
+       "stxv		42, 48(%3)	\n\t"
+       "stxv		45, 64(%3)	\n\t"
+       "stxv		44, 80(%3)	\n\t"
+       "stxv		47, 96(%3)	\n\t"
+       "stxv		46, 112(%3)	\n\t"
+
+       "addi		%3, %3, 128	\n\t"
+
+       "stxv		49, 0(%4)	\n\t"
+       "stxv		48, 16(%4)	\n\t"
+       "stxv		51, 32(%4)	\n\t"
+       "stxv		50, 48(%4)	\n\t"
+       "stxv		53, 64(%4)	\n\t"
+       "stxv		52, 80(%4)	\n\t"
+       "stxv		55, 96(%4)	\n\t"
+       "stxv		54, 112(%4)	\n\t"
+
+       "addi		%4, %4, 128	\n\t"
+
+       "stxv		57, 0(%4)	\n\t"
+       "stxv		56, 16(%4)	\n\t"
+       "stxv		59, 32(%4)	\n\t"
+       "stxv		58, 48(%4)	\n\t"
+       "stxv		61, 64(%4)	\n\t"
+       "stxv		60, 80(%4)	\n\t"
+       "stxv		63, 96(%4)	\n\t"
+       "stxv		62, 112(%4)	\n\t"
+
+       "addi		%4, %4, 128	\n\t"
+
+#if defined(DOUBLE)
+       "addic.		%2, %2, -16	\n\t"
+#else
+       "addic.		%2, %2, -32	\n\t"
+#endif
+       "bgt		one%=		\n"
+
+     "#n=%2 x=%0=%3 y=%1=%4"
+     :
+       "+m" (*x),
+       "+m" (*y),
+       "+r" (n),	// 2
+       "+b" (x),	// 3
+       "+b" (y)		// 4
+     :
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47",
+       "vs48","vs49","vs50","vs51","vs52","vs53","vs54","vs55",
+       "vs56","vs57","vs58","vs59","vs60","vs61","vs62","vs63"
+     );
+}

kernel/power/drot.c

@@ -39,9 +39,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #pragma GCC optimize "O1"
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #include "drot_microk_power8.c"
+#elif defined(POWER10)
+#include "drot_microk_power10.c"
 #endif
 #endif
 
@@ -115,12 +117,30 @@ int CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y, FLOAT
 	if ( (inc_x == 1) && (inc_y == 1) )
 	{
 
+#if defined(POWER10)
+		if ( n >= 16 )
+		{
+			BLASLONG align = ((32 - ((uintptr_t)y & (uintptr_t)0x1F)) >> 3) & 0x3;
+			for (i = 0; i < align; i++) {
+				temp  = c*x[i] + s*y[i] ;
+				y[i]  = c*y[i] - s*x[i] ;
+				x[i]  = temp ;
+			}
+		}
+		BLASLONG n1 = (n-i) & -16;
+		if ( n1 > 0 )
+		{
+			drot_kernel_16(n1,&x[i], &y[i], c, s);
+			i+=n1;
+		}
+#else
 		BLASLONG n1 = n & -16;
 		if ( n1 > 0 )
 		{
 			drot_kernel_16(n1, x1, y1, c, s);
 			i=n1;
 		}
+#endif
 
 		while(i < n)
 		{

kernel/power/drot_microk_power10.c

@@ -0,0 +1,148 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+
+#define HAVE_KERNEL_16 1
+
+static void drot_kernel_16 (long n, double *x, double *y, double c, double s)
+{
+  __asm__
+    (
+       XXSPLTD_S(36,%x5,0)	// load c to both dwords
+       XXSPLTD_S(37,%x6,0)	// load s to both dwords
+       "lxvp            32, 0(%3)       \n\t"   // load x
+       "lxvp            34, 32(%3)      \n\t"
+       "lxvp            48, 0(%4)       \n\t"   // load y
+       "lxvp            50, 32(%4)      \n\t"
+
+       "addic.		%2, %2, -8	\n\t"
+       "ble		two%=		\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "xvmuldp		40, 32, 36	\n\t"	// c * x
+       "xvmuldp		41, 33, 36	\n\t"
+       "xvmuldp		42, 34, 36	\n\t"
+       "xvmuldp		43, 35, 36	\n\t"
+
+       "xvmuldp		44, 32, 37	\n\t"	// s * x
+       "xvmuldp		45, 33, 37	\n\t"
+       "xvmuldp		46, 34, 37	\n\t"
+       "xvmuldp		47, 35, 37	\n\t"
+
+       "lxvp            32, 64(%3)       \n\t"   // load x
+       "lxvp            34, 96(%3)      \n\t"
+       "xvmuldp		52, 48, 36	\n\t"	// c * y
+       "xvmuldp		53, 49, 36	\n\t"
+       "xvmuldp		54, 50, 36	\n\t"
+       "xvmuldp		55, 51, 36	\n\t"
+
+       "xvmuldp		38, 48, 37	\n\t"	// s * y
+       "xvmuldp		39, 49, 37	\n\t"
+       "xvmuldp		56, 50, 37	\n\t"
+       "xvmuldp		57, 51, 37	\n\t"
+
+       "lxvp            48, 64(%4)       \n\t"   // load y
+       "lxvp            50, 96(%4)      \n\t"
+
+       "xvadddp		40, 40, 38	\n\t"	// c * x + s * y
+       "xvadddp		41, 41, 39	\n\t"	// c * x + s * y
+       "xvadddp		42, 42, 56	\n\t"	// c * x + s * y
+       "xvadddp		43, 43, 57	\n\t"	// c * x + s * y
+
+       "stxvp           40, 0(%3)       \n\t"   // store x
+       "stxvp           42, 32(%3)      \n\t"
+
+       "xvsubdp         52, 52, 44      \n\t"   // c * y - s * x
+       "xvsubdp         53, 53, 45      \n\t"   // c * y - s * x
+       "xvsubdp         54, 54, 46      \n\t"   // c * y - s * x
+       "xvsubdp         55, 55, 47      \n\t"   // c * y - s * x
+
+       "stxvp           52, 0(%4)       \n\t"   // store y
+       "stxvp           54, 32(%4)      \n\t"
+
+       "addi		%3, %3, 64	\n\t"
+       "addi		%4, %4, 64	\n\t"
+
+       "addic.		%2, %2, -8	\n\t"
+       "bgt		one%=		\n"
+
+     "two%=:				\n\t"
+
+       "xvmuldp		40, 32, 36	\n\t"	// c * x
+       "xvmuldp		41, 33, 36	\n\t"
+       "xvmuldp		42, 34, 36	\n\t"
+       "xvmuldp		43, 35, 36	\n\t"
+
+       "xvmuldp         52, 48, 36      \n\t"   // c * y
+       "xvmuldp         53, 49, 36      \n\t"
+       "xvmuldp         54, 50, 36      \n\t"
+       "xvmuldp         55, 51, 36      \n\t"
+
+       "xvmuldp         44, 32, 37      \n\t"   // s * x
+       "xvmuldp         45, 33, 37      \n\t"
+       "xvmuldp         46, 34, 37      \n\t"
+       "xvmuldp         47, 35, 37      \n\t"
+
+       "xvmuldp         38, 48, 37     \n\t"   // s * y
+       "xvmuldp         39, 49, 37     \n\t"
+       "xvmuldp         56, 50, 37     \n\t"
+       "xvmuldp         57, 51, 37     \n\t"
+
+       "xvadddp         40, 40, 38     \n\t"   // c * x + s * y
+       "xvadddp         41, 41, 39     \n\t"   // c * x + s * y
+       "xvadddp         42, 42, 56     \n\t"   // c * x + s * y
+       "xvadddp         43, 43, 57     \n\t"   // c * x + s * y
+
+       "stxvp           40, 0(%3)       \n\t"   // store x
+       "stxvp           42, 32(%3)      \n\t"
+       "xvsubdp         52, 52, 44      \n\t"   // c * y - s * x
+       "xvsubdp         53, 53, 45      \n\t"   // c * y - s * x
+       "xvsubdp         54, 54, 46      \n\t"   // c * y - s * x
+       "xvsubdp         55, 55, 47      \n\t"   // c * y - s * x
+
+       "stxvp           52, 0(%4)       \n\t"   // store y
+       "stxvp           54, 32(%4)      \n\t"
+
+     "#n=%2 x=%0=%3 y=%1=%4 c=%5 s=%6\n"
+     :
+       "+m" (*x),
+       "+m" (*y),
+       "+r" (n),	// 2
+       "+b" (x),	// 3
+       "+b" (y) 	// 4
+     :
+       "d" (c),		// 5 
+       "d" (s)		// 6 
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47",
+       "vs48","vs49","vs50","vs51","vs52","vs53","vs54","vs55",
+       "vs56","vs57"
+     );
+}

kernel/power/dscal.c

@@ -35,9 +35,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "common.h"
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #include "dscal_microk_power8.c"
+#elif defined(POWER10)
+#include "dscal_microk_power10.c"
 #endif
 #endif
 
@@ -100,12 +102,28 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 		if ( da == 0.0 )
 		{		
 
+#if defined(POWER10)
+			if ( n >= 16 )
+			{
+				BLASLONG align = ((32 - ((uintptr_t)x & (uintptr_t)0x1F)) >> 3) & 0x3;
+				for (j = 0; j < align; j++) {
+					x[j] = 0.0;
+				}
+			}
+			BLASLONG n1 = (n-j) & -16;
+			if ( n1 > 0 )
+			{
+				dscal_kernel_8_zero(n1, &x[j]);
+				j+=n1;
+			}
+#else
 			BLASLONG n1 = n & -16;
 			if ( n1 > 0 )
 			{
 				dscal_kernel_8_zero(n1, x);
 				j=n1;
 			}
+#endif
 
 			while(j < n)
 			{
@@ -118,12 +136,28 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 		else
 		{
 
+#if defined(POWER10)
+			if ( n >= 16 )
+			{
+				BLASLONG align = ((32 - ((uintptr_t)x & (uintptr_t)0x1F)) >> 3) & 0x3;
+				for (j = 0; j < align; j++) {
+					x[j] = da * x[j];
+				}
+			}
+			BLASLONG n1 = (n-j) & -16;
+			if ( n1 > 0 )
+			{
+				dscal_kernel_8(n1, &x[j], da);
+				j+=n1;
+			}
+#else
 			BLASLONG n1 = n & -16;
 			if ( n1 > 0 )
 			{
 				dscal_kernel_8(n1, x, da);
 				j=n1;
 			}
+#endif
 			while(j < n)
 			{
 

kernel/power/dscal_microk_power10.c

@@ -0,0 +1,134 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+
+#define HAVE_KERNEL_8 1
+
+static void dscal_kernel_8 (long n, double *x, double alpha)
+{
+  __asm__
+    (
+       "dcbt		0, %2		\n\t"
+
+       XXSPLTD_S(48,%x3,0)
+
+       "lxvp		32, 0(%2)	\n\t"
+       "lxvp		34, 32(%2)	\n\t"
+       "lxvp		36, 64(%2)	\n\t"
+       "lxvp		38, 96(%2)	\n\t"
+
+       "addic.		%1, %1, -16	\n\t"
+       "ble		two%=		\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "xvmuldp		40, 32, 48	\n\t"
+       "xvmuldp		41, 33, 48	\n\t"
+       "xvmuldp		42, 34, 48	\n\t"
+       "xvmuldp		43, 35, 48	\n\t"
+       "lxvp		32, 128(%2)	\n\t"
+       "lxvp		34, 160(%2)	\n\t"
+       "xvmuldp		44, 36, 48	\n\t"
+       "xvmuldp		45, 37, 48	\n\t"
+       "xvmuldp		46, 38, 48	\n\t"
+       "xvmuldp		47, 39, 48	\n\t"
+       "lxvp		36, 192(%2)	\n\t"
+       "lxvp		38, 224(%2)	\n\t"
+
+       "stxvp		40, 0(%2)	\n\t"
+       "stxvp		42, 32(%2)	\n\t"
+       "stxvp		44, 64(%2)	\n\t"
+       "stxvp		46, 96(%2)	\n\t"
+
+       "addi		%2, %2, 128	\n\t"
+
+       "addic.		%1, %1, -16	\n\t"
+       "bgt		one%=		\n"
+
+     "two%=:				\n\t"
+
+       "xvmuldp		40, 32, 48	\n\t"
+       "xvmuldp		41, 33, 48	\n\t"
+       "xvmuldp		42, 34, 48	\n\t"
+       "xvmuldp		43, 35, 48	\n\t"
+
+       "xvmuldp		44, 36, 48	\n\t"
+       "xvmuldp		45, 37, 48	\n\t"
+       "xvmuldp		46, 38, 48	\n\t"
+       "xvmuldp		47, 39, 48	\n\t"
+
+       "stxvp		40, 0(%2)	\n\t"
+       "stxvp		42, 32(%2)	\n\t"
+       "stxvp		44, 64(%2)	\n\t"
+       "stxvp		46, 96(%2)	\n\t"
+
+     "#n=%1 alpha=%3 x=%0=%2"
+     :
+       "+m" (*x),
+       "+r" (n),	// 1
+       "+b" (x)		// 2
+     :
+       "d" (alpha)	// 3
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47","vs48"
+     );
+}
+
+
+static void dscal_kernel_8_zero (long n, double *x)
+{
+
+  __asm__
+    (
+       "xxlxor		32, 32, 32	\n\t"
+       "xxlxor		33, 33, 33	\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "stxvp		32, 0(%2)	\n\t"
+       "stxvp		32, 32(%2)	\n\t"
+       "stxvp		32, 64(%2)	\n\t"
+       "stxvp		32, 96(%2)	\n\t"
+
+       "addi		%2, %2, 128	\n\t"
+
+       "addic.		%1, %1, -16	\n\t"
+       "bgt		one%=		\n"
+
+     "#n=%1 x=%0=%2 "
+     :
+       "=m" (*x),
+       "+r" (n),	// 1
+       "+b" (x)	// 2
+     :
+     :
+       "cr0","vs32","vs33"
+     );
+}

kernel/power/dswap.c

@@ -35,9 +35,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "common.h"
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #include "dswap_microk_power8.c"
+#elif defined(POWER10)
+#include "swap_microk_power10.c"
 #endif
 #endif
 
@@ -115,12 +117,30 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT dummy3, FLOAT *x,
 	if ( (inc_x == 1) && (inc_y == 1 ))
 	{
 
+#if defined(POWER10)
+		if ( n >= 32 )
+		{
+			BLASLONG align = ((32 - ((uintptr_t)y & (uintptr_t)0x1F)) >> 3) & 0x3;
+			for (i = 0; i < align; i++) {
+				temp = y[i];
+				y[i] = x[i];
+				x[i] = temp;
+			}
+		}
+		BLASLONG n1 = (n-i) & -32;
+		if ( n1 > 0 )
+		{
+			dswap_kernel_32(n1,&x[i], &y[i]);
+			i+=n1;
+		}
+#else
 		BLASLONG n1 = n & -32;
 		if ( n1 > 0 )
 		{
 			dswap_kernel_32(n1, x, y);
 			i=n1;
 		}
+#endif
 
 		while(i < n)
 		{

kernel/power/srot.c

@@ -39,9 +39,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #pragma GCC optimize "O1"
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #include "srot_microk_power8.c"
+#elif defined(POWER10)
+#include "srot_microk_power10.c"
 #endif
 #endif
 
@@ -115,6 +117,23 @@ int CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y, FLOAT
 	if ( (inc_x == 1) && (inc_y == 1) )
 	{
 
+#if defined(POWER10)
+		if ( n >= 16 )
+		{
+			BLASLONG align = ((32 - ((uintptr_t)y & (uintptr_t)0x1F)) >> 2) & 0x7;
+			for (i = 0; i < align; i++) {
+				temp  = c*x[i] + s*y[i] ;
+				y[i]  = c*y[i] - s*x[i] ;
+				x[i]  = temp ;
+			}
+		}
+		BLASLONG n1 = (n-i) & -16;
+		if ( n1 > 0 )
+		{
+			srot_kernel_16(n1, &x1[i], &y1[i], c, s);
+			i+=n1;
+		}
+#else
 		BLASLONG n1 = n & -16;
 		if ( n1 > 0 )
 		{
@@ -122,6 +141,7 @@ int CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x, FLOAT *y, BLASLONG inc_y, FLOAT
 			i=n1;
 		}
 
+#endif
 		while(i < n)
 		{
 			temp  = c*x[i] + s*y[i] ;

kernel/power/srot_microk_power10.c

@@ -0,0 +1,151 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+
+#define HAVE_KERNEL_16 1
+
+static void srot_kernel_16 (long n, float *x, float *y, float c, float s)
+{
+  __asm__
+    (
+       "xscvdpspn       36, %x5        \n\t"   // load c to all words
+       "xxspltw         36, 36, 0       \n\t"
+
+       "xscvdpspn       37, %x6        \n\t"   // load s to all words
+       "xxspltw         37, 37, 0       \n\t"
+       "lxvp            32, 0(%3)       \n\t"   // load x
+       "lxvp            34, 32(%3)      \n\t"
+       "lxvp            48, 0(%4)       \n\t"   // load y
+       "lxvp            50, 32(%4)      \n\t"
+
+       "addic.		%2, %2, -16	\n\t"
+       "ble		two%=		\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "xvmulsp		40, 32, 36	\n\t"	// c * x
+       "xvmulsp		41, 33, 36	\n\t"
+       "xvmulsp		42, 34, 36	\n\t"
+       "xvmulsp		43, 35, 36	\n\t"
+
+       "xvmulsp		44, 32, 37	\n\t"	// s * x
+       "xvmulsp		45, 33, 37	\n\t"
+       "xvmulsp		46, 34, 37	\n\t"
+       "xvmulsp		47, 35, 37	\n\t"
+
+       "lxvp            32, 64(%3)       \n\t"   // load x
+       "lxvp            34, 96(%3)      \n\t"
+       "xvmulsp		52, 48, 36	\n\t"	// c * y
+       "xvmulsp		53, 49, 36	\n\t"
+       "xvmulsp		54, 50, 36	\n\t"
+       "xvmulsp		55, 51, 36	\n\t"
+
+       "xvmulsp		38, 48, 37	\n\t"	// s * y
+       "xvmulsp		39, 49, 37	\n\t"
+       "xvmulsp		56, 50, 37	\n\t"
+       "xvmulsp		57, 51, 37	\n\t"
+
+       "lxvp            48, 64(%4)       \n\t"   // load y
+       "lxvp            50, 96(%4)      \n\t"
+
+       "xvaddsp		40, 40, 38	\n\t"	// c * x + s * y
+       "xvaddsp		41, 41, 39	\n\t"	// c * x + s * y
+       "xvaddsp		42, 42, 56	\n\t"	// c * x + s * y
+       "xvaddsp		43, 43, 57	\n\t"	// c * x + s * y
+
+       "stxvp           40, 0(%3)       \n\t"   // store x
+       "stxvp           42, 32(%3)      \n\t"
+
+       "xvsubsp         52, 52, 44      \n\t"   // c * y - s * x
+       "xvsubsp         53, 53, 45      \n\t"   // c * y - s * x
+       "xvsubsp         54, 54, 46      \n\t"   // c * y - s * x
+       "xvsubsp         55, 55, 47      \n\t"   // c * y - s * x
+
+       "stxvp           52, 0(%4)       \n\t"   // store y
+       "stxvp           54, 32(%4)      \n\t"
+
+       "addi		%3, %3, 64	\n\t"
+       "addi		%4, %4, 64	\n\t"
+
+       "addic.		%2, %2, -16	\n\t"
+       "bgt		one%=		\n"
+
+     "two%=:				\n\t"
+
+       "xvmulsp		40, 32, 36	\n\t"	// c * x
+       "xvmulsp		41, 33, 36	\n\t"
+       "xvmulsp		42, 34, 36	\n\t"
+       "xvmulsp		43, 35, 36	\n\t"
+
+       "xvmulsp         52, 48, 36      \n\t"   // c * y
+       "xvmulsp         53, 49, 36      \n\t"
+       "xvmulsp         54, 50, 36      \n\t"
+       "xvmulsp         55, 51, 36      \n\t"
+
+       "xvmulsp         44, 32, 37      \n\t"   // s * x
+       "xvmulsp         45, 33, 37      \n\t"
+       "xvmulsp         46, 34, 37      \n\t"
+       "xvmulsp         47, 35, 37      \n\t"
+
+       "xvmulsp         38, 48, 37     \n\t"   // s * y
+       "xvmulsp         39, 49, 37     \n\t"
+       "xvmulsp         56, 50, 37     \n\t"
+       "xvmulsp         57, 51, 37     \n\t"
+
+       "xvaddsp         40, 40, 38     \n\t"   // c * x + s * y
+       "xvaddsp         41, 41, 39     \n\t"   // c * x + s * y
+       "xvaddsp         42, 42, 56     \n\t"   // c * x + s * y
+       "xvaddsp         43, 43, 57     \n\t"   // c * x + s * y
+
+       "stxvp           40, 0(%3)       \n\t"   // store x
+       "stxvp           42, 32(%3)      \n\t"
+       "xvsubsp         52, 52, 44      \n\t"   // c * y - s * x
+       "xvsubsp         53, 53, 45      \n\t"   // c * y - s * x
+       "xvsubsp         54, 54, 46      \n\t"   // c * y - s * x
+       "xvsubsp         55, 55, 47      \n\t"   // c * y - s * x
+
+       "stxvp           52, 0(%4)       \n\t"   // store y
+       "stxvp           54, 32(%4)      \n\t"
+
+     "#n=%2 x=%0=%3 y=%1=%4 c=%5 s=%6\n"
+     :
+       "+m" (*x),
+       "+m" (*y),
+       "+r" (n),	// 2
+       "+b" (x),	// 3
+       "+b" (y) 	// 4
+     :
+       "f" (c),		// 5 
+       "f" (s)		// 6 
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47",
+       "vs48","vs49","vs50","vs51","vs52","vs53","vs54","vs55",
+       "vs56","vs57"
+     );
+}

kernel/power/sscal.c

@@ -35,9 +35,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "common.h"
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #include "sscal_microk_power8.c"
+#elif defined(POWER10)
+#include "sscal_microk_power10.c"
 #endif
 #endif
 
@@ -102,12 +104,28 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 		if ( da == 0.0 )
 		{		
 
+#if defined(POWER10)
+			if ( n >= 32 )
+			{
+				BLASLONG align = ((32 - ((uintptr_t)x & (uintptr_t)0x1F)) >> 2) & 0x7;
+				for (j = 0; j < align; j++) {
+					x[j] = 0.0;
+				}
+			}
+			BLASLONG n1 = (n-j) & -32;
+			if ( n1 > 0 )
+			{
+				sscal_kernel_16_zero(n1, &x[j]);
+				j+=n1;
+			}
+#else
 			BLASLONG n1 = n & -32;
 			if ( n1 > 0 )
 			{
 				sscal_kernel_16_zero(n1, x);
 				j=n1;
 			}
+#endif
 
 			while(j < n)
 			{
@@ -120,12 +138,28 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da, FLOAT *x, BLAS
 		else
 		{
 
+#if defined(POWER10)
+			if ( n >= 32 )
+			{
+				BLASLONG align = ((32 - ((uintptr_t)x & (uintptr_t)0x1F)) >> 2) & 0x7;
+				for (j = 0; j < align; j++) {
+					x[j] = da * x[j];
+				}
+			}
+			BLASLONG n1 = (n-j) & -32;
+			if ( n1 > 0 )
+			{
+				sscal_kernel_16(n1, &x[j], da);
+				j+=n1;
+			}
+#else
 			BLASLONG n1 = n & -32;
 			if ( n1 > 0 )
 			{
 				sscal_kernel_16(n1, x, da);
 				j=n1;
 			}
+#endif
 			while(j < n)
 			{
 

kernel/power/sscal_microk_power10.c

@@ -0,0 +1,135 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+
+#define HAVE_KERNEL_16 1
+
+static void sscal_kernel_16 (long n, float *x, float alpha)
+{
+  __asm__
+    (
+       "dcbt		0, %2		\n\t"
+
+       "xscvdpspn	48, %x3	\n\t"
+       "xxspltw		48, 48, 0	\n\t"
+
+       "lxvp		32, 0(%2)	\n\t"
+       "lxvp		34, 32(%2)	\n\t"
+       "lxvp		36, 64(%2)	\n\t"
+       "lxvp		38, 96(%2)	\n\t"
+
+       "addic.		%1, %1, -32	\n\t"
+       "ble		two%=		\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "xvmulsp		40, 32, 48	\n\t"
+       "xvmulsp		41, 33, 48	\n\t"
+       "xvmulsp		42, 34, 48	\n\t"
+       "xvmulsp		43, 35, 48	\n\t"
+       "lxvp		32, 128(%2)	\n\t"
+       "lxvp		34, 160(%2)	\n\t"
+       "xvmulsp		44, 36, 48	\n\t"
+       "xvmulsp		45, 37, 48	\n\t"
+       "xvmulsp		46, 38, 48	\n\t"
+       "xvmulsp		47, 39, 48	\n\t"
+       "lxvp		36, 192(%2)	\n\t"
+       "lxvp		38, 224(%2)	\n\t"
+
+       "stxvp		40, 0(%2)	\n\t"
+       "stxvp		42, 32(%2)	\n\t"
+       "stxvp		44, 64(%2)	\n\t"
+       "stxvp		46, 96(%2)	\n\t"
+
+       "addi		%2, %2, 128	\n\t"
+
+       "addic.		%1, %1, -32	\n\t"
+       "bgt		one%=		\n"
+
+     "two%=:				\n\t"
+
+       "xvmulsp		40, 32, 48	\n\t"
+       "xvmulsp		41, 33, 48	\n\t"
+       "xvmulsp		42, 34, 48	\n\t"
+       "xvmulsp		43, 35, 48	\n\t"
+
+       "xvmulsp		44, 36, 48	\n\t"
+       "xvmulsp		45, 37, 48	\n\t"
+       "xvmulsp		46, 38, 48	\n\t"
+       "xvmulsp		47, 39, 48	\n\t"
+
+       "stxvp		40, 0(%2)	\n\t"
+       "stxvp		42, 32(%2)	\n\t"
+       "stxvp		44, 64(%2)	\n\t"
+       "stxvp		46, 96(%2)	\n\t"
+
+     "#n=%1 alpha=%3 x=%0=%2"
+     :
+       "+m" (*x),
+       "+r" (n),	// 1
+       "+b" (x)		// 2
+     :
+       "f" (alpha)	// 3
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47","vs48"
+     );
+}
+
+
+static void sscal_kernel_16_zero (long n, float *x)
+{
+
+  __asm__
+    (
+       "xxlxor		32, 32, 32	\n\t"
+       "xxlxor		33, 33, 33	\n\t"
+
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "stxvp		32, 0(%2)	\n\t"
+       "stxvp		32, 32(%2)	\n\t"
+       "stxvp		32, 64(%2)	\n\t"
+       "stxvp		32, 96(%2)	\n\t"
+
+       "addi		%2, %2, 128	\n\t"
+
+       "addic.		%1, %1, -32	\n\t"
+       "bgt		one%=		\n"
+
+     "#n=%1 x=%0=%2 "
+     :
+       "=m" (*x),
+       "+r" (n),	// 1
+       "+b" (x)	// 2
+     :
+     :
+       "cr0","vs32","vs33"
+     );
+}

kernel/power/sswap.c

@@ -35,9 +35,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "common.h"
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #include "sswap_microk_power8.c"
+#elif defined(POWER10)
+#include "swap_microk_power10.c"
 #endif
 #endif
 
@@ -115,12 +117,30 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT dummy3, FLOAT *x,
 	if ( (inc_x == 1) && (inc_y == 1 ))
 	{
 
+#if defined(POWER10)
+		if ( n >= 64 )
+		{
+			BLASLONG align = ((32 - ((uintptr_t)y & (uintptr_t)0x1F)) >> 2) & 0x7;
+			for (i = 0; i < align; i++) {
+				temp = y[i];
+				y[i] = x[i];
+				x[i] = temp;
+			}
+		}
+		BLASLONG n1 = (n-i) & -64;
+		if ( n1 > 0 )
+		{
+			sswap_kernel_32(n1,&x[i], &y[i]);
+			i+=n1;
+		}
+#else
 		BLASLONG n1 = n & -32;
 		if ( n1 > 0 )
 		{
 			sswap_kernel_32(n1, x, y);
 			i=n1;
 		}
+#endif
 
 		while(i < n)
 		{

kernel/power/swap_microk_power10.c

@@ -0,0 +1,105 @@
+/***************************************************************************
+Copyright (c) 2021, 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.
+*****************************************************************************/
+#define HAVE_KERNEL_32 1
+
+#if defined(DOUBLE)
+static void dswap_kernel_32 (long n, double *x, double *y)
+#else
+static void sswap_kernel_32 (long n, float *x, float *y)
+#endif
+{
+  __asm__
+    (
+       ".align	5		\n"
+     "one%=:				\n\t"
+
+       "lxvp		32, 0(%4)	\n\t"
+       "lxvp		34, 32(%4)	\n\t"
+       "lxvp		36, 64(%4)	\n\t"
+       "lxvp		38, 96(%4)	\n\t"
+
+       "lxvp		40, 128(%4)	\n\t"
+       "lxvp		42, 160(%4)	\n\t"
+       "lxvp		44, 192(%4)	\n\t"
+       "lxvp		46, 224(%4)	\n\t"
+
+       "lxvp		48, 0(%3)	\n\t"
+       "lxvp		50, 32(%3)	\n\t"
+       "lxvp		52, 64(%3)	\n\t"
+       "lxvp		54, 96(%3)	\n\t"
+
+       "lxvp		56, 128(%3)	\n\t"
+       "lxvp		58, 160(%3)	\n\t"
+       "lxvp		60, 192(%3)	\n\t"
+       "lxvp		62, 224(%3)	\n\t"
+
+       "stxvp		32, 0(%3)	\n\t"
+       "stxvp		34, 32(%3)	\n\t"
+       "stxvp		36, 64(%3)	\n\t"
+       "stxvp		38, 96(%3)	\n\t"
+
+       "stxvp		40, 128(%3)	\n\t"
+       "stxvp		42, 160(%3)	\n\t"
+       "stxvp		44, 192(%3)	\n\t"
+       "stxvp		46, 224(%3)	\n\t"
+
+       "stxvp		48, 0(%4)	\n\t"
+       "stxvp		50, 32(%4)	\n\t"
+       "stxvp		52, 64(%4)	\n\t"
+       "stxvp		54, 96(%4)	\n\t"
+
+       "stxvp		56, 128(%4)	\n\t"
+       "stxvp		58, 160(%4)	\n\t"
+       "stxvp		60, 192(%4)	\n\t"
+       "stxvp		62, 224(%4)	\n\t"
+
+       "addi		%4, %4, 256	\n\t"
+       "addi		%3, %3, 256	\n\t"
+
+#if defined(DOUBLE)
+       "addic.		%2, %2, -32	\n\t"
+#else
+       "addic.		%2, %2, -64	\n\t"
+#endif
+       "bgt		one%=		\n"
+
+     "#n=%2 x=%0=%3 y=%1=%4"
+     :
+       "+m" (*x),
+       "+m" (*y),
+       "+r" (n),	// 2
+       "+b" (x),	// 3
+       "+b" (y)		// 4
+     :
+     :
+       "cr0",
+       "vs32","vs33","vs34","vs35","vs36","vs37","vs38","vs39",
+       "vs40","vs41","vs42","vs43","vs44","vs45","vs46","vs47",
+       "vs48","vs49","vs50","vs51","vs52","vs53","vs54","vs55",
+       "vs56","vs57","vs58","vs59","vs60","vs61","vs62","vs63"
+     );
+}

kernel/power/zscal.c

@@ -38,11 +38,17 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #pragma GCC optimize "O1"
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #if defined(DOUBLE)
 #include "zscal_microk_power8.c"
 #endif
+#elif defined(POWER10)
+#if defined(DOUBLE)
+#include "zscal_microk_power8.c"
+#else
+#include "cscal_microk_power10.c"
+#endif
 #endif
 #endif
 
@@ -145,7 +151,11 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r,FLOAT da_i, F
 	{
 
 
+#if defined(DOUBLE)
 		n1 = n & -8;
+#else
+		n1 = n & -16;
+#endif
 		if ( n1 > 0 )
 		{
 			zscal_kernel_8(n1, x, da_r, da_i);

kernel/power/zswap.c

@@ -36,9 +36,11 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "common.h"
 
 
-#if defined(POWER8) || defined(POWER9) || defined(POWER10)
 #if defined(__VEC__) || defined(__ALTIVEC__)
+#if defined(POWER8) || defined(POWER9)
 #include "zswap_microk_power8.c"
+#elif defined(POWER10)
+#include "cswap_microk_power10.c"
 #endif
 #endif
 

kernel/x86_64/KERNEL

@@ -489,3 +489,6 @@ XGEMM3MKERNEL    =  xgemm3m_kernel_2x2.S
 
 SSUMKERNEL = ../arm/sum.c
 DSUMKERNEL = ../arm/sum.c
+
+SOMATCOPY_RT = omatcopy_rt.c
+DOMATCOPY_RT = omatcopy_rt.c

kernel/x86_64/KERNEL.ZEN

@@ -97,3 +97,5 @@ ZTRSMKERNEL_RT	= ../generic/trsm_kernel_RT.c
 CGEMM3MKERNEL    =  cgemm3m_kernel_8x4_haswell.c
 ZGEMM3MKERNEL    =  zgemm3m_kernel_4x4_haswell.c
 
+SROTKERNEL = srot.c
+DROTKERNEL = drot.c

kernel/x86_64/dasum.c

@@ -6,7 +6,7 @@
 
 #if defined(SKYLAKEX)
 #include "dasum_microk_skylakex-2.c"
-#elif defined(HASWELL)
+#elif defined(HASWELL) || defined(ZEN)
 #include "dasum_microk_haswell-2.c"
 #endif
 
@@ -93,7 +93,6 @@ FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x)
 #if defined(SMP)
     int nthreads;
     FLOAT dummy_alpha;
-    FLOAT * dummy_b;
 #endif
     FLOAT sumf = 0.0;
 
@@ -115,7 +114,7 @@ FLOAT CNAME(BLASLONG n, FLOAT *x, BLASLONG inc_x)
 #else
         mode = BLAS_DOUBLE | BLAS_REAL;
 #endif
-        blas_level1_thread_with_return_value(mode, n, 0, 0, &dummy_alpha, x, inc_x, dummy_b, 0, result, 0, (void *)asum_thread_function, nthreads);
+        blas_level1_thread_with_return_value(mode, n, 0, 0, &dummy_alpha, x, inc_x, NULL, 0, result, 0, (void *)asum_thread_function, nthreads);
         ptr = (FLOAT *)result;
         for (i = 0; i < nthreads; i++) {
             sumf += (*ptr);

kernel/x86_64/drot.c

@@ -2,7 +2,7 @@
 
 #if defined(SKYLAKEX)
 #include "drot_microk_skylakex-2.c"
-#elif defined(HASWELL)
+#elif defined(HASWELL) || defined(ZEN)
 #include "drot_microk_haswell-2.c"
 #endif
 

kernel/x86_64/omatcopy_rt.c

@@ -0,0 +1,373 @@
+/***************************************************************************
+Copyright (c) 2021, 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 "common.h"
+
+#ifdef HAVE_AVX
+
+#define ROWS_OF_BLOCK 384
+
+  /* +r: %0 = src, %1 = dst, %2 = src_ld, %3 = dst_ld, %4 = dst_tmp */
+/* m: %5 = num_rows, %6 = alpha */
+/* xmm15 = alpha */
+#define TRANS_4x4(a1_no,a2_no,a3_no,a4_no,t1_no,t2_no,t3_no,t4_no)\
+  "vunpcklps %%xmm"#a2_no",%%xmm"#a1_no",%%xmm"#t1_no"; vunpckhps %%xmm"#a2_no",%%xmm"#a1_no",%%xmm"#t2_no";"\
+  "vunpcklps %%xmm"#a4_no",%%xmm"#a3_no",%%xmm"#t3_no"; vunpckhps %%xmm"#a4_no",%%xmm"#a3_no",%%xmm"#t4_no";"\
+  "vunpcklpd %%xmm"#t3_no",%%xmm"#t1_no",%%xmm"#a1_no"; vunpckhpd %%xmm"#t3_no",%%xmm"#t1_no",%%xmm"#a2_no";"\
+  "vunpcklpd %%xmm"#t4_no",%%xmm"#t2_no",%%xmm"#a3_no"; vunpckhpd %%xmm"#t4_no",%%xmm"#t2_no",%%xmm"#a4_no";"
+
+#define TRANS_4x8(a1_no,a2_no,a3_no,a4_no,t1_no,t2_no,t3_no,t4_no)\
+  "vunpcklps %%ymm"#a2_no",%%ymm"#a1_no",%%ymm"#t1_no"; vunpckhps %%ymm"#a2_no",%%ymm"#a1_no",%%ymm"#t2_no";"\
+  "vunpcklps %%ymm"#a4_no",%%ymm"#a3_no",%%ymm"#t3_no"; vunpckhps %%ymm"#a4_no",%%ymm"#a3_no",%%ymm"#t4_no";"\
+  "vunpcklpd %%ymm"#t3_no",%%ymm"#t1_no",%%ymm"#a1_no"; vunpckhpd %%ymm"#t3_no",%%ymm"#t1_no",%%ymm"#a2_no";"\
+  "vunpcklpd %%ymm"#t4_no",%%ymm"#t2_no",%%ymm"#a3_no"; vunpckhpd %%ymm"#t4_no",%%ymm"#t2_no",%%ymm"#a4_no";"
+
+#define SAVE_4x4(b1_no,b2_no,b3_no,b4_no)\
+  "vmovups %%xmm"#b1_no",(%4); vmovups %%xmm"#b2_no",(%4,%3,1); leaq (%4,%3,2),%4;"\
+  "vmovups %%xmm"#b3_no",(%4); vmovups %%xmm"#b4_no",(%4,%3,1); leaq (%4,%3,2),%4;"
+
+#define SAVE_4x8(b1_no,b2_no,b3_no,b4_no) SAVE_4x4(b1_no,b2_no,b3_no,b4_no)\
+  "vextractf128 $1,%%ymm"#b1_no",(%4); vextractf128 $1,%%ymm"#b2_no",(%4,%3,1); leaq (%4,%3,2),%4;"\
+  "vextractf128 $1,%%ymm"#b3_no",(%4); vextractf128 $1,%%ymm"#b4_no",(%4,%3,1); leaq (%4,%3,2),%4;"
+
+#define COPY_4x16 "movq %1,%4; addq $16,%1;"\
+  "vmulps (%0),%%ymm15,%%ymm0; vmulps 32(%0),%%ymm15,%%ymm4; vmulps (%0,%2,1),%%ymm15,%%ymm1; vmulps 32(%0,%2,1),%%ymm15,%%ymm5; leaq (%0,%2,2),%0;"\
+  "vmulps (%0),%%ymm15,%%ymm2; vmulps 32(%0),%%ymm15,%%ymm6; vmulps (%0,%2,1),%%ymm15,%%ymm3; vmulps 32(%0,%2,1),%%ymm15,%%ymm7; leaq (%0,%2,2),%0;"\
+  TRANS_4x8(0,1,2,3,8,9,10,11) SAVE_4x8(0,1,2,3)\
+  TRANS_4x8(4,5,6,7,8,9,10,11) SAVE_4x8(4,5,6,7)
+
+#define COPY_4x8 "movq %1,%4; addq $16,%1;"\
+  "vmulps (%0),%%ymm15,%%ymm0; vmulps (%0,%2,1),%%ymm15,%%ymm1; leaq (%0,%2,2),%0;"\
+  "vmulps (%0),%%ymm15,%%ymm2; vmulps (%0,%2,1),%%ymm15,%%ymm3; leaq (%0,%2,2),%0;"\
+  TRANS_4x8(0,1,2,3,8,9,10,11) SAVE_4x8(0,1,2,3)
+
+#define COPY_4x4 "movq %1,%4; addq $16,%1;"\
+  "vmulps (%0),%%xmm15,%%xmm0; vmulps (%0,%2,1),%%xmm15,%%xmm1; leaq (%0,%2,2),%0;"\
+  "vmulps (%0),%%xmm15,%%xmm2; vmulps (%0,%2,1),%%xmm15,%%xmm3; leaq (%0,%2,2),%0;"\
+  TRANS_4x4(0,1,2,3,8,9,10,11) SAVE_4x4(0,1,2,3)
+
+#define COPY_4x2 \
+  "vmovsd (%0),%%xmm0; vmovhpd (%0,%2,1),%%xmm0,%%xmm0; vmulps %%xmm15,%%xmm0,%%xmm0; leaq (%0,%2,2),%0;"\
+  "vmovsd (%0),%%xmm1; vmovhpd (%0,%2,1),%%xmm1,%%xmm1; vmulps %%xmm15,%%xmm1,%%xmm1; leaq (%0,%2,2),%0;"\
+  "vpermilps $216,%%xmm0,%%xmm0; vpermilps $216,%%xmm1,%%xmm1; vunpcklpd %%xmm1,%%xmm0,%%xmm2; vunpckhpd %%xmm1,%%xmm0,%%xmm3;"\
+  "vmovups %%xmm2,(%1); vmovups %%xmm3,(%1,%3,1); addq $16,%1;"
+
+#define COPY_4x1 \
+  "vmovss (%0),%%xmm0; vinsertps $16,(%0,%2,1),%%xmm0,%%xmm0; leaq (%0,%2,2),%0;"\
+  "vinsertps $32,(%0),%%xmm0,%%xmm0; vinsertps $48,(%0,%2,1),%%xmm0,%%xmm0; leaq (%0,%2,2),%0;"\
+  "vmulps %%xmm15,%%xmm0,%%xmm0; vmovups %%xmm0,(%1); addq $16,%1;"
+
+#define SAVE_2x4(c1_no,c2_no,t1_no,t2_no) \
+  "vunpcklps %%xmm"#c2_no",%%xmm"#c1_no",%%xmm"#t1_no"; vmulps %%xmm15,%%xmm"#t1_no",%%xmm"#t1_no";"\
+  "vmovsd %%xmm"#t1_no",(%4); vmovhpd %%xmm"#t1_no",(%4,%3,1); leaq (%4,%3,2),%4;"\
+  "vunpckhps %%xmm"#c2_no",%%xmm"#c1_no",%%xmm"#t2_no"; vmulps %%xmm15,%%xmm"#t2_no",%%xmm"#t2_no";"\
+  "vmovsd %%xmm"#t2_no",(%4); vmovhpd %%xmm"#t2_no",(%4,%3,1); leaq (%4,%3,2),%4;"
+
+#define COPY_2x16 "movq %1,%4; addq $8,%1;"\
+  "vmovups (%0),%%ymm0; vmovups 32(%0),%%ymm2; vmovups (%0,%2,1),%%ymm1; vmovups 32(%0,%2,1),%%ymm3; leaq (%0,%2,2),%0;"\
+  "vextractf128 $1,%%ymm0,%%xmm4; vextractf128 $1,%%ymm2,%%xmm6; vextractf128 $1,%%ymm1,%%xmm5; vextractf128 $1,%%ymm3,%%xmm7;"\
+  SAVE_2x4(0,1,8,9) SAVE_2x4(4,5,8,9) SAVE_2x4(2,3,8,9) SAVE_2x4(6,7,8,9)
+
+#define COPY_2x8 "movq %1,%4; addq $8,%1;"\
+  "vmovups (%0),%%ymm0; vmovups (%0,%2,1),%%ymm1; leaq (%0,%2,2),%0;"\
+  "vextractf128 $1,%%ymm0,%%xmm2; vextractf128 $1,%%ymm1,%%xmm3;"\
+  SAVE_2x4(0,1,4,5) SAVE_2x4(2,3,4,5)
+
+#define COPY_2x4 "movq %1,%4; addq $8,%1;"\
+  "vmovups (%0),%%xmm0; vmovups (%0,%2,1),%%xmm1; leaq (%0,%2,2),%0;"\
+  SAVE_2x4(0,1,4,5)
+
+#define COPY_2x2 \
+  "vmovsd (%0),%%xmm0; vmovhpd (%0,%2,1),%%xmm0,%%xmm0; vmulps %%xmm15,%%xmm0,%%xmm0; leaq (%0,%2,2),%0; vpermilps $216,%%xmm0,%%xmm0;"\
+  "vmovsd %%xmm0,(%1); vmovhpd %%xmm0,(%1,%3,1); addq $8,%1;"
+
+#define COPY_2x1 \
+  "vmovss (%0),%%xmm0; vinsertps $16,(%0,%2,1),%%xmm0,%%xmm0; vmulps %%xmm15,%%xmm0,%%xmm0; leaq (%0,%2,2),%0; vmovsd %%xmm0,(%1); addq $8,%1;"
+
+#define SAVE_1x4(c1_no)\
+  "vmulps %%xmm15,%%xmm"#c1_no",%%xmm"#c1_no"; vmovss %%xmm"#c1_no",(%4); vextractps $1,%%xmm"#c1_no",(%4,%3,1); leaq (%4,%3,2),%4;"\
+  "vextractps $2,%%xmm"#c1_no",(%4); vextractps $3,%%xmm"#c1_no",(%4,%3,1); leaq (%4,%3,2),%4;"
+
+#define COPY_1x16 "movq %1,%4; addq $4,%1;"\
+  "vmovups (%0),%%xmm1;" SAVE_1x4(1) "vmovups 16(%0),%%xmm2;" SAVE_1x4(2)\
+  "vmovups 32(%0),%%xmm1;" SAVE_1x4(1) "vmovups 48(%0),%%xmm2;" SAVE_1x4(2) "addq %2,%0;"
+
+#define COPY_1x8 "movq %1,%4; addq $4,%1;"\
+  "vmovups (%0),%%xmm1;" SAVE_1x4(1) "vmovups 16(%0),%%xmm2;" SAVE_1x4(2) "addq %2,%0;"
+
+#define COPY_1x4 "movq %1,%4; addq $4,%1; vmovups (%0),%%xmm1;" SAVE_1x4(1) "addq %2,%0;"
+
+#define COPY_1x2 "vmovsd (%0),%%xmm1; addq %2,%0; vmulps %%xmm15,%%xmm1,%%xmm1; vmovss %%xmm1,(%1); vextractps $1,%%xmm1,(%1,%3,1); addq $4,%1;"
+
+#define COPY_1x1 "vmulss (%0),%%xmm15,%%xmm1; vmovss %%xmm1,(%1); addq %2,%0; addq $4,%1;"
+
+#define COMPUTE(ndim){\
+  src = src_base; dst = dst_base;\
+  __asm__ __volatile__(\
+    "vbroadcastss %6,%%ymm15; movq %5,%%r11; cmpq $4,%%r11; jb "#ndim"32f;"\
+    #ndim"31:\n\t"\
+    COPY_4x##ndim "subq $4,%%r11; cmpq $4,%%r11; jnb "#ndim"31b;"\
+    #ndim"32:\n\t"\
+    "cmpq $2,%%r11; jb "#ndim"33f;"\
+    COPY_2x##ndim "subq $2,%%r11;"\
+    #ndim"33:\n\t"\
+    "testq %%r11,%%r11; jz "#ndim"34f;"\
+    COPY_1x##ndim "subq $1,%%r11;"\
+    #ndim"34:\n\t"\
+    :"+r"(src),"+r"(dst),"+r"(src_ld_bytes),"+r"(dst_ld_bytes),"+r"(dst_tmp):"m"(num_rows),"m"(ALPHA):"r11","cc","memory"\
+    ,"xmm0","xmm1","xmm2","xmm3","xmm4","xmm5","xmm6","xmm7","xmm8","xmm9","xmm10","xmm11","xmm12","xmm13","xmm14","xmm15");\
+}
+int CNAME(BLASLONG rows, BLASLONG cols, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *b, BLASLONG ldb){
+  float *src, *dst, *dst_tmp, *src_base, *dst_base;
+  uint64_t src_ld_bytes = (uint64_t)lda * sizeof(float), dst_ld_bytes = (uint64_t)ldb * sizeof(float), num_rows = 0;
+  BLASLONG cols_left, rows_done; float ALPHA = alpha;
+  if(ALPHA==0.0){
+    dst_base = b;
+    for(cols_left=cols;cols_left>0;cols_left--) {memset(dst_base,0,rows*sizeof(float)); dst_base += ldb;}
+    return 0;
+  }
+  for(rows_done=0;rows_done<rows;rows_done+=num_rows){
+    num_rows = rows-rows_done;
+    if(num_rows > ROWS_OF_BLOCK) num_rows = ROWS_OF_BLOCK;
+    cols_left = cols; src_base = a + (int64_t)lda * (int64_t)rows_done; dst_base = b + rows_done;
+    if(ldb%1024>3 && ldb%1024<1021) for(;cols_left>15;cols_left-=16){COMPUTE(16) src_base += 16; dst_base += 16 * ldb;}
+    for(;cols_left>7;cols_left-=8){COMPUTE(8) src_base += 8; dst_base += 8 * ldb;}
+    for(;cols_left>3;cols_left-=4){COMPUTE(4) src_base += 4; dst_base += 4 * ldb;}
+    for(;cols_left>1;cols_left-=2){COMPUTE(2) src_base += 2; dst_base += 2 * ldb;}
+    if(cols_left>0){COMPUTE(1) src_base ++; dst_base += ldb;}
+  }
+  return 0;
+}
+
+#else
+
+int CNAME(BLASLONG rows, BLASLONG cols, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *b, BLASLONG ldb)
+{
+    BLASLONG i, j;
+    FLOAT *a_offset, *a_offset1, *a_offset2, *a_offset3, *a_offset4;
+    FLOAT *b_offset, *b_offset1, *b_offset2, *b_offset3, *b_offset4;
+
+    if (rows <= 0)  return 0;
+    if (cols <= 0)  return 0;
+
+    a_offset = a;
+    b_offset = b;
+
+    i = (rows >> 2);
+    if (i > 0) {
+        do {
+            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_offset2 = b_offset1 + ldb;
+            b_offset3 = b_offset2 + ldb;
+            b_offset4 = b_offset3 + ldb;
+            b_offset += 4;
+
+            j = (cols >> 2);
+            if (j > 0) {
+                do {
+                /* Column 1 of MAT_B */
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha; // Row 1 of MAT_A
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+                *(b_offset3 + 0) = *(a_offset1 + 2)*alpha;
+                *(b_offset4 + 0) = *(a_offset1 + 3)*alpha;
+
+                /* Column 2 of MAT_B */
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha; // Row 2 of MAT_A
+                *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+                *(b_offset3 + 1) = *(a_offset2 + 2)*alpha;
+                *(b_offset4 + 1) = *(a_offset2 + 3)*alpha;
+
+                /* Column 3 of MAT_B */
+                *(b_offset1 + 2) = *(a_offset3 + 0)*alpha; // Row 3 of MAT_A
+                *(b_offset2 + 2) = *(a_offset3 + 1)*alpha;
+                *(b_offset3 + 2) = *(a_offset3 + 2)*alpha;
+                *(b_offset4 + 2) = *(a_offset3 + 3)*alpha;
+
+                /* Column 4 of MAT_B */
+                *(b_offset1 + 3) = *(a_offset4 + 0)*alpha; // Row 4 of MAT_A
+                *(b_offset2 + 3) = *(a_offset4 + 1)*alpha;
+                *(b_offset3 + 3) = *(a_offset4 + 2)*alpha;
+                *(b_offset4 + 3) = *(a_offset4 + 3)*alpha;
+
+                a_offset1 += 4;
+                a_offset2 += 4;
+                a_offset3 += 4;
+                a_offset4 += 4;
+                b_offset1 += ldb * 4;
+                b_offset2 += ldb * 4;
+                b_offset3 += ldb * 4;
+                b_offset4 += ldb * 4;
+                
+                j--;
+                } while (j > 0);
+            } // if(j > 0)
+
+
+            if (cols & 2) {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+                *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+
+                *(b_offset1 + 2) = *(a_offset3 + 0)*alpha;
+                *(b_offset2 + 2) = *(a_offset3 + 1)*alpha;
+
+                *(b_offset1 + 3) = *(a_offset4 + 0)*alpha;
+                *(b_offset2 + 3) = *(a_offset4 + 1)*alpha;
+
+                a_offset1 += 2;
+                a_offset2 += 2;
+                a_offset3 += 2;
+                a_offset4 += 2;
+
+                b_offset1 += ldb*2;
+
+            }
+
+            if (cols & 1) {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+
+                *(b_offset1 + 2) = *(a_offset3 + 0)*alpha;
+
+                *(b_offset1 + 3) = *(a_offset4 + 0)*alpha;
+            }
+
+            i--;
+        } while (i > 0);
+    }
+
+
+    if (rows & 2) {
+        a_offset1 = a_offset;
+        a_offset2 = a_offset1 + lda;
+        a_offset += 2 * lda;
+
+        b_offset1 = b_offset;
+        b_offset2 = b_offset1 + ldb;
+        b_offset3 = b_offset2 + ldb;
+        b_offset4 = b_offset3 + ldb;
+        b_offset += 2;
+
+        j = (cols >> 2);
+        if (j > 0){
+            do {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+                *(b_offset3 + 0) = *(a_offset1 + 2)*alpha;
+                *(b_offset4 + 0) = *(a_offset1 + 3)*alpha;
+
+                *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+                *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+                *(b_offset3 + 1) = *(a_offset2 + 2)*alpha;
+                *(b_offset4 + 1) = *(a_offset2 + 3)*alpha;
+                
+                a_offset1 += 4;
+                a_offset2 += 4;
+                b_offset1 += ldb * 4;
+                b_offset2 += ldb * 4;
+                b_offset3 += ldb * 4;
+                b_offset4 += ldb * 4;
+
+                j--;
+            } while (j > 0);
+        }
+
+
+        if (cols & 2){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+            *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+
+            *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+            *(b_offset2 + 1) = *(a_offset2 + 1)*alpha;
+
+            a_offset1 += 2;
+            a_offset2 += 2;
+            b_offset1 += ldb*2;
+
+        }
+
+
+        if (cols & 1){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+            *(b_offset1 + 1) = *(a_offset2 + 0)*alpha;
+        }
+    } // if (rows & 2)
+
+
+    if (rows & 1) {
+        a_offset1 = a_offset;
+        a_offset += lda;
+
+        b_offset1 = b_offset;
+        b_offset2 = b_offset1 + ldb;
+        b_offset3 = b_offset2 + ldb;
+        b_offset4 = b_offset3 + ldb;
+
+        j = (cols >> 2);
+        if (j > 0){
+            do {
+                *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+                *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+                *(b_offset3 + 0) = *(a_offset1 + 2)*alpha;
+                *(b_offset4 + 0) = *(a_offset1 + 3)*alpha;
+                
+                a_offset1 += 4;
+                b_offset1 += ldb * 4;
+                b_offset2 += ldb * 4;
+                b_offset3 += ldb * 4;
+                b_offset4 += ldb * 4;
+                
+                j--;
+            } while (j > 0);
+        }
+
+        if (cols & 2){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+            *(b_offset2 + 0) = *(a_offset1 + 1)*alpha;
+
+            a_offset1 += 2;
+            b_offset1 += ldb * 2;
+        }
+
+        if (cols & 1){
+            *(b_offset1 + 0) = *(a_offset1 + 0)*alpha;
+        }
+    }
+
+    return 0;
+}
+
+#endif

kernel/x86_64/sasum.c

@@ -11,7 +11,7 @@
 
 #if defined(SKYLAKEX)
 #include "sasum_microk_skylakex-2.c"
-#elif defined(HASWELL)
+#elif defined(HASWELL) || defined(ZEN)
 #include "sasum_microk_haswell-2.c"
 #endif
 

kernel/x86_64/sbgemm_block_microk_cooperlake.c

@@ -0,0 +1,426 @@
+#include "sbgemm.h"
+
+#include <immintrin.h>
+// Walk around those intrinsics that missed by compiler
+#define MM256_LOADU_EPI16(addr)   \
+            _mm256_maskz_loadu_epi16(~0, (addr))
+#define MM256_STOREU_EPI16(addr, reg)  \
+            _mm256_mask_storeu_epi16((addr), ~0, (reg))
+
+#include <stdio.h>
+void print_block(BLASLONG m, BLASLONG n, bfloat16 * mat)
+{
+    printf("---- BLOCK %ld x %ld ----\n", m, n);
+    for (BLASLONG i=0; i<m; i++) {
+        for (BLASLONG j=0; j<n; j++) {
+            printf("%-4X  ", *(mat + i*n +j));
+        }
+        printf("\n");
+    }
+    printf("---- End of BLOCK ----\n");
+}
+
+void COL_MAJOR_INCOPY_KERNEL_Kx32(BLASLONG k, bfloat16 * A, BLASLONG lda, bfloat16 * block_A)
+{
+    BLASLONG tag_k_2x = k & (~1);
+
+    __m512i array512_0, array512_1, array512_2, array512_3;
+
+    BLASLONG idx_src_base0, idx_src_base1;
+    BLASLONG idx_target_base0, idx_target_base1;
+
+    BLASLONG LDA_2x = 2*lda;
+    BLASLONG BF16_BLOCK_T_M_2x = 2*32;
+    idx_src_base0 = 0;
+    idx_src_base1 = lda;
+    idx_target_base0 = 0;
+    idx_target_base1 = 32;
+    for (BLASLONG idx_k = 0; idx_k < tag_k_2x; idx_k += 2) {
+        array512_0 = _mm512_loadu_si512(&A[idx_src_base0]);
+        array512_1 = _mm512_loadu_si512(&A[idx_src_base1]);
+        array512_2 = _mm512_unpacklo_epi16(array512_0, array512_1);
+        array512_3 = _mm512_unpackhi_epi16(array512_0, array512_1);
+        _mm512_storeu_si512(&block_A[idx_target_base0], array512_2);
+        _mm512_storeu_si512(&block_A[idx_target_base1], array512_3);
+
+        idx_src_base0 += LDA_2x;
+        idx_src_base1 += LDA_2x;
+        idx_target_base0 += BF16_BLOCK_T_M_2x;
+        idx_target_base1 += BF16_BLOCK_T_M_2x;
+    }
+
+    if (tag_k_2x != k) {
+        __m512i ZERO512 = _mm512_setzero_si512();
+        array512_0 = _mm512_loadu_si512(&A[idx_src_base0]);
+        array512_2 = _mm512_unpacklo_epi16(array512_0, ZERO512);
+        array512_3 = _mm512_unpackhi_epi16(array512_0, ZERO512);
+        _mm512_storeu_si512(&block_A[idx_target_base0], array512_2);
+        _mm512_storeu_si512(&block_A[idx_target_base1], array512_3);
+   }
+
+#ifdef DEBUG_PROFILE
+   print_block(BF16_BLOCK_THRES_K, BF16_BLOCK_THRES_M, block_A);
+#endif
+}
+
+void COL_MAJOR_INCOPY_KERNEL_Kx32m(BLASLONG k, BLASLONG m, bfloat16 * A, BLASLONG lda, bfloat16 * block_A)
+{
+    BLASLONG tag_k_2x = k & (~1);
+    unsigned int tail_mask_value = (((unsigned int)0xffffffff) >> (32-m));
+    __mmask32 tail_mask = *((__mmask32*) &tail_mask_value);
+
+    __m512i array512_0, array512_1, array512_2, array512_3;
+
+    BLASLONG idx_src_base0, idx_src_base1;
+    BLASLONG idx_target_base0, idx_target_base1;
+
+    BLASLONG LDA_2x = 2*lda;
+    BLASLONG BF16_BLOCK_T_M_2x = 2*32;
+    idx_src_base0 = 0;
+    idx_src_base1 = lda;
+    idx_target_base0 = 0;
+    idx_target_base1 = 32;
+    for (BLASLONG idx_k = 0; idx_k < tag_k_2x; idx_k += 2) {
+        array512_0 = _mm512_maskz_loadu_epi16(tail_mask, &A[idx_src_base0]);
+        array512_1 = _mm512_maskz_loadu_epi16(tail_mask, &A[idx_src_base1]);
+        array512_2 = _mm512_unpacklo_epi16(array512_0, array512_1);
+        array512_3 = _mm512_unpackhi_epi16(array512_0, array512_1);
+        _mm512_storeu_si512(&block_A[idx_target_base0], array512_2);
+        _mm512_storeu_si512(&block_A[idx_target_base1], array512_3);
+
+        idx_src_base0 += LDA_2x;
+        idx_src_base1 += LDA_2x;
+        idx_target_base0 += BF16_BLOCK_T_M_2x;
+        idx_target_base1 += BF16_BLOCK_T_M_2x;
+    }
+
+    if (tag_k_2x != k) {
+        __m512i ZERO512 = _mm512_setzero_si512();
+        array512_0 = _mm512_maskz_loadu_epi16(tail_mask, &A[idx_src_base0]);
+        array512_2 = _mm512_unpacklo_epi16(array512_0, ZERO512);
+        array512_3 = _mm512_unpackhi_epi16(array512_0, ZERO512);
+        _mm512_storeu_si512(&block_A[idx_target_base0], array512_2);
+        _mm512_storeu_si512(&block_A[idx_target_base1], array512_3);
+    }
+
+#ifdef DEBUG_PROFILE
+   print_block(BF16_BLOCK_THRES_K, BF16_BLOCK_THRES_M, block_A);
+#endif
+}
+
+void COL_MAJOR_INCOPY_KERNEL_Kx16(BLASLONG k, BLASLONG m, bfloat16 * A, BLASLONG lda, bfloat16 * block_A)
+{
+    BLASLONG tag_k_2x = k & (~1);
+
+    __m256i array256_0, array256_1, array256_2, array256_3;
+
+    BLASLONG idx_src_base0, idx_src_base1;
+    BLASLONG idx_target_base0;
+
+    BLASLONG LDA_2x = 2*lda;
+    idx_src_base0 = 0;
+    idx_src_base1 = lda;
+    idx_target_base0 = 0;
+    for (BLASLONG idx_k = 0; idx_k < tag_k_2x; idx_k += 2) {
+        array256_0 = MM256_LOADU_EPI16(&A[idx_src_base0]);
+        array256_1 = MM256_LOADU_EPI16(&A[idx_src_base1]);
+        array256_2 = _mm256_unpacklo_epi16(array256_0, array256_1);
+        array256_3 = _mm256_unpackhi_epi16(array256_0, array256_1);
+        // Store in one row of block_B
+        MM256_STOREU_EPI16(&block_A[idx_target_base0],      array256_2);
+        MM256_STOREU_EPI16(&block_A[idx_target_base0 + 16], array256_3);
+
+        idx_src_base0 += LDA_2x;
+        idx_src_base1 += LDA_2x;
+        idx_target_base0 += 32;
+    }
+
+    if (tag_k_2x != k) {
+        __m256i ZERO256 = _mm256_setzero_si256();
+        array256_0 = MM256_LOADU_EPI16(&A[idx_src_base0]);
+        array256_2 = _mm256_unpacklo_epi16(array256_0, ZERO256);
+        array256_3 = _mm256_unpackhi_epi16(array256_0, ZERO256);
+        // Store in one row of block_B
+        MM256_STOREU_EPI16(&block_A[idx_target_base0],      array256_2);
+        MM256_STOREU_EPI16(&block_A[idx_target_base0 + 16], array256_3);
+    }
+
+#ifdef DEBUG_PROFILE
+   print_block(BF16_BLOCK_THRES_K, BF16_BLOCK_THRES_M, block_A);
+#endif
+}
+
+void COL_MAJOR_INCOPY_KERNEL_Kx16m(BLASLONG k, BLASLONG m, bfloat16 * A, BLASLONG lda, bfloat16 * block_A)
+{
+    BLASLONG tag_k_2x = k & (~1);
+    unsigned short tail_mask_value = (((unsigned short)0xffff) >> (16-m));
+    __mmask16 tail_mask = *((__mmask16*) &tail_mask_value);
+
+    __m256i array256_0, array256_1, array256_2, array256_3;
+
+    BLASLONG idx_src_base0, idx_src_base1;
+    BLASLONG idx_target_base0;
+
+    BLASLONG LDA_2x = 2*lda;
+    idx_src_base0 = 0;
+    idx_src_base1 = lda;
+    idx_target_base0 = 0;
+    for (BLASLONG idx_k = 0; idx_k < tag_k_2x; idx_k += 2) {
+        array256_0 = _mm256_maskz_loadu_epi16(tail_mask, &A[idx_src_base0]);
+        array256_1 = _mm256_maskz_loadu_epi16(tail_mask, &A[idx_src_base1]);
+        array256_2 = _mm256_unpacklo_epi16(array256_0, array256_1);
+        array256_3 = _mm256_unpackhi_epi16(array256_0, array256_1);
+        // Store in one row of block_B
+        MM256_STOREU_EPI16(&block_A[idx_target_base0],      array256_2);
+        MM256_STOREU_EPI16(&block_A[idx_target_base0 + 16], array256_3);
+
+        idx_src_base0 += LDA_2x;
+        idx_src_base1 += LDA_2x;
+        idx_target_base0 += 32;
+    }
+
+    if (tag_k_2x != k) {
+        __m256i ZERO256 = _mm256_setzero_si256();
+        array256_0 = _mm256_maskz_loadu_epi16(tail_mask, &A[idx_src_base0]);
+        array256_2 = _mm256_unpacklo_epi16(array256_0, ZERO256);
+        array256_3 = _mm256_unpackhi_epi16(array256_0, ZERO256);
+        // Store in one row of block_B
+        MM256_STOREU_EPI16(&block_A[idx_target_base0],      array256_2);
+        MM256_STOREU_EPI16(&block_A[idx_target_base0 + 16], array256_3);
+    }
+
+#ifdef DEBUG_PROFILE
+   print_block(BF16_BLOCK_THRES_K, BF16_BLOCK_THRES_M, block_A);
+#endif
+}
+
+void COL_MAJOR_ONCOPY_KERNEL_8x32(BLASLONG k, bfloat16 * B, BLASLONG ldb, bfloat16 * block_B)
+{
+    BLASLONG tag_k_32x = k & (~31);
+    BLASLONG idx_src_base0, idx_src_base1, idx_src_base2, idx_src_base3, idx_src_base4, idx_src_base5, idx_src_base6, idx_src_base7;
+    BLASLONG idx_target_base0;
+
+    idx_src_base0 = 0;
+    idx_src_base1 = 1*ldb;
+    idx_src_base2 = 2*ldb;
+    idx_src_base3 = 3*ldb;
+    idx_src_base4 = 4*ldb;
+    idx_src_base5 = 5*ldb;
+    idx_src_base6 = 6*ldb;
+    idx_src_base7 = 7*ldb;
+    idx_target_base0 = 0;
+
+    for (BLASLONG idx_k = 0; idx_k < tag_k_32x; idx_k += 32) {
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*0], _mm512_loadu_si512(&B[idx_src_base0+idx_k]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*1], _mm512_loadu_si512(&B[idx_src_base1+idx_k]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*2], _mm512_loadu_si512(&B[idx_src_base2+idx_k]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*3], _mm512_loadu_si512(&B[idx_src_base3+idx_k]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*4], _mm512_loadu_si512(&B[idx_src_base4+idx_k]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*5], _mm512_loadu_si512(&B[idx_src_base5+idx_k]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*6], _mm512_loadu_si512(&B[idx_src_base6+idx_k]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*7], _mm512_loadu_si512(&B[idx_src_base7+idx_k]));
+        idx_target_base0 += 32*8;
+    }
+
+    if (tag_k_32x != k) {
+        unsigned int tail_mask_value = (((unsigned int)0xffffffff) >> (32-(k-tag_k_32x)));
+        __mmask32 tail_mask = *((__mmask32*) &tail_mask_value);
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*0], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base0+tag_k_32x]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*1], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base1+tag_k_32x]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*2], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base2+tag_k_32x]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*3], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base3+tag_k_32x]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*4], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base4+tag_k_32x]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*5], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base5+tag_k_32x]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*6], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base6+tag_k_32x]));
+        _mm512_storeu_si512(&block_B[idx_target_base0+ 32*7], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base7+tag_k_32x]));
+    }
+
+#ifdef DEBUG_PROFILE
+   print_block(BF16_BLOCK_THRES_N, BF16_BLOCK_THRES_K, block_B);
+#endif
+}
+
+void COL_MAJOR_ONCOPY_KERNEL_Nx32(BLASLONG n, BLASLONG k, bfloat16 * B, BLASLONG ldb, bfloat16 * block_B)
+{
+    BLASLONG tag_k_32x = k & (~31);
+    BLASLONG tag_n_2x  = n & (~1);
+    BLASLONG idx_src_base0;
+    BLASLONG idx_target_base0;
+
+    BLASLONG LDB_2x = 2*ldb;
+
+    idx_target_base0 = 0;
+
+    for (BLASLONG idx_k = 0; idx_k < tag_k_32x; idx_k += 32) {
+        idx_src_base0 = 0;
+        for (BLASLONG idx_n = 0; idx_n < tag_n_2x; idx_n += 2) {
+            _mm512_storeu_si512(&block_B[idx_target_base0+ 32*0], _mm512_loadu_si512(&B[idx_src_base0       + idx_k]));
+            _mm512_storeu_si512(&block_B[idx_target_base0+ 32*1], _mm512_loadu_si512(&B[idx_src_base0 + ldb + idx_k]));
+            idx_src_base0 += LDB_2x;
+            idx_target_base0 += 64;
+        }
+        
+        if (tag_n_2x != n) {
+            _mm512_storeu_si512(&block_B[idx_target_base0],  _mm512_loadu_si512(&B[idx_src_base0       + idx_k]));
+            idx_target_base0 += 32;
+        }
+    }
+
+    if (tag_k_32x != k) {
+        unsigned int tail_mask_value = (((unsigned int)0xffffffff) >> (32-(k-tag_k_32x)));
+        __mmask32 tail_mask = *((__mmask32*) &tail_mask_value);
+        idx_src_base0 = 0;
+        for (BLASLONG idx_n = 0; idx_n < tag_n_2x; idx_n += 2) {
+            _mm512_storeu_si512(&block_B[idx_target_base0+ 32*0], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base0       + tag_k_32x]));
+            _mm512_storeu_si512(&block_B[idx_target_base0+ 32*1], _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base0 + ldb + tag_k_32x]));
+            idx_src_base0 += LDB_2x;
+            idx_target_base0 += 64;
+        }
+        
+        if (tag_n_2x != n) {
+            _mm512_storeu_si512(&block_B[idx_target_base0],  _mm512_maskz_loadu_epi16(tail_mask, &B[idx_src_base0       + tag_k_32x]));
+        }
+    }
+
+#ifdef DEBUG_PROFILE
+   print_block(BF16_BLOCK_THRES_N, BF16_BLOCK_THRES_K, block_B);
+#endif
+}
+
+// Scale matrix C while beta is not ZERO or ONE
+void sbgemm_scal_operation(OPENBLAS_CONST enum CBLAS_ORDER Order, OPENBLAS_CONST blasint M, OPENBLAS_CONST blasint N, OPENBLAS_CONST float beta, float *C, OPENBLAS_CONST blasint ldc)
+{
+    BLASLONG tag_n_Nx = N & (~3);
+    BLASLONG tag_n_Mx = M & (~15);
+
+    BLASLONG LDC4x = ldc*4;
+    BLASLONG idx_base_0 = 0;
+    BLASLONG idx_base_1 = ldc;
+    BLASLONG idx_base_2 = ldc*2;
+    BLASLONG idx_base_3 = ldc*3;
+
+    unsigned short tail_mask_value = (((unsigned short)0xffff) >> (16-M+tag_n_Mx));
+    __mmask16 tail_mask = *((__mmask16*) &tail_mask_value);
+
+    __m512 array_512_0, array_512_1, array_512_2, array_512_3;
+
+    __m512  BETAVECTOR  = _mm512_set1_ps(beta);
+
+    if (Order == CblasColMajor) {
+        for (BLASLONG idx_n = 0; idx_n < tag_n_Nx; idx_n += 4) {
+            for (BLASLONG idx_m = 0; idx_m < tag_n_Mx; idx_m += 16) {
+                array_512_0 = _mm512_loadu_ps(&C[idx_base_0+idx_m]);
+                array_512_1 = _mm512_loadu_ps(&C[idx_base_1+idx_m]);
+                array_512_2 = _mm512_loadu_ps(&C[idx_base_2+idx_m]);
+                array_512_3 = _mm512_loadu_ps(&C[idx_base_3+idx_m]);
+
+                array_512_0 = _mm512_mul_ps(BETAVECTOR, array_512_0);
+                array_512_1 = _mm512_mul_ps(BETAVECTOR, array_512_1);
+                array_512_2 = _mm512_mul_ps(BETAVECTOR, array_512_2);
+                array_512_3 = _mm512_mul_ps(BETAVECTOR, array_512_3);
+
+                _mm512_storeu_ps(&C[idx_base_0+idx_m], array_512_0);
+                _mm512_storeu_ps(&C[idx_base_1+idx_m], array_512_1);
+                _mm512_storeu_ps(&C[idx_base_2+idx_m], array_512_2);
+                _mm512_storeu_ps(&C[idx_base_3+idx_m], array_512_3);
+            }
+
+            if (tag_n_Mx != M) {
+                array_512_0 = _mm512_maskz_loadu_ps(tail_mask, &C[idx_base_0+tag_n_Mx]);
+                array_512_1 = _mm512_maskz_loadu_ps(tail_mask, &C[idx_base_1+tag_n_Mx]);
+                array_512_2 = _mm512_maskz_loadu_ps(tail_mask, &C[idx_base_2+tag_n_Mx]);
+                array_512_3 = _mm512_maskz_loadu_ps(tail_mask, &C[idx_base_3+tag_n_Mx]);
+
+                array_512_0 = _mm512_mul_ps(BETAVECTOR, array_512_0);
+                array_512_1 = _mm512_mul_ps(BETAVECTOR, array_512_1);
+                array_512_2 = _mm512_mul_ps(BETAVECTOR, array_512_2);
+                array_512_3 = _mm512_mul_ps(BETAVECTOR, array_512_3);
+
+                _mm512_mask_storeu_ps(&C[idx_base_0+tag_n_Mx], tail_mask, array_512_0);
+                _mm512_mask_storeu_ps(&C[idx_base_1+tag_n_Mx], tail_mask, array_512_1);
+                _mm512_mask_storeu_ps(&C[idx_base_2+tag_n_Mx], tail_mask, array_512_2);
+                _mm512_mask_storeu_ps(&C[idx_base_3+tag_n_Mx], tail_mask, array_512_3);
+            }
+
+            idx_base_0 += LDC4x;
+            idx_base_1 += LDC4x;
+            idx_base_2 += LDC4x;
+            idx_base_3 += LDC4x;
+        }
+
+        if (tag_n_Nx != N) {
+            for (BLASLONG idx_n = tag_n_Nx; idx_n < N; idx_n++) {
+                for (BLASLONG idx_m = 0; idx_m < tag_n_Mx; idx_m += 16) {
+                    array_512_0 = _mm512_loadu_ps(&C[idx_base_0+idx_m]);
+                    array_512_0 = _mm512_mul_ps(BETAVECTOR, array_512_0);
+                    _mm512_storeu_ps(&C[idx_base_0+idx_m], array_512_0);
+                }
+
+                if (tag_n_Mx != M) {
+                    array_512_0 = _mm512_maskz_loadu_ps(tail_mask, &C[idx_base_0+tag_n_Mx]);
+                    array_512_0 = _mm512_mul_ps(BETAVECTOR, array_512_0);
+                    _mm512_mask_storeu_ps(&C[idx_base_0+tag_n_Mx], tail_mask, array_512_0);
+                }
+                idx_base_0 += ldc;
+            }
+        }
+    } else {
+
+    }
+}
+
+// Scale matrix C while beta is not ZERO or ONE
+void sbgemm_zero_operation(OPENBLAS_CONST enum CBLAS_ORDER Order, OPENBLAS_CONST blasint M, OPENBLAS_CONST blasint N, float *C, OPENBLAS_CONST blasint ldc)
+{
+    BLASLONG tag_n_Nx = N & (~3);
+    BLASLONG tag_n_Mx = M & (~15);
+
+    BLASLONG LDC4x = ldc*4;
+    BLASLONG idx_base_0 = 0;
+    BLASLONG idx_base_1 = ldc;
+    BLASLONG idx_base_2 = ldc*2;
+    BLASLONG idx_base_3 = ldc*3;
+
+    unsigned short tail_mask_value = (((unsigned short)0xffff) >> (16-M+tag_n_Mx));
+    __mmask16 tail_mask = *((__mmask16*) &tail_mask_value);
+
+    __m512  ZEROVECTOR  = _mm512_setzero_ps();
+
+    if (Order == CblasColMajor) {
+        for (BLASLONG idx_n = 0; idx_n < tag_n_Nx; idx_n += 4) {
+            for (BLASLONG idx_m = 0; idx_m < tag_n_Mx; idx_m += 16) {
+                _mm512_storeu_ps(&C[idx_base_0+idx_m], ZEROVECTOR);
+                _mm512_storeu_ps(&C[idx_base_1+idx_m], ZEROVECTOR);
+                _mm512_storeu_ps(&C[idx_base_2+idx_m], ZEROVECTOR);
+                _mm512_storeu_ps(&C[idx_base_3+idx_m], ZEROVECTOR);
+            }
+
+            if (tag_n_Mx != M) {
+                _mm512_mask_storeu_ps(&C[idx_base_0+tag_n_Mx], tail_mask, ZEROVECTOR);
+                _mm512_mask_storeu_ps(&C[idx_base_1+tag_n_Mx], tail_mask, ZEROVECTOR);
+                _mm512_mask_storeu_ps(&C[idx_base_2+tag_n_Mx], tail_mask, ZEROVECTOR);
+                _mm512_mask_storeu_ps(&C[idx_base_3+tag_n_Mx], tail_mask, ZEROVECTOR);
+            }
+
+            idx_base_0 += LDC4x;
+            idx_base_1 += LDC4x;
+            idx_base_2 += LDC4x;
+            idx_base_3 += LDC4x;
+        }
+
+        if (tag_n_Nx != N) {
+            for (BLASLONG idx_n = tag_n_Nx; idx_n < N; idx_n++) {
+                for (BLASLONG idx_m = 0; idx_m < tag_n_Mx; idx_m += 16) {
+                    _mm512_storeu_ps(&C[idx_base_0+idx_m], ZEROVECTOR);
+                }
+
+                if (tag_n_Mx != M) {
+                    _mm512_mask_storeu_ps(&C[idx_base_0+tag_n_Mx], tail_mask, ZEROVECTOR);
+                }
+                idx_base_0 += ldc;
+            }
+        }
+    } else {
+
+    }
+}

kernel/x86_64/sbgemm_microk_cooperlake_template.c

@@ -0,0 +1,625 @@
+#include "sbgemm.h"
+#include "bf16_common_macros.h"
+#include <immintrin.h>
+
+#undef STORE16_COMPLETE_RESULT
+#undef STORE16_MASK_COMPLETE_RESULT
+#undef SBGEMM_BLOCK_KERNEL_32x8x32
+#undef SBGEMM_BLOCK_KERNEL_16x8x32
+#undef SBGEMM_BLOCK_KERNEL_32xNx32
+#undef SBGEMM_BLOCK_KERNEL_16xNx32
+#undef SBGEMM_BLOCKING_KERNEL_2
+
+#ifndef ONE_ALPHA      // ALPHA is not ONE
+    #define STORE16_COMPLETE_RESULT       STORE16_COMPLETE_RESULT_ALPHA_ONE
+    #define STORE16_MASK_COMPLETE_RESULT  STORE16_MASK_COMPLETE_RESULT_ALPHA_ONE
+    #define SBGEMM_BLOCK_KERNEL_32x8x32   sbgemm_block_kernel_32x8x32_alpha
+    #define SBGEMM_BLOCK_KERNEL_16x8x32   sbgemm_block_kernel_16x8x32_alpha
+    #define SBGEMM_BLOCK_KERNEL_32xNx32   sbgemm_block_kernel_32xNx32_alpha
+    #define SBGEMM_BLOCK_KERNEL_16xNx32   sbgemm_block_kernel_16xNx32_alpha
+    #define SBGEMM_BLOCKING_KERNEL_2      sbgemm_blocking_kernel_2_alpha
+#else                  // ALPHA is ONE
+    #define STORE16_COMPLETE_RESULT       STORE16_COMPLETE_RESULT_ONE_ONE
+    #define STORE16_MASK_COMPLETE_RESULT  STORE16_MASK_COMPLETE_RESULT_ONE_ONE
+    #define SBGEMM_BLOCK_KERNEL_32x8x32   sbgemm_block_kernel_32x8x32_one
+    #define SBGEMM_BLOCK_KERNEL_16x8x32   sbgemm_block_kernel_16x8x32_one
+    #define SBGEMM_BLOCK_KERNEL_32xNx32   sbgemm_block_kernel_32xNx32_one
+    #define SBGEMM_BLOCK_KERNEL_16xNx32   sbgemm_block_kernel_16xNx32_one
+    #define SBGEMM_BLOCKING_KERNEL_2      sbgemm_blocking_kernel_2_one
+#endif
+
+
+// SBGEMM Kernel for 16<M<=32, N=8, K can be any number, but the processing will take 32 as a base
+#ifndef ONE_ALPHA      // ALPHA is not ONE
+void sbgemm_block_kernel_32x8x32_alpha(BLASLONG m, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#else                  // ALPHA is ONE
+void sbgemm_block_kernel_32x8x32_one(BLASLONG m, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#endif
+{
+    int  SHUFFLE_MAGIC_NO = 0x39;
+    BLASLONG tag_k_32x = k & (~31);
+    BLASLONG idxA_base = 0;
+    BLASLONG idxB_base = 0;
+    BLASLONG width = 32;
+
+#ifndef ONE_ALPHA
+    __m512  ALPHAVECTOR = _mm512_set1_ps(alpha);
+#endif
+
+    __m512i arrayA_512_0, arrayA_512_1;
+    __m512i arrayB_512_0, arrayB_512_1, arrayB_512_2, arrayB_512_3, arrayB_512_4, arrayB_512_5, arrayB_512_6, arrayB_512_7;
+    __m512  result_512_0, result_512_1, result_512_2, result_512_3, result_512_4, result_512_5, result_512_6, result_512_7,
+            result_512_8, result_512_9, result_512_10, result_512_11, result_512_12, result_512_13, result_512_14, result_512_15;
+    __m512  result_512_tmp_0, result_512_tmp_1, result_512_tmp_2, result_512_tmp_3;
+
+    __m512i M512_EPI32_8      = _mm512_set1_epi32(8);
+    __m512i shuffle_idx_base0 = _mm512_set_epi32(23, 22, 21, 20, 7, 6, 5, 4, 19, 18, 17, 16, 3, 2, 1, 0);
+    __m512i shuffle_idx_base1 = _mm512_add_epi32(shuffle_idx_base0, M512_EPI32_8);
+
+    result_512_0  = _mm512_setzero_ps();
+    result_512_1  = _mm512_setzero_ps();
+    result_512_2  = _mm512_setzero_ps();
+    result_512_3  = _mm512_setzero_ps();
+    result_512_4  = _mm512_setzero_ps();
+    result_512_5  = _mm512_setzero_ps();
+    result_512_6  = _mm512_setzero_ps();
+    result_512_7  = _mm512_setzero_ps();
+    result_512_8  = _mm512_setzero_ps();
+    result_512_9  = _mm512_setzero_ps();
+    result_512_10 = _mm512_setzero_ps();
+    result_512_11 = _mm512_setzero_ps();
+    result_512_12 = _mm512_setzero_ps();
+    result_512_13 = _mm512_setzero_ps();
+    result_512_14 = _mm512_setzero_ps();
+    result_512_15 = _mm512_setzero_ps();
+
+    for (BLASLONG idx_k = 0; idx_k < k; idx_k += 32) {
+        // Load B with unroll 8
+        idxB_base = idx_k << 3;
+        arrayB_512_0 = _mm512_loadu_si512(&B[idxB_base + 32*0]);
+        arrayB_512_1 = _mm512_loadu_si512(&B[idxB_base + 32*1]);
+        arrayB_512_2 = _mm512_loadu_si512(&B[idxB_base + 32*2]);
+        arrayB_512_3 = _mm512_loadu_si512(&B[idxB_base + 32*3]);
+        arrayB_512_4 = _mm512_loadu_si512(&B[idxB_base + 32*4]);
+        arrayB_512_5 = _mm512_loadu_si512(&B[idxB_base + 32*5]);
+        arrayB_512_6 = _mm512_loadu_si512(&B[idxB_base + 32*6]);
+        arrayB_512_7 = _mm512_loadu_si512(&B[idxB_base + 32*7]);
+
+        if (idx_k == tag_k_32x) {width = k - tag_k_32x;}
+
+        for (BLASLONG idx = 0; idx < width;) {
+            // Each two rows are a group for 32-pair bf16 elements
+            idxA_base = idx << 5;
+            arrayA_512_0 = _mm512_loadu_si512(&A[idxA_base]);
+            arrayA_512_1 = _mm512_loadu_si512(&A[idxA_base + 32]);
+
+            result_512_0  = _mm512_dpbf16_ps(result_512_0,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_0)));
+            result_512_1  = _mm512_dpbf16_ps(result_512_1,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_1)));
+            result_512_2  = _mm512_dpbf16_ps(result_512_2,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_2)));
+            result_512_3  = _mm512_dpbf16_ps(result_512_3,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_3)));
+            result_512_4  = _mm512_dpbf16_ps(result_512_4,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_4)));
+            result_512_5  = _mm512_dpbf16_ps(result_512_5,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_5)));
+            result_512_6  = _mm512_dpbf16_ps(result_512_6,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_6)));
+            result_512_7  = _mm512_dpbf16_ps(result_512_7,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_7)));
+            result_512_8  = _mm512_dpbf16_ps(result_512_8,  (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_0)));
+            result_512_9  = _mm512_dpbf16_ps(result_512_9,  (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_1)));
+            result_512_10 = _mm512_dpbf16_ps(result_512_10, (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_2)));
+            result_512_11 = _mm512_dpbf16_ps(result_512_11, (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_3)));
+            result_512_12 = _mm512_dpbf16_ps(result_512_12, (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_4)));
+            result_512_13 = _mm512_dpbf16_ps(result_512_13, (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_5)));
+            result_512_14 = _mm512_dpbf16_ps(result_512_14, (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_6)));
+            result_512_15 = _mm512_dpbf16_ps(result_512_15, (__m512bh) arrayA_512_1, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_7)));
+
+            arrayB_512_0 = _mm512_shuffle_epi32(arrayB_512_0, SHUFFLE_MAGIC_NO);
+            arrayB_512_1 = _mm512_shuffle_epi32(arrayB_512_1, SHUFFLE_MAGIC_NO);
+            arrayB_512_2 = _mm512_shuffle_epi32(arrayB_512_2, SHUFFLE_MAGIC_NO);
+            arrayB_512_3 = _mm512_shuffle_epi32(arrayB_512_3, SHUFFLE_MAGIC_NO);
+            arrayB_512_4 = _mm512_shuffle_epi32(arrayB_512_4, SHUFFLE_MAGIC_NO);
+            arrayB_512_5 = _mm512_shuffle_epi32(arrayB_512_5, SHUFFLE_MAGIC_NO);
+            arrayB_512_6 = _mm512_shuffle_epi32(arrayB_512_6, SHUFFLE_MAGIC_NO);
+            arrayB_512_7 = _mm512_shuffle_epi32(arrayB_512_7, SHUFFLE_MAGIC_NO);
+
+            idx += 2;
+            // Every 4 loops we need to switch to next 128 bits of arrayB registers
+            if ((idx & (~7)) == idx) {
+                arrayB_512_0 = _mm512_shuffle_i32x4(arrayB_512_0, arrayB_512_0, SHUFFLE_MAGIC_NO);
+                arrayB_512_1 = _mm512_shuffle_i32x4(arrayB_512_1, arrayB_512_1, SHUFFLE_MAGIC_NO);
+                arrayB_512_2 = _mm512_shuffle_i32x4(arrayB_512_2, arrayB_512_2, SHUFFLE_MAGIC_NO);
+                arrayB_512_3 = _mm512_shuffle_i32x4(arrayB_512_3, arrayB_512_3, SHUFFLE_MAGIC_NO);
+                arrayB_512_4 = _mm512_shuffle_i32x4(arrayB_512_4, arrayB_512_4, SHUFFLE_MAGIC_NO);
+                arrayB_512_5 = _mm512_shuffle_i32x4(arrayB_512_5, arrayB_512_5, SHUFFLE_MAGIC_NO);
+                arrayB_512_6 = _mm512_shuffle_i32x4(arrayB_512_6, arrayB_512_6, SHUFFLE_MAGIC_NO);
+                arrayB_512_7 = _mm512_shuffle_i32x4(arrayB_512_7, arrayB_512_7, SHUFFLE_MAGIC_NO);
+            }
+        }
+    }
+
+    if (m != 32) {
+        unsigned short tail_mask_value = (((unsigned short)0xffff) >> (32-m));
+        __mmask16 tail_mask = *((__mmask16*) &tail_mask_value);
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_0, shuffle_idx_base0, result_512_8);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_0, shuffle_idx_base1, result_512_8);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_1, shuffle_idx_base0, result_512_9);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_1, shuffle_idx_base1, result_512_9);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*0]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*0+16]), tail_mask)
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*1]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*1+16]), tail_mask)
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_2, shuffle_idx_base0, result_512_10);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_2, shuffle_idx_base1, result_512_10);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_3, shuffle_idx_base0, result_512_11);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_3, shuffle_idx_base1, result_512_11);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*2]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*2+16]), tail_mask)
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*3]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*3+16]), tail_mask)
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_4, shuffle_idx_base0, result_512_12);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_4, shuffle_idx_base1, result_512_12);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_5, shuffle_idx_base0, result_512_13);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_5, shuffle_idx_base1, result_512_13);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*4]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*4+16]), tail_mask)
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*5]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*5+16]), tail_mask)
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_6, shuffle_idx_base0, result_512_14);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_6, shuffle_idx_base1, result_512_14);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_7, shuffle_idx_base0, result_512_15);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_7, shuffle_idx_base1, result_512_15);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*6]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*6+16]), tail_mask)
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*7]))
+        STORE16_MASK_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*7+16]), tail_mask)
+    } else {
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_0, shuffle_idx_base0, result_512_8);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_0, shuffle_idx_base1, result_512_8);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_1, shuffle_idx_base0, result_512_9);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_1, shuffle_idx_base1, result_512_9);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*0]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*0+16]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*1]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*1+16]))
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_2, shuffle_idx_base0, result_512_10);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_2, shuffle_idx_base1, result_512_10);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_3, shuffle_idx_base0, result_512_11);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_3, shuffle_idx_base1, result_512_11);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*2]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*2+16]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*3]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*3+16]))
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_4, shuffle_idx_base0, result_512_12);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_4, shuffle_idx_base1, result_512_12);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_5, shuffle_idx_base0, result_512_13);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_5, shuffle_idx_base1, result_512_13);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*4]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*4+16]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*5]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*5+16]))
+        result_512_tmp_0 = _mm512_permutex2var_ps(result_512_6, shuffle_idx_base0, result_512_14);
+        result_512_tmp_1 = _mm512_permutex2var_ps(result_512_6, shuffle_idx_base1, result_512_14);
+        result_512_tmp_2 = _mm512_permutex2var_ps(result_512_7, shuffle_idx_base0, result_512_15);
+        result_512_tmp_3 = _mm512_permutex2var_ps(result_512_7, shuffle_idx_base1, result_512_15);
+        STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*6]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*6+16]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_2, (&C[ldc*7]))
+        STORE16_COMPLETE_RESULT(result_512_tmp_3, (&C[ldc*7+16]))
+    }
+}
+
+// SBGEMM Kernel for M<=16, N=8, K can be any number, but the processing will take 32 as a base
+#ifndef ONE_ALPHA      // ALPHA is not ONE
+void sbgemm_block_kernel_16x8x32_alpha(BLASLONG m, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#else                  // ALPHA is ONE
+void sbgemm_block_kernel_16x8x32_one(BLASLONG m, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#endif
+{
+    int  SHUFFLE_MAGIC_NO = 0x39;
+    BLASLONG tag_k_32x = k & (~31);
+    BLASLONG idxB_base = 0;
+    BLASLONG width = 32;
+
+#ifndef ONE_ALPHA
+    __m512  ALPHAVECTOR = _mm512_set1_ps(alpha);
+#endif
+
+    __m512i arrayA_512_0;
+    __m512i arrayB_512_0, arrayB_512_1, arrayB_512_2, arrayB_512_3, arrayB_512_4, arrayB_512_5, arrayB_512_6, arrayB_512_7;
+    __m512  result_512_0, result_512_1, result_512_2, result_512_3, result_512_4, result_512_5, result_512_6, result_512_7;
+
+    result_512_0  = _mm512_setzero_ps();
+    result_512_1  = _mm512_setzero_ps();
+    result_512_2  = _mm512_setzero_ps();
+    result_512_3  = _mm512_setzero_ps();
+    result_512_4  = _mm512_setzero_ps();
+    result_512_5  = _mm512_setzero_ps();
+    result_512_6  = _mm512_setzero_ps();
+    result_512_7  = _mm512_setzero_ps();
+
+    for (BLASLONG idx_k = 0; idx_k < k; idx_k += 32) {
+        // Load B with unroll 8
+        idxB_base = idx_k << 3;
+        arrayB_512_0 = _mm512_loadu_si512(&B[idxB_base + 32*0]);
+        arrayB_512_1 = _mm512_loadu_si512(&B[idxB_base + 32*1]);
+        arrayB_512_2 = _mm512_loadu_si512(&B[idxB_base + 32*2]);
+        arrayB_512_3 = _mm512_loadu_si512(&B[idxB_base + 32*3]);
+        arrayB_512_4 = _mm512_loadu_si512(&B[idxB_base + 32*4]);
+        arrayB_512_5 = _mm512_loadu_si512(&B[idxB_base + 32*5]);
+        arrayB_512_6 = _mm512_loadu_si512(&B[idxB_base + 32*6]);
+        arrayB_512_7 = _mm512_loadu_si512(&B[idxB_base + 32*7]);
+
+        if (idx_k == tag_k_32x) {width = k - tag_k_32x;}
+
+        for (BLASLONG idx = 0; idx < width;) {
+            // Each two rows are a group for 32-pair bf16 elements
+            // Load two rows into a 512 register
+            arrayA_512_0 = _mm512_loadu_si512(&A[idx<<4]);
+
+            result_512_0  = _mm512_dpbf16_ps(result_512_0,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_0)));
+            result_512_1  = _mm512_dpbf16_ps(result_512_1,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_1)));
+            result_512_2  = _mm512_dpbf16_ps(result_512_2,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_2)));
+            result_512_3  = _mm512_dpbf16_ps(result_512_3,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_3)));
+            result_512_4  = _mm512_dpbf16_ps(result_512_4,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_4)));
+            result_512_5  = _mm512_dpbf16_ps(result_512_5,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_5)));
+            result_512_6  = _mm512_dpbf16_ps(result_512_6,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_6)));
+            result_512_7  = _mm512_dpbf16_ps(result_512_7,  (__m512bh) arrayA_512_0, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512_7)));
+
+            arrayB_512_0 = _mm512_shuffle_epi32(arrayB_512_0, SHUFFLE_MAGIC_NO);
+            arrayB_512_1 = _mm512_shuffle_epi32(arrayB_512_1, SHUFFLE_MAGIC_NO);
+            arrayB_512_2 = _mm512_shuffle_epi32(arrayB_512_2, SHUFFLE_MAGIC_NO);
+            arrayB_512_3 = _mm512_shuffle_epi32(arrayB_512_3, SHUFFLE_MAGIC_NO);
+            arrayB_512_4 = _mm512_shuffle_epi32(arrayB_512_4, SHUFFLE_MAGIC_NO);
+            arrayB_512_5 = _mm512_shuffle_epi32(arrayB_512_5, SHUFFLE_MAGIC_NO);
+            arrayB_512_6 = _mm512_shuffle_epi32(arrayB_512_6, SHUFFLE_MAGIC_NO);
+            arrayB_512_7 = _mm512_shuffle_epi32(arrayB_512_7, SHUFFLE_MAGIC_NO);
+
+            idx += 2;
+            // Every 4 loops we need to switch to next 128 bits of arrayB registers
+            if ((idx & (~7)) == idx) {
+                arrayB_512_0 = _mm512_shuffle_i32x4(arrayB_512_0, arrayB_512_0, SHUFFLE_MAGIC_NO);
+                arrayB_512_1 = _mm512_shuffle_i32x4(arrayB_512_1, arrayB_512_1, SHUFFLE_MAGIC_NO);
+                arrayB_512_2 = _mm512_shuffle_i32x4(arrayB_512_2, arrayB_512_2, SHUFFLE_MAGIC_NO);
+                arrayB_512_3 = _mm512_shuffle_i32x4(arrayB_512_3, arrayB_512_3, SHUFFLE_MAGIC_NO);
+                arrayB_512_4 = _mm512_shuffle_i32x4(arrayB_512_4, arrayB_512_4, SHUFFLE_MAGIC_NO);
+                arrayB_512_5 = _mm512_shuffle_i32x4(arrayB_512_5, arrayB_512_5, SHUFFLE_MAGIC_NO);
+                arrayB_512_6 = _mm512_shuffle_i32x4(arrayB_512_6, arrayB_512_6, SHUFFLE_MAGIC_NO);
+                arrayB_512_7 = _mm512_shuffle_i32x4(arrayB_512_7, arrayB_512_7, SHUFFLE_MAGIC_NO);
+            }
+        }
+    }
+
+    if (m != 16) {
+        unsigned short tail_mask_value = (((unsigned short)0xffff) >> (16-m));
+        __mmask16 tail_mask = *((__mmask16*) &tail_mask_value);
+
+        result_512_0 = _mm512_shuffle_f32x4(result_512_0, result_512_0, 0xd8);
+        result_512_1 = _mm512_shuffle_f32x4(result_512_1, result_512_1, 0xd8);
+        result_512_2 = _mm512_shuffle_f32x4(result_512_2, result_512_2, 0xd8);
+        result_512_3 = _mm512_shuffle_f32x4(result_512_3, result_512_3, 0xd8);
+        STORE16_MASK_COMPLETE_RESULT(result_512_0, (&C[ldc*0]), tail_mask)
+        STORE16_MASK_COMPLETE_RESULT(result_512_1, (&C[ldc*1]), tail_mask)
+        STORE16_MASK_COMPLETE_RESULT(result_512_2, (&C[ldc*2]), tail_mask)
+        STORE16_MASK_COMPLETE_RESULT(result_512_3, (&C[ldc*3]), tail_mask)
+        result_512_4 = _mm512_shuffle_f32x4(result_512_4, result_512_4, 0xd8);
+        result_512_5 = _mm512_shuffle_f32x4(result_512_5, result_512_5, 0xd8);
+        result_512_6 = _mm512_shuffle_f32x4(result_512_6, result_512_6, 0xd8);
+        result_512_7 = _mm512_shuffle_f32x4(result_512_7, result_512_7, 0xd8);
+        STORE16_MASK_COMPLETE_RESULT(result_512_4, (&C[ldc*4]), tail_mask)
+        STORE16_MASK_COMPLETE_RESULT(result_512_5, (&C[ldc*5]), tail_mask)
+        STORE16_MASK_COMPLETE_RESULT(result_512_6, (&C[ldc*6]), tail_mask)
+        STORE16_MASK_COMPLETE_RESULT(result_512_7, (&C[ldc*7]), tail_mask)
+    } else {
+        result_512_0 = _mm512_shuffle_f32x4(result_512_0, result_512_0, 0xd8);
+        result_512_1 = _mm512_shuffle_f32x4(result_512_1, result_512_1, 0xd8);
+        result_512_2 = _mm512_shuffle_f32x4(result_512_2, result_512_2, 0xd8);
+        result_512_3 = _mm512_shuffle_f32x4(result_512_3, result_512_3, 0xd8);
+        STORE16_COMPLETE_RESULT(result_512_0, (&C[ldc*0]))
+        STORE16_COMPLETE_RESULT(result_512_1, (&C[ldc*1]))
+        STORE16_COMPLETE_RESULT(result_512_2, (&C[ldc*2]))
+        STORE16_COMPLETE_RESULT(result_512_3, (&C[ldc*3]))
+        result_512_4 = _mm512_shuffle_f32x4(result_512_4, result_512_4, 0xd8);
+        result_512_5 = _mm512_shuffle_f32x4(result_512_5, result_512_5, 0xd8);
+        result_512_6 = _mm512_shuffle_f32x4(result_512_6, result_512_6, 0xd8);
+        result_512_7 = _mm512_shuffle_f32x4(result_512_7, result_512_7, 0xd8);
+        STORE16_COMPLETE_RESULT(result_512_4, (&C[ldc*4]))
+        STORE16_COMPLETE_RESULT(result_512_5, (&C[ldc*5]))
+        STORE16_COMPLETE_RESULT(result_512_6, (&C[ldc*6]))
+        STORE16_COMPLETE_RESULT(result_512_7, (&C[ldc*7]))
+    }
+}
+
+// SBGEMM Kernel for 16<M<=32, N<8, K can be any number, but the processing will take 32 as a base
+#ifndef ONE_ALPHA      // ALPHA is not ONE
+void sbgemm_block_kernel_32xNx32_alpha(BLASLONG m, BLASLONG n, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#else                  // ALPHA is ONE
+void sbgemm_block_kernel_32xNx32_one(BLASLONG m, BLASLONG n, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#endif
+{
+    int  SHUFFLE_MAGIC_NO = 0x39;
+    BLASLONG tag_k_32x = k & (~31);
+    BLASLONG idxA_base = 0;
+    BLASLONG idxB_base = 0;
+    BLASLONG width = 32;
+
+#ifndef ONE_ALPHA
+    __m512  ALPHAVECTOR = _mm512_set1_ps(alpha);
+#endif
+
+    __m512i arrayA_512[2];
+    __m512i arrayB_512[8];
+    __m512  result_512[16];
+    __m512  result_512_tmp_0, result_512_tmp_1;
+
+    __m512i M512_EPI32_8      = _mm512_set1_epi32(8);
+    __m512i shuffle_idx_base0 = _mm512_set_epi32(23, 22, 21, 20, 7, 6, 5, 4, 19, 18, 17, 16, 3, 2, 1, 0);
+    __m512i shuffle_idx_base1 = _mm512_add_epi32(shuffle_idx_base0, M512_EPI32_8);
+
+    for (int i = 0; i < 15; i += 2) {
+        result_512[i]    = _mm512_setzero_ps();
+        result_512[i+1]  = _mm512_setzero_ps();
+    }
+
+    for (BLASLONG idx_k = 0; idx_k < k; idx_k += 32) {
+        // Load B with unroll n
+        for (int i = 0; i < n; i ++) {
+            arrayB_512[i] = _mm512_loadu_si512(&B[idxB_base]);
+            idxB_base += 32;
+        }
+
+        if (idx_k == tag_k_32x) {width = k - tag_k_32x;}
+
+        for (BLASLONG idx = 0; idx < width;) {
+            // Each two rows are a group for 32-pair bf16 elements
+            idxA_base = idx << 5;
+            arrayA_512[0] = _mm512_loadu_si512(&A[idxA_base]);
+            arrayA_512[1] = _mm512_loadu_si512(&A[idxA_base + 32]);
+
+            for (int i = 0; i < n; i++) {
+                result_512[i]   = _mm512_dpbf16_ps(result_512[i]  , (__m512bh) arrayA_512[0], (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512[i])));
+                result_512[i+8] = _mm512_dpbf16_ps(result_512[i+8], (__m512bh) arrayA_512[1], (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512[i])));
+                arrayB_512[i]   = _mm512_shuffle_epi32(arrayB_512[i], SHUFFLE_MAGIC_NO);
+            }
+
+            idx += 2;
+            // Every 4 loops we need to switch to next 128 bits of arrayB registers
+            if ((idx & (~7)) == idx) {
+                for (int i = 0; i < n; i++) {
+                    arrayB_512[i] = _mm512_shuffle_i32x4(arrayB_512[i], arrayB_512[i], SHUFFLE_MAGIC_NO);
+                }
+            }
+        }
+    }
+
+    if (m != 32) {
+        unsigned short tail_mask_value = (((unsigned short)0xffff) >> (32-m));
+        __mmask16 tail_mask = *((__mmask16*) &tail_mask_value);
+        for (int i = 0; i < n; i++) {
+            result_512_tmp_0 = _mm512_permutex2var_ps(result_512[i], shuffle_idx_base0, result_512[i+8]);
+            result_512_tmp_1 = _mm512_permutex2var_ps(result_512[i], shuffle_idx_base1, result_512[i+8]);
+            STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*i]))
+            STORE16_MASK_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*i+16]), tail_mask)
+        }
+    } else {
+        for (int i = 0; i < n; i++) {
+            result_512_tmp_0 = _mm512_permutex2var_ps(result_512[i], shuffle_idx_base0, result_512[i+8]);
+            result_512_tmp_1 = _mm512_permutex2var_ps(result_512[i], shuffle_idx_base1, result_512[i+8]);
+            STORE16_COMPLETE_RESULT(result_512_tmp_0, (&C[ldc*i]))
+            STORE16_COMPLETE_RESULT(result_512_tmp_1, (&C[ldc*i+16]))
+        }
+    }
+}
+
+// SBGEMM Kernel for 16<=M, N<8, K can be any number, but the processing will take 32 as a base
+#ifndef ONE_ALPHA      // ALPHA is not ONE
+void sbgemm_block_kernel_16xNx32_alpha(BLASLONG m, BLASLONG n, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#else                  // ALPHA is ONE
+void sbgemm_block_kernel_16xNx32_one(BLASLONG m, BLASLONG n, BLASLONG k, float alpha, bfloat16 *A, bfloat16 *B, float *C, int ldc)
+#endif
+{
+    int  SHUFFLE_MAGIC_NO = 0x39;
+    BLASLONG tag_k_32x = k & (~31);
+    BLASLONG idxB_base = 0;
+    BLASLONG width = 32;
+
+#ifndef ONE_ALPHA
+    __m512  ALPHAVECTOR = _mm512_set1_ps(alpha);
+#endif
+
+    __m512i arrayA_512;
+    __m512i arrayB_512[8];
+    __m512  result_512[8];
+
+    for (int i = 0; i < 8; i += 2) {
+        result_512[i]    = _mm512_setzero_ps();
+        result_512[i+1]  = _mm512_setzero_ps();
+    }
+
+    for (BLASLONG idx_k = 0; idx_k < k; idx_k += 32) {
+        // Load B with unroll n
+        for (int i = 0; i < n; i ++) {
+            arrayB_512[i] = _mm512_loadu_si512(&B[idxB_base]);
+            idxB_base += 32;
+        }
+
+        if (idx_k == tag_k_32x) {width = k - tag_k_32x;}
+
+        for (BLASLONG idx = 0; idx < width;) {
+            // Each two rows are a group for 32-pair bf16 elements
+            // Load two rows into a 512 register
+            arrayA_512 = _mm512_loadu_si512(&A[idx<<4]);
+
+            for (int i = 0; i < n; i ++) {
+                result_512[i]  = _mm512_dpbf16_ps(result_512[i],  (__m512bh) arrayA_512, (__m512bh) _mm512_broadcastd_epi32(_mm512_castsi512_si128(arrayB_512[i])));
+                arrayB_512[i] = _mm512_shuffle_epi32(arrayB_512[i], SHUFFLE_MAGIC_NO);
+            }
+
+            idx += 2;
+            // Every 4 loops we need to switch to next 128 bits of arrayB registers
+            if ((idx & (~7)) == idx) {
+                for (int i = 0; i < n; i++) {
+                    arrayB_512[i] = _mm512_shuffle_i32x4(arrayB_512[i], arrayB_512[i], SHUFFLE_MAGIC_NO);
+                }
+            }
+        }
+    }
+
+    if (m != 16) {
+        unsigned short tail_mask_value = (((unsigned short)0xffff) >> (16-m));
+        __mmask16 tail_mask = *((__mmask16*) &tail_mask_value);
+        for (int i = 0; i < n; i++) {
+            result_512[i] = _mm512_shuffle_f32x4(result_512[i], result_512[i], 0xd8);
+            STORE16_MASK_COMPLETE_RESULT(result_512[i], (&C[ldc*i]), tail_mask)
+        }
+    } else {
+        for (int i = 0; i < n; i++) {
+            result_512[i] = _mm512_shuffle_f32x4(result_512[i], result_512[i], 0xd8);
+            STORE16_COMPLETE_RESULT(result_512[i], (&C[ldc*i]))
+        }
+    }
+}
+#ifndef ONE_ALPHA      // ALPHA is not ONE
+void sbgemm_blocking_kernel_2_alpha(blasint M, blasint N, blasint K, float alpha, bfloat16 *A, blasint lda, bfloat16 *B, blasint ldb, float *C, blasint ldc, bfloat16 * block_A, bfloat16 * block_B)
+#else                  // ALPHA is ONE
+void sbgemm_blocking_kernel_2_one(blasint M, blasint N, blasint K, float alpha, bfloat16 *A, blasint lda, bfloat16 *B, blasint ldb, float *C, blasint ldc, bfloat16 * block_A, bfloat16 * block_B)
+#endif
+{
+    BLASLONG m_step, n_step, k_step, k_step_round32;
+    BLASLONG tag_m_Nx = M & (~(BF16_BLOCK_THRES_M-1));
+
+    BLASLONG n_from, n_to;
+    BLASLONG tag_n_Nx;
+
+    n_from = 0;
+    n_to = (BF16_BLOCK_THRES_N > N) ? N : BF16_BLOCK_THRES_N;
+    tag_n_Nx = n_to & (~(BF16_BLOCK_STEP_N-1));
+
+    k_step = (K > BF16_BLOCK_THRES_K) ? BF16_BLOCK_THRES_K : K;
+    k_step_round32 = k_step & (~31);
+    k_step_round32 = (k_step > k_step_round32) ? (k_step_round32 + 32) : k_step_round32;
+
+    if (M >= BF16_BLOCK_THRES_M) {
+        while (n_from < N) {
+            for (BLASLONG idx_k = 0; idx_k < K;) {
+                // Use Kx32 kernel when BF16_BLOCK_THRES_M==32, Kx16 kernel when BF16_BLOCK_THRES_M==16, ...
+                COL_MAJOR_INCOPY_KERNEL_Kx32(k_step, &A(idx_k, 0), lda, block_A);
+                // TODO: MT
+                for (BLASLONG idx_n = n_from; idx_n < tag_n_Nx; idx_n += BF16_BLOCK_STEP_N) {
+                    // Use 8x32 kernel when BF16_BLOCK_THRES_N==8, 4x32 kernel when BF16_BLOCK_THRES_N==4, ...
+                    COL_MAJOR_ONCOPY_KERNEL_8x32(k_step, &B(idx_n, idx_k), ldb, block_B + (idx_n-n_from)*k_step_round32);
+                    SBGEMM_BLOCK_KERNEL_32x8x32(32, k_step, alpha, block_A, block_B + (idx_n-n_from)*k_step_round32, &C(idx_n, 0), ldc);
+                }
+
+                if (tag_n_Nx != n_to) {
+                    n_step = n_to - tag_n_Nx;
+                    COL_MAJOR_ONCOPY_KERNEL_Nx32(n_step, k_step, &B(tag_n_Nx, idx_k), ldb, block_B + (tag_n_Nx-n_from)*k_step_round32);
+                    SBGEMM_BLOCK_KERNEL_32xNx32(32, n_step, k_step, alpha, block_A, block_B + (tag_n_Nx-n_from)*k_step_round32, &C(tag_n_Nx, 0), ldc);
+                }
+
+                for (BLASLONG idx_m = BF16_BLOCK_THRES_M; idx_m < tag_m_Nx; idx_m += BF16_BLOCK_THRES_M) {
+                    COL_MAJOR_INCOPY_KERNEL_Kx32(k_step, &A(idx_k, idx_m), lda, block_A);
+                    for (BLASLONG idx_n = n_from; idx_n < tag_n_Nx; idx_n += BF16_BLOCK_STEP_N) {
+                        SBGEMM_BLOCK_KERNEL_32x8x32(32, k_step, alpha, block_A, block_B + (idx_n-n_from)*k_step_round32, &C(idx_n, idx_m), ldc);
+                    }
+
+                    if (tag_n_Nx != n_to) {
+                        n_step = n_to - tag_n_Nx;
+                        SBGEMM_BLOCK_KERNEL_32xNx32(32, n_step, k_step, alpha, block_A, block_B + (tag_n_Nx-n_from)*k_step_round32, &C(tag_n_Nx, idx_m), ldc);
+                    }
+                }
+
+                if (tag_m_Nx != M) {
+                    m_step = M - tag_m_Nx;
+                    if (m_step > 16) {
+                        COL_MAJOR_INCOPY_KERNEL_Kx32m(k_step, m_step, &A(idx_k, tag_m_Nx), lda, block_A);
+                        for (BLASLONG idx_n = n_from; idx_n < tag_n_Nx; idx_n += BF16_BLOCK_STEP_N) {
+                            SBGEMM_BLOCK_KERNEL_32x8x32(m_step, k_step, alpha, block_A, block_B + (idx_n-n_from)*k_step_round32, &C(idx_n, tag_m_Nx), ldc);
+                        }
+
+                        if (tag_n_Nx != n_to) {
+                            n_step = n_to - tag_n_Nx;
+                            SBGEMM_BLOCK_KERNEL_32xNx32(m_step, n_step, k_step, alpha, block_A, block_B + (tag_n_Nx-n_from)*k_step_round32, &C(tag_n_Nx, tag_m_Nx), ldc);
+                        }
+                    } else if (m_step == 16) {
+                        COL_MAJOR_INCOPY_KERNEL_Kx16(k_step, m_step, &A(idx_k, tag_m_Nx), lda, block_A);
+                        for (BLASLONG idx_n = n_from; idx_n < tag_n_Nx; idx_n += BF16_BLOCK_STEP_N) {
+                            SBGEMM_BLOCK_KERNEL_16x8x32(m_step, k_step, alpha, block_A, block_B + (idx_n-n_from)*k_step_round32, &C(idx_n, tag_m_Nx), ldc);
+                        }
+
+                        if (tag_n_Nx != n_to) {
+                            n_step = n_to - tag_n_Nx;
+                            SBGEMM_BLOCK_KERNEL_16xNx32(m_step, n_step, k_step, alpha, block_A, block_B + (tag_n_Nx-n_from)*k_step_round32, &C(tag_n_Nx, tag_m_Nx), ldc);
+                        }
+                    } else {
+                        COL_MAJOR_INCOPY_KERNEL_Kx16m(k_step, m_step, &A(idx_k, tag_m_Nx), lda, block_A);
+                        for (BLASLONG idx_n = n_from; idx_n < tag_n_Nx; idx_n += BF16_BLOCK_STEP_N) {
+                            SBGEMM_BLOCK_KERNEL_16x8x32(m_step, k_step, alpha, block_A, block_B + (idx_n-n_from)*k_step_round32, &C(idx_n, tag_m_Nx), ldc);
+                        }
+
+                        if (tag_n_Nx != n_to) {
+                            n_step = n_to - tag_n_Nx;
+                            SBGEMM_BLOCK_KERNEL_16xNx32(m_step, n_step, k_step, alpha, block_A, block_B + (tag_n_Nx-n_from)*k_step_round32, &C(tag_n_Nx, tag_m_Nx), ldc);
+                        }
+                    }
+                }
+
+                idx_k += k_step;
+                k_step = K - idx_k;
+                k_step = (k_step > BF16_BLOCK_THRES_K) ? BF16_BLOCK_THRES_K : k_step;
+                k_step_round32 = k_step & (~31);
+                k_step_round32 = (k_step > k_step_round32) ? (k_step_round32 + 32) : k_step_round32;
+            }
+
+            n_from = n_to;
+            n_to += BF16_BLOCK_THRES_N;
+            n_to = (n_to > N) ? N : n_to;
+            tag_n_Nx = n_to & (~(BF16_BLOCK_STEP_N-1));           
+        }
+    } else {
+        m_step = M - tag_m_Nx;
+        while (n_from < N) {
+            for (BLASLONG idx_k = 0; idx_k < K;) {
+                // Use Kx32 kernel when BF16_BLOCK_THRES_M==32, Kx16 kernel when BF16_BLOCK_THRES_M==16, ...
+                COL_MAJOR_INCOPY_KERNEL_Kx32m(k_step, m_step, &A(idx_k, 0), lda, block_A);
+                // TODO: MT
+                for (BLASLONG idx_n = n_from; idx_n < tag_n_Nx; idx_n += BF16_BLOCK_STEP_N) {
+                    // Use 8x32 kernel when BF16_BLOCK_THRES_N==8, 4x32 kernel when BF16_BLOCK_THRES_N==4, ...
+                    COL_MAJOR_ONCOPY_KERNEL_8x32(k_step, &B(idx_n, idx_k), ldb, block_B + (idx_n-n_from)*k_step_round32);
+                    SBGEMM_BLOCK_KERNEL_32x8x32(m_step, k_step, alpha, block_A, block_B + (idx_n-n_from)*k_step_round32, &C(idx_n, 0), ldc);
+                }
+
+                if (tag_n_Nx != n_to) {
+                    n_step = n_to - tag_n_Nx;
+                    COL_MAJOR_ONCOPY_KERNEL_Nx32(n_step, k_step, &B(tag_n_Nx, idx_k), ldb, block_B + (tag_n_Nx-n_from)*k_step_round32);
+                    SBGEMM_BLOCK_KERNEL_32xNx32(m_step, n_step, k_step, alpha, block_A, block_B + (tag_n_Nx-n_from)*k_step_round32, &C(tag_n_Nx, 0), ldc);
+                }
+
+                idx_k += k_step;
+                k_step = K - idx_k;
+                k_step = (k_step > BF16_BLOCK_THRES_K) ? BF16_BLOCK_THRES_K : k_step;
+                k_step_round32 = k_step & (~31);
+                k_step_round32 = (k_step > k_step_round32) ? (k_step_round32 + 32) : k_step_round32;
+            }
+            n_from = n_to;
+            n_to += BF16_BLOCK_THRES_N;
+            n_to = (n_to > N) ? N : n_to;
+            tag_n_Nx = n_to & (~(BF16_BLOCK_STEP_N-1));
+        }
+    }
+}
+
+#ifndef ONE_ALPHA      // ALPHA is not ONE
+void sbgemm_internal_kernel_alpha(OPENBLAS_CONST enum CBLAS_ORDER Order, OPENBLAS_CONST enum CBLAS_TRANSPOSE TransA, OPENBLAS_CONST enum CBLAS_TRANSPOSE TransB, OPENBLAS_CONST blasint M, OPENBLAS_CONST blasint N, OPENBLAS_CONST blasint K,
+		 OPENBLAS_CONST float alpha, OPENBLAS_CONST bfloat16 *A, OPENBLAS_CONST blasint lda, OPENBLAS_CONST bfloat16 *B, OPENBLAS_CONST blasint ldb, float *C, OPENBLAS_CONST blasint ldc)
+#else                  // ALPHA is ONE
+void sbgemm_internal_kernel_one(OPENBLAS_CONST enum CBLAS_ORDER Order, OPENBLAS_CONST enum CBLAS_TRANSPOSE TransA, OPENBLAS_CONST enum CBLAS_TRANSPOSE TransB, OPENBLAS_CONST blasint M, OPENBLAS_CONST blasint N, OPENBLAS_CONST blasint K,
+		 OPENBLAS_CONST float alpha, OPENBLAS_CONST bfloat16 *A, OPENBLAS_CONST blasint lda, OPENBLAS_CONST bfloat16 *B, OPENBLAS_CONST blasint ldb, float *C, OPENBLAS_CONST blasint ldc)
+#endif
+{
+    bfloat16 block_A[BF16_BLOCK_THRES_K * BF16_BLOCK_THRES_M];
+    bfloat16 block_B[BF16_BLOCK_THRES_N * BF16_BLOCK_THRES_K];
+
+    // TODO: assume no trans for both A and B, to complement these scenarios later
+    if (Order == CblasColMajor) {
+        SBGEMM_BLOCKING_KERNEL_2(M, N, K, alpha, A, lda, B, ldb, C, ldc, block_A, block_B);
+    } else {
+        
+    }
+}

kernel/x86_64/srot.c

@@ -2,7 +2,7 @@
 
 #if defined(SKYLAKEX)
 #include "srot_microk_skylakex-2.c"
-#elif defined(HASWELL)
+#elif defined(HASWELL) || defined(ZEN)
 #include "srot_microk_haswell-2.c"
 #endif
 

lapack-netlib/SRC/chgeqz.f

@@ -320,12 +320,13 @@
      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
       COMPLEX            ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
      $                   CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
-     $                   U12, X
+     $                   U12, X, ABI12, Y
 *     ..
 *     .. External Functions ..
+      COMPLEX            CLADIV
       LOGICAL            LSAME
       REAL               CLANHS, SLAMCH
-      EXTERNAL           LSAME, CLANHS, SLAMCH
+      EXTERNAL           CLADIV, LLSAME, CLANHS, SLAMCH
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           CLARTG, CLASET, CROT, CSCAL, XERBLA
@@ -729,22 +730,34 @@
             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
      $             ( BSCALE*T( ILAST, ILAST ) )
             ABI22 = AD22 - U12*AD21
+            ABI12 = AD12 - U12*AD11
 *
-            T1 = HALF*( AD11+ABI22 )
-            RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
-            TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) +
-     $             AIMAG( T1-ABI22 )*AIMAG( RTDISC )
-            IF( TEMP.LE.ZERO ) THEN
-               SHIFT = T1 + RTDISC
-            ELSE
-               SHIFT = T1 - RTDISC
+            SHIFT = ABI22
+            CTEMP = SQRT( ABI12 )*SQRT( AD21 )
+            TEMP = ABS1( CTEMP )
+            IF( CTEMP.NE.ZERO ) THEN
+               X = HALF*( AD11-SHIFT )
+               TEMP2 = ABS1( X )
+               TEMP = MAX( TEMP, ABS1( X ) )
+               Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
+               IF( TEMP2.GT.ZERO ) THEN
+                  IF( REAL( X / TEMP2 )*REAL( Y )+
+     $                AIMAG( X / TEMP2 )*AIMAG( Y ).LT.ZERO )Y = -Y
+               END IF
+               SHIFT = SHIFT - CTEMP*CLADIV( CTEMP, ( X+Y ) )
             END IF
          ELSE
 *
 *           Exceptional shift.  Chosen for no particularly good reason.
 *
-            ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
-     $                        (BSCALE*T(ILAST-1,ILAST-1))
+            IF( ( IITER / 20 )*20.EQ.IITER .AND. 
+     $         BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
+               ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
+     $            ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
+            ELSE
+               ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
+     $            ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
+            END IF
             SHIFT = ESHIFT
          END IF
 *

lapack-netlib/SRC/zhgeqz.f

@@ -320,12 +320,13 @@
      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
       COMPLEX*16         ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
      $                   CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
-     $                   U12, X
+     $                   U12, X, ABI12, Y
 *     ..
 *     .. External Functions ..
+      COMPLEX*16         ZLADIV
       LOGICAL            LSAME
       DOUBLE PRECISION   DLAMCH, ZLANHS
-      EXTERNAL           LSAME, DLAMCH, ZLANHS
+      EXTERNAL           ZLADIV, LSAME, DLAMCH, ZLANHS
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL
@@ -730,22 +731,34 @@
             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
      $             ( BSCALE*T( ILAST, ILAST ) )
             ABI22 = AD22 - U12*AD21
+            ABI12 = AD12 - U12*AD11
 *
-            T1 = HALF*( AD11+ABI22 )
-            RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
-            TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) +
-     $             DIMAG( T1-ABI22 )*DIMAG( RTDISC )
-            IF( TEMP.LE.ZERO ) THEN
-               SHIFT = T1 + RTDISC
-            ELSE
-               SHIFT = T1 - RTDISC
+            SHIFT = ABI22
+            CTEMP = SQRT( ABI12 )*SQRT( AD21 )
+            TEMP = ABS1( CTEMP )
+            IF( CTEMP.NE.ZERO ) THEN
+               X = HALF*( AD11-SHIFT )
+               TEMP2 = ABS1( X )
+               TEMP = MAX( TEMP, ABS1( X ) )
+               Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
+               IF( TEMP2.GT.ZERO ) THEN
+                  IF( DBLE( X / TEMP2 )*DBLE( Y )+
+     $                DIMAG( X / TEMP2 )*DIMAG( Y ).LT.ZERO )Y = -Y
+               END IF
+               SHIFT = SHIFT - CTEMP*ZLADIV( CTEMP, ( X+Y ) )
             END IF
          ELSE
 *
 *           Exceptional shift.  Chosen for no particularly good reason.
 *
-            ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
-     $                        (BSCALE*T(ILAST-1,ILAST-1))
+            IF( ( IITER / 20 )*20.EQ.IITER .AND. 
+     $         BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
+               ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
+     $            ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
+            ELSE
+               ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
+     $            ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
+            END IF
             SHIFT = ESHIFT
          END IF
 *

param.h

@@ -2399,6 +2399,9 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define GEMM_DEFAULT_OFFSET_B 65536
 #define GEMM_DEFAULT_ALIGN 0x0ffffUL
 
+#define SWITCH_RATIO            16
+#define GEMM_PREFERED_SIZE      16
+
 #define SGEMM_DEFAULT_UNROLL_M 16
 #define SGEMM_DEFAULT_UNROLL_N 8
 #define DGEMM_DEFAULT_UNROLL_M 16
@@ -2435,10 +2438,18 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define GEMM_DEFAULT_OFFSET_B 65536
 #define GEMM_DEFAULT_ALIGN 0x0ffffUL
 
+#define SWITCH_RATIO            16
+#define GEMM_PREFERED_SIZE      16
+
 #define SGEMM_DEFAULT_UNROLL_M 16
 #define SGEMM_DEFAULT_UNROLL_N 8
+#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
+#define DGEMM_DEFAULT_UNROLL_M 16
+#define DGEMM_DEFAULT_UNROLL_N 4
+#else
 #define DGEMM_DEFAULT_UNROLL_M 8
 #define DGEMM_DEFAULT_UNROLL_N 8
+#endif
 #define CGEMM_DEFAULT_UNROLL_M 8
 #define CGEMM_DEFAULT_UNROLL_N 4
 #define ZGEMM_DEFAULT_UNROLL_M 8