From: Georgi Gerganov Date: Sat, 21 Jun 2025 06:21:28 +0000 (+0300) Subject: tests : cleanup old tests (#1282) X-Git-Tag: upstream/0.0.2309~92 X-Git-Url: https://git.djapps.eu/?a=commitdiff_plain;h=b2a092a7000391461fb48764485e7105e78fed42;p=pkg%2Fggml%2Fsources%2Fggml tests : cleanup old tests (#1282) ggml-ci --- diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5263cfd5..2651a172 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -159,84 +159,6 @@ endif() # undefine NDEBUG so asserts don't get disabled in tests add_definitions(-UNDEBUG) -# -# test-vec0 - -set(TEST_TARGET test-vec0) -add_executable(${TEST_TARGET} ${TEST_TARGET}.c) -target_link_libraries(${TEST_TARGET} PRIVATE ggml) - -# -# test-vec1 (x86) -if (${CMAKE_SYSTEM_PROCESSOR} MATCHES "x86" AND "${CMAKE_C_FLAGS}" MATCHES "avx") - set(TEST_TARGET test-vec1) - add_executable(${TEST_TARGET} ${TEST_TARGET}.c) - target_link_libraries(${TEST_TARGET} PRIVATE ggml) -endif() - -# -# test-vec2 (arm) -if (${CMAKE_SYSTEM_PROCESSOR} MATCHES "arm") - set(TEST_TARGET test-vec2) - add_executable(${TEST_TARGET} ${TEST_TARGET}.c) - target_link_libraries(${TEST_TARGET} PRIVATE ggml) -endif() - -# -# test-mul-mat1 (arm) - -if (${CMAKE_SYSTEM_PROCESSOR} MATCHES "arm" AND NOT GGML_NO_ACCELERATE) - set(TEST_TARGET test-mul-mat1) - add_executable(${TEST_TARGET} ${TEST_TARGET}.c) - target_link_libraries(${TEST_TARGET} PRIVATE ggml ${GGML_EXTRA_LIBS}) - target_compile_options(${TEST_TARGET} PRIVATE ${GGML_EXTRA_FLAGS}) -endif() - -# -# test-blas0 (arm) - -if (${CMAKE_SYSTEM_PROCESSOR} MATCHES "arm" AND NOT GGML_NO_ACCELERATE) - set(TEST_TARGET test-blas0) - add_executable(${TEST_TARGET} ${TEST_TARGET}.c) - target_link_libraries(${TEST_TARGET} PRIVATE ggml ${GGML_EXTRA_LIBS}) - target_compile_options(${TEST_TARGET} PRIVATE ${GGML_EXTRA_FLAGS}) - add_test(NAME ${TEST_TARGET} COMMAND $ 128 128 128) - set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") -endif() - -# -# test-mul-mat2 - -set(TEST_TARGET test-mul-mat2) -add_executable(${TEST_TARGET} ${TEST_TARGET}.c) -target_link_libraries(${TEST_TARGET} PRIVATE ggml) -add_test(NAME ${TEST_TARGET} COMMAND $) -set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - -if (MATH_LIBRARY) - target_link_libraries(test-mul-mat2 PRIVATE ${MATH_LIBRARY}) -endif() - -# -# test0 - -set(TEST_TARGET test0) -add_executable(${TEST_TARGET} ${TEST_TARGET}.c) -target_link_libraries(${TEST_TARGET} PRIVATE ggml) -add_test(NAME ${TEST_TARGET} COMMAND $) -set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - -# -# test-svd0 (x86) - -if (${CMAKE_SYSTEM_PROCESSOR} MATCHES "x86" AND GGML_OPENBLAS) - set(TEST_TARGET test-svd0) - add_executable(${TEST_TARGET} ${TEST_TARGET}.c) - target_link_libraries(${TEST_TARGET} PRIVATE ggml ${GGML_EXTRA_LIBS}) - target_compile_options(${TEST_TARGET} PRIVATE ${GGML_EXTRA_FLAGS}) -endif() - - # # test-backend-ops @@ -274,19 +196,6 @@ if (NOT GGML_BACKEND_DL) add_test(NAME ${TEST_TARGET} COMMAND $) set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - # - # test-mul-mat0 - - set(TEST_TARGET test-mul-mat0) - add_executable(${TEST_TARGET} ${TEST_TARGET}.c) - target_link_libraries(${TEST_TARGET} PRIVATE ggml ${GGML_EXTRA_LIBS}) - if (MSVC) - target_link_options(${TEST_TARGET} PRIVATE "/STACK: 8388608") # 8MB - endif() - target_compile_options(${TEST_TARGET} PRIVATE ${GGML_EXTRA_FLAGS}) - add_test(NAME ${TEST_TARGET} COMMAND $) - set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - # # test-pool @@ -391,7 +300,6 @@ if (NOT GGML_BACKEND_DL) add_test(NAME ${TEST_TARGET} COMMAND $) set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - # # test-conv2d @@ -401,7 +309,6 @@ if (NOT GGML_BACKEND_DL) add_test(NAME ${TEST_TARGET} COMMAND $) set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - # # test-conv2d-dw @@ -411,16 +318,6 @@ if (NOT GGML_BACKEND_DL) add_test(NAME ${TEST_TARGET} COMMAND $) set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - - # - # test-mul-mat - - set(TEST_TARGET test-mul-mat) - add_executable(${TEST_TARGET} ${TEST_TARGET}.cpp) - target_link_libraries(${TEST_TARGET} PRIVATE ggml) - add_test(NAME ${TEST_TARGET} COMMAND $) - set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw") - # # test-cont diff --git a/tests/test-blas0.c b/tests/test-blas0.c deleted file mode 100644 index fcdfcb90..00000000 --- a/tests/test-blas0.c +++ /dev/null @@ -1,270 +0,0 @@ -#include "ggml.h" -#include "ggml-cpu.h" - -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -#include - -uint64_t get_time_us(void) { - struct timeval tv; - gettimeofday(&tv, NULL); - return tv.tv_sec * 1000000 + tv.tv_usec; -} - -// -// naive implementation -// - -void mul_mat_f32_0( - const float * restrict src0, // M x K - const float * restrict src1, // N x K (transposed) - float * dst, - int m, int n, int k) { - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - float sum = 0; - for (int l = 0; l < k; l++) { - sum += src0[i*k + l] * src1[j*k + l]; - } - dst[j*m + i] = sum; - } - } -} - -int main(int argc, const char ** argv) { - if (argc < 4) { - printf("Usage: %s M N K\n", argv[0]); - return 1; - } - - const int n_threads = 1; - - int M = atoi(argv[1]); - int N = atoi(argv[2]); - int K = atoi(argv[3]); - - srand(time(NULL)); - - if (M == 0) M = rand() % 1000 + 1; - if (N == 0) N = rand() % 1000 + 1; - if (K == 0) K = rand() % 1000 + 1; - - printf("M = %d, N = %d, K = %d\n", M, N, K); - - float * src0 = malloc(sizeof(float)*M*K); - float * src1 = malloc(sizeof(float)*N*K); - float * dst0 = malloc(sizeof(float)*M*N); // naive - float * dst1 = malloc(sizeof(float)*M*N); // blas - - struct ggml_init_params params = { - .mem_size = 2048ul*1024*1024, - .mem_buffer = NULL, - .no_alloc = false, - }; - - struct ggml_context * ctx0 = ggml_init(params); - - struct ggml_tensor * s0_f32 = ggml_new_tensor_2d(ctx0, GGML_TYPE_F32, K, M); - struct ggml_tensor * s1_f32 = ggml_new_tensor_2d(ctx0, GGML_TYPE_F32, K, N); - - struct ggml_tensor * s0_f16 = ggml_new_tensor_2d(ctx0, GGML_TYPE_F16, K, M); - struct ggml_tensor * s1_f16 = ggml_new_tensor_2d(ctx0, GGML_TYPE_F16, K, N); - - for (int j = 0; j < M; j++) { - for (int i = 0; i < K; i++) { - //src0[j*K + i] = j; - src0[j*K + i] = 1e-3*(rand() % 1000); - } - } - - for (int j = 0; j < N; j++) { - for (int i = 0; i < K; i++) { - //src1[j*K + i] = j + 1; - src1[j*K + i] = 1e-3*(rand() % 1000); - } - } - - // copy src0 to s0_f32 - { - float * p_f32 = s0_f32->data; - ggml_fp16_t * p_f16 = s0_f16->data; - for (int i = 0; i < M; i++) { - for (int j = 0; j < K; j++) { - p_f32[i*K + j] = src0[i*K + j]; - p_f16[i*K + j] = ggml_fp32_to_fp16(src0[i*K + j]); - } - } - } - - // copy src1 to s1_f32 - { - float * p_f32 = s1_f32->data; - ggml_fp16_t * p_f16 = s1_f16->data; - for (int i = 0; i < N; i++) { - for (int j = 0; j < K; j++) { - p_f32[i*K + j] = src1[i*K + j]; - p_f16[i*K + j] = ggml_fp32_to_fp16(src1[i*K + j]); - } - } - } - - const clock_t start = clock(); - const uint64_t start_us = get_time_us(); - - double iM = 1.0/M; - mul_mat_f32_0(src0, src1, dst0, M, N, K); - - // Use BLAS sgemm from Accelerate framework - cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, N, M, K, 1.0f, src1, K, src0, K, 0.0f, dst1, M); - - struct ggml_tensor * dst2 = NULL; - struct ggml_tensor * dst3 = NULL; - - { - dst2 = ggml_mul_mat(ctx0, s0_f32, s1_f32); - - struct ggml_cgraph * gf = ggml_new_graph(ctx0); - ggml_build_forward_expand(gf, dst2); - ggml_graph_compute_with_ctx(ctx0, gf, n_threads); - } - - { - dst3 = ggml_mul_mat(ctx0, s0_f16, s1_f32); - - struct ggml_cgraph * gf = ggml_new_graph(ctx0); - ggml_build_forward_expand(gf, dst3); - ggml_graph_compute_with_ctx(ctx0, gf, n_threads); - } - - bool ok_blas = true; - bool ok_ggml_f32 = true; - bool ok_ggml_f16 = true; - - // check BLAS - for (int i = 0; i < M*N; i++) { - if (fabs(dst0[i] - dst1[i])/fabs(dst0[i]) > 0.0001) { - printf("dst0[%d] = %f, dst1[%d] = %f\n", i, dst0[i], i, dst1[i]); - ok_blas = false; - } - } - - // check ggml (f32) - { - float * p = dst2->data; - for (int i = 0; i < M*N; i++) { - if (fabs(dst0[i] - p[i])/fabs(dst0[i]) > 0.0001) { - printf("dst0[%d] = %f, dst2[%d] = %f\n", i, dst0[i], i, p[i]); - ok_ggml_f32 = false; - } - } - } - - // check ggml (f16) - { - float * p = dst3->data; - for (int i = 0; i < M*N; i++) { - if (fabs(dst0[i] - p[i])/fabs(dst0[i]) > 0.01) { - printf("dst0[%d] = %f, dst3[%d] = %f\n", i, dst0[i], i, p[i]); - ok_ggml_f16 = false; - } - } - } - - { - const clock_t end = clock(); - const uint64_t end_us = get_time_us(); - printf("%s: elapsed ticks: %ld\n", __func__, end - start); - } - -#if 0 - // print src0 - printf("src0:\n"); - for (int i = 0; i < M; i++) { - for (int j = 0; j < K; j++) { - printf("%4.1f ", src0[i*K+j]); - } - printf("\n"); - } - - // print src1 - printf("src1:\n"); - for (int i = 0; i < N; i++) { - for (int j = 0; j < K; j++) { - printf("%4.1f ", src1[i*K+j]); - } - printf("\n"); - } - - printf("\n"); - printf("dst0 (naive):\n"); - for (int j = 0; j < N; j++) { - for (int i = 0; i < M; i++) { - printf("%4.1f ", dst0[j*M+i]); - } - printf("\n"); - } - - printf("\n"); - printf("dst1 (BLAS):\n"); - for (int j = 0; j < N; j++) { - for (int i = 0; i < M; i++) { - printf("%4.1f ", dst1[j*M+i]); - } - printf("\n"); - } - - printf("\n"); - printf("dst2 (ggml f32):\n"); - for (int j = 0; j < N; j++) { - for (int i = 0; i < M; i++) { - printf("%4.1f ", ((float *)dst2->data)[j*M+i]); - } - printf("\n"); - } - - printf("\n"); - printf("dst3 (ggml f16):\n"); - for (int j = 0; j < N; j++) { - for (int i = 0; i < M; i++) { - printf("%4.1f ", ((float *)dst3->data)[j*M+i]); - } - printf("\n"); - } - - printf("\n"); -#endif - - free(src0); - free(src1); - free(dst0); - free(dst1); - - ggml_free(ctx0); - - printf("ok_blas = %d\n", ok_blas); - if (!ok_blas) { - printf("ERROR: BLAS failed\n"); - } - - printf("ok_ggml_f32 = %d\n", ok_ggml_f32); - if (!ok_ggml_f32) { - printf("ERROR: ggml failed\n"); - } - - printf("ok_ggml_f16 = %d\n", ok_ggml_f16); - if (!ok_ggml_f16) { - printf("ERROR: ggml failed\n"); - } - - return (ok_blas && ok_ggml_f32 && ok_ggml_f16) ? 0 : 1; -} diff --git a/tests/test-mul-mat.cpp b/tests/test-mul-mat.cpp deleted file mode 100644 index 578d3e78..00000000 --- a/tests/test-mul-mat.cpp +++ /dev/null @@ -1,347 +0,0 @@ -#include "ggml.h" -#include "ggml-cpu.h" -#include "ggml-alloc.h" -#include "ggml-backend.h" - -#ifdef GGML_USE_CUDA -#include "ggml-cuda.h" -#endif - -#ifdef GGML_USE_METAL -#include "ggml-metal.h" -#endif - -#include -#include -#include -#include -#include -#include -#include -#include - -struct test_model { - struct ggml_tensor * a; - struct ggml_tensor * b; - ggml_backend_t backend = NULL; - ggml_backend_buffer_t buffer; - struct ggml_context * ctx; -}; - -void load_model(test_model & model, float* a, float* b, int M, int N, int K, bool use_gpu = false) { - size_t buffer_size = 0; - { - buffer_size += (M * N) * ggml_type_size(GGML_TYPE_F32); // tensor a - buffer_size += (N * K) * ggml_type_size(GGML_TYPE_F32); // tensor b - buffer_size += 1024; // overhead - } - - printf("%s: ggml tensor size = %d bytes\n", __func__, (int) sizeof(ggml_tensor)); - printf("%s: backend buffer size = %d bytes\n", __func__, (int) buffer_size); - - int num_tensors = 2; - struct ggml_init_params params { - /*.mem_size =*/ ggml_tensor_overhead() * num_tensors, - /*.mem_buffer =*/ NULL, - /*.no_alloc =*/ true, - }; - - // initialize the backend -#ifdef GGML_USE_CUDA - if (use_gpu) { - fprintf(stderr, "%s: using CUDA backend\n", __func__); - model.backend = ggml_backend_cuda_init(0); - if (!model.backend) { - fprintf(stderr, "%s: ggml_backend_cuda_init() failed\n", __func__); - } - } -#endif - -#ifdef GGML_USE_METAL - if (use_gpu) { - fprintf(stderr, "%s: using Metal backend\n", __func__); - model.backend = ggml_backend_metal_init(); - if (!model.backend) { - fprintf(stderr, "%s: ggml_backend_metal_init() failed\n", __func__); - } - } -#endif - - if(!model.backend) { - // fallback to CPU backend - model.backend = ggml_backend_cpu_init(); - } - - model.buffer = ggml_backend_alloc_buffer(model.backend, buffer_size); - - // create context - model.ctx = ggml_init(params); - - // create tensors - model.a = ggml_new_tensor_2d(model.ctx, GGML_TYPE_F32, K, M); - printf("Matrix A: [%i, %i]\n", K, M); - model.b = ggml_new_tensor_2d(model.ctx, GGML_TYPE_F32, K, N); - printf("Matrix B: [%i, %i]\n", K, N); - - // create a allocator - struct ggml_tallocr alloc = ggml_tallocr_new(model.buffer); - - // alloc memory - ggml_tallocr_alloc(&alloc, model.a); - - // load data to buffer - if(ggml_backend_is_cpu(model.backend) -#ifdef GGML_USE_METAL - || ggml_backend_is_metal(model.backend) -#endif - ) { - memcpy(model.a->data, a, ggml_nbytes(model.a)); - } else { - ggml_backend_tensor_set(model.a, a, 0, ggml_nbytes(model.a)); // cuda requires copy the data directly to device - } - - // alloc memory - ggml_tallocr_alloc(&alloc, model.b); - - if(ggml_backend_is_cpu(model.backend) -#ifdef GGML_USE_METAL - || ggml_backend_is_metal(model.backend) -#endif - ) { - memcpy(model.b->data, b, ggml_nbytes(model.b)); - } else { - ggml_backend_tensor_set(model.b, b, 0, ggml_nbytes(model.b)); // cuda requires copy the data directly to device - } -} - -struct ggml_cgraph * build_graph(const test_model& model) { - static size_t buf_size = ggml_tensor_overhead()*GGML_DEFAULT_GRAPH_SIZE + ggml_graph_overhead(); - static std::vector buf(buf_size); - - struct ggml_init_params params0 = { - /*.mem_size =*/ buf_size, - /*.mem_buffer =*/ buf.data(), - /*.no_alloc =*/ true, // the tensors will be allocated later by ggml_gallocr_alloc_graph() - }; - - // create a temporally context to build the graph - struct ggml_context * ctx0 = ggml_init(params0); - - struct ggml_cgraph * gf = ggml_new_graph(ctx0); - - // zT = x @ yT - struct ggml_tensor * result = ggml_mul_mat(ctx0, model.a, ggml_cont(ctx0, model.b)); - - // z = (zT)T - ggml_build_forward_expand(gf, ggml_cont(ctx0, ggml_transpose(ctx0, result))); - - // delete the temporally context used to build the graph - ggml_free(ctx0); - return gf; -} - -struct ggml_tensor* compute(const test_model & model, ggml_gallocr_t allocr) { - struct ggml_cgraph * gf = build_graph(model); - - // allocate tensors - ggml_gallocr_alloc_graph(allocr, gf); - int n_threads = 1; - - if (ggml_backend_is_cpu(model.backend)) { - ggml_backend_cpu_set_n_threads(model.backend, n_threads); - } - - - ggml_backend_graph_compute(model.backend, gf); - - //ggml_graph_print(gf); - - // in this case, the output tensor is the last one in the graph - return ggml_graph_node(gf, -1); -} - - -static void ggml_vec_dot_f16(const int n, float * s, float * x, float * y) { - float sumf = 0.0; - for (int i = 0; i < n; ++i) { - sumf += x[i] * y[i]; - } - *s = sumf; -} - -static void gemm_f16_out_f32(int m, int n, int k, - float * A, - float * B, - float * C, - const int ith, const int nth) { - // does not seem to make a difference - int m0, m1, n0, n1; - // patches per thread - if (m > n) { - n0 = 0; - n1 = n; - - // total patches in dst - const int np = m; - - // patches per thread - const int dp = (np + nth - 1)/nth; - - // patch range for this thread - m0 = dp*ith; - m1 = std::min(m0 + dp, np); - } else { - m0 = 0; - m1 = m; - - // total patches in dst - const int np = n; - - // patches per thread - const int dp = (np + nth - 1)/nth; - - // patch range for this thread - n0 = dp*ith; - n1 = std::min(n0 + dp, np); - } - - // block-tiling attempt - int64_t blck_n = 16; - int64_t blck_m = 16; - - for (int j = n0; j < n1; j+=blck_n) { - for (int i = m0; i < m1; i+=blck_m) { - // printf("i j k => %d %d %d\n", i, j, K); - for (int ii = i; ii < i + blck_m && ii < m1; ii++) { - for (int jj = j; jj < j + blck_n && jj < n1; jj++) { - ggml_vec_dot_f16(k, - C + ii*n + jj, - A + ii * k, - B + jj * k); - } - } - } - } -} - - -void perform_gemm_test(float* a, float* b, float* expected, int M, int N, int K) { - printf("\nPerforming gemm_f16_out_f32 test:\n"); - - std::vector gemm_out(M * N); - gemm_f16_out_f32(M, N, K, a, b, gemm_out.data(), 0, 1); - - for (int i = 0; i < M; i++) { - for (int j = 0; j < N; j++) { - printf("%.1ff,", gemm_out[i * N + j]); - } - printf("\n"); - } - - bool passed = true; - - for(int i = 0; i < M * N; i++) { - if(gemm_out[i] != expected[i]) { - passed = false; - break; - } - } - - printf("gemm_mult (%i): %s\n", (M * N), passed ? "\033[32mPASSED\033[0m" : "\033[31mFAILED\033[0m"); -} - -int main(void) -{ - ggml_time_init(); - const int M = 4, N = 16, K = 36; // a conv2d expected matrix multiplication - - // matrix A (4 X 36) - float matrixA[M * K] = { - 2.0f, 9.0f, 2.0f, 10.0f, 6.0f, 4.0f, 3.0f, 6.0f, 3.0f, 6.0f, 9.0f, 7.0f, 8.0f, 8.0f, 3.0f, 3.0f, 10.0f, 5.0f, 2.0f, 10.0f, 7.0f, 10.0f, 9.0f, 3.0f, 6.0f, 6.0f, 5.0f, 10.0f, 2.0f, 3.0f, 6.0f, 1.0f, 9.0f, 4.0f, 10.0f, 4.0f, - 10.0f, 7.0f, 8.0f, 10.0f, 10.0f, 8.0f, 7.0f, 10.0f, 4.0f, 6.0f, 8.0f, 7.0f, 7.0f, 6.0f, 9.0f, 3.0f, 6.0f, 5.0f, 5.0f, 2.0f, 7.0f, 2.0f, 7.0f, 4.0f, 4.0f, 6.0f, 6.0f, 4.0f, 3.0f, 9.0f, 3.0f, 6.0f, 4.0f, 7.0f, 2.0f, 9.0f, - 7.0f, 3.0f, 2.0f, 5.0f, 7.0f, 3.0f, 10.0f, 2.0f, 6.0f, 1.0f, 4.0f, 7.0f, 5.0f, 10.0f, 3.0f, 10.0f, 4.0f, 5.0f, 5.0f, 1.0f, 6.0f, 10.0f, 7.0f, 4.0f, 5.0f, 3.0f, 9.0f, 9.0f, 8.0f, 6.0f, 9.0f, 2.0f, 3.0f, 6.0f, 8.0f, 5.0f, - 5.0f, 5.0f, 5.0f, 5.0f, 3.0f, 10.0f, 4.0f, 1.0f, 8.0f, 8.0f, 9.0f, 8.0f, 4.0f, 1.0f, 4.0f, 9.0f, 3.0f, 6.0f, 3.0f, 1.0f, 4.0f, 8.0f, 3.0f, 10.0f, 8.0f, 6.0f, 4.0f, 5.0f, 4.0f, 3.0f, 2.0f, 2.0f, 4.0f, 3.0f, 6.0f, 4.0f, - }; - - // matrix B (16 X 36) - float matrixB[N * K] = { - 9.0f, 7.0f, 1.0f, 3.0f, 5.0f, 9.0f, 7.0f, 6.0f, 1.0f, 10.0f, 1.0f, 1.0f, 7.0f, 2.0f, 4.0f, 9.0f, 10.0f, 4.0f, 5.0f, 5.0f, 7.0f, 1.0f, 7.0f, 7.0f, 2.0f, 9.0f, 5.0f, 10.0f, 7.0f, 4.0f, 8.0f, 9.0f, 9.0f, 3.0f, 10.0f, 2.0f, - 4.0f, 6.0f, 10.0f, 9.0f, 5.0f, 1.0f, 8.0f, 7.0f, 4.0f, 7.0f, 2.0f, 6.0f, 5.0f, 3.0f, 1.0f, 10.0f, 8.0f, 4.0f, 8.0f, 3.0f, 7.0f, 1.0f, 2.0f, 7.0f, 6.0f, 8.0f, 6.0f, 5.0f, 2.0f, 3.0f, 1.0f, 1.0f, 2.0f, 5.0f, 7.0f, 1.0f, - 8.0f, 2.0f, 8.0f, 8.0f, 8.0f, 8.0f, 4.0f, 4.0f, 6.0f, 10.0f, 10.0f, 9.0f, 2.0f, 9.0f, 3.0f, 7.0f, 7.0f, 1.0f, 4.0f, 9.0f, 1.0f, 2.0f, 3.0f, 6.0f, 1.0f, 10.0f, 5.0f, 8.0f, 9.0f, 4.0f, 6.0f, 2.0f, 3.0f, 1.0f, 2.0f, 7.0f, - 5.0f, 1.0f, 7.0f, 2.0f, 9.0f, 10.0f, 9.0f, 5.0f, 2.0f, 5.0f, 4.0f, 10.0f, 9.0f, 9.0f, 1.0f, 9.0f, 8.0f, 8.0f, 9.0f, 4.0f, 9.0f, 4.0f, 8.0f, 2.0f, 1.0f, 8.0f, 4.0f, 5.0f, 10.0f, 7.0f, 6.0f, 2.0f, 1.0f, 10.0f, 10.0f, 7.0f, - 9.0f, 4.0f, 5.0f, 9.0f, 5.0f, 10.0f, 10.0f, 3.0f, 6.0f, 6.0f, 4.0f, 4.0f, 4.0f, 8.0f, 5.0f, 4.0f, 9.0f, 1.0f, 9.0f, 9.0f, 1.0f, 7.0f, 9.0f, 2.0f, 10.0f, 9.0f, 10.0f, 8.0f, 3.0f, 3.0f, 9.0f, 3.0f, 9.0f, 10.0f, 1.0f, 8.0f, - 9.0f, 2.0f, 6.0f, 9.0f, 7.0f, 2.0f, 3.0f, 5.0f, 3.0f, 6.0f, 9.0f, 7.0f, 3.0f, 7.0f, 6.0f, 4.0f, 10.0f, 3.0f, 5.0f, 7.0f, 2.0f, 9.0f, 3.0f, 2.0f, 2.0f, 10.0f, 8.0f, 7.0f, 3.0f, 10.0f, 6.0f, 3.0f, 1.0f, 1.0f, 4.0f, 10.0f, - 2.0f, 9.0f, 2.0f, 10.0f, 6.0f, 4.0f, 3.0f, 6.0f, 3.0f, 6.0f, 9.0f, 7.0f, 8.0f, 8.0f, 3.0f, 3.0f, 10.0f, 5.0f, 2.0f, 10.0f, 7.0f, 10.0f, 9.0f, 3.0f, 6.0f, 6.0f, 5.0f, 10.0f, 2.0f, 3.0f, 6.0f, 1.0f, 9.0f, 4.0f, 10.0f, 4.0f, - 10.0f, 7.0f, 8.0f, 10.0f, 10.0f, 8.0f, 7.0f, 10.0f, 4.0f, 6.0f, 8.0f, 7.0f, 7.0f, 6.0f, 9.0f, 3.0f, 6.0f, 5.0f, 5.0f, 2.0f, 7.0f, 2.0f, 7.0f, 4.0f, 4.0f, 6.0f, 6.0f, 4.0f, 3.0f, 9.0f, 3.0f, 6.0f, 4.0f, 7.0f, 2.0f, 9.0f, - 7.0f, 3.0f, 2.0f, 5.0f, 7.0f, 3.0f, 10.0f, 2.0f, 6.0f, 1.0f, 4.0f, 7.0f, 5.0f, 10.0f, 3.0f, 10.0f, 4.0f, 5.0f, 5.0f, 1.0f, 6.0f, 10.0f, 7.0f, 4.0f, 5.0f, 3.0f, 9.0f, 9.0f, 8.0f, 6.0f, 9.0f, 2.0f, 3.0f, 6.0f, 8.0f, 5.0f, - 5.0f, 5.0f, 5.0f, 5.0f, 3.0f, 10.0f, 4.0f, 1.0f, 8.0f, 8.0f, 9.0f, 8.0f, 4.0f, 1.0f, 4.0f, 9.0f, 3.0f, 6.0f, 3.0f, 1.0f, 4.0f, 8.0f, 3.0f, 10.0f, 8.0f, 6.0f, 4.0f, 5.0f, 4.0f, 3.0f, 2.0f, 2.0f, 4.0f, 3.0f, 6.0f, 4.0f, - 6.0f, 2.0f, 3.0f, 3.0f, 3.0f, 7.0f, 5.0f, 1.0f, 8.0f, 1.0f, 4.0f, 5.0f, 1.0f, 1.0f, 6.0f, 4.0f, 2.0f, 1.0f, 7.0f, 8.0f, 6.0f, 1.0f, 1.0f, 5.0f, 6.0f, 5.0f, 10.0f, 6.0f, 7.0f, 5.0f, 9.0f, 3.0f, 2.0f, 7.0f, 9.0f, 4.0f, - 2.0f, 5.0f, 9.0f, 5.0f, 10.0f, 3.0f, 1.0f, 8.0f, 1.0f, 7.0f, 1.0f, 8.0f, 1.0f, 6.0f, 7.0f, 8.0f, 4.0f, 9.0f, 5.0f, 10.0f, 3.0f, 7.0f, 6.0f, 8.0f, 8.0f, 5.0f, 6.0f, 8.0f, 10.0f, 9.0f, 4.0f, 1.0f, 3.0f, 3.0f, 4.0f, 7.0f, - 8.0f, 2.0f, 6.0f, 6.0f, 5.0f, 1.0f, 3.0f, 7.0f, 1.0f, 7.0f, 2.0f, 2.0f, 2.0f, 8.0f, 4.0f, 1.0f, 1.0f, 5.0f, 9.0f, 4.0f, 1.0f, 2.0f, 3.0f, 10.0f, 1.0f, 4.0f, 9.0f, 9.0f, 6.0f, 8.0f, 8.0f, 1.0f, 9.0f, 10.0f, 4.0f, 1.0f, - 8.0f, 5.0f, 8.0f, 9.0f, 4.0f, 8.0f, 2.0f, 1.0f, 1.0f, 9.0f, 4.0f, 5.0f, 6.0f, 1.0f, 2.0f, 5.0f, 6.0f, 7.0f, 3.0f, 1.0f, 4.0f, 6.0f, 7.0f, 7.0f, 7.0f, 8.0f, 7.0f, 8.0f, 8.0f, 2.0f, 10.0f, 2.0f, 7.0f, 3.0f, 8.0f, 3.0f, - 8.0f, 7.0f, 6.0f, 2.0f, 4.0f, 10.0f, 10.0f, 6.0f, 10.0f, 3.0f, 7.0f, 6.0f, 4.0f, 3.0f, 5.0f, 5.0f, 5.0f, 3.0f, 8.0f, 10.0f, 3.0f, 4.0f, 8.0f, 4.0f, 2.0f, 6.0f, 8.0f, 9.0f, 6.0f, 9.0f, 4.0f, 3.0f, 5.0f, 2.0f, 2.0f, 6.0f, - 10.0f, 6.0f, 2.0f, 1.0f, 7.0f, 5.0f, 6.0f, 4.0f, 1.0f, 9.0f, 10.0f, 2.0f, 4.0f, 5.0f, 8.0f, 5.0f, 7.0f, 4.0f, 7.0f, 6.0f, 3.0f, 9.0f, 2.0f, 1.0f, 4.0f, 2.0f, 6.0f, 6.0f, 3.0f, 3.0f, 2.0f, 8.0f, 5.0f, 9.0f, 3.0f, 4.0f, - }; - - // matrix C (4 x 16) - float expected_result[M * N] = { - 1224.0f, 1023.0f, 1158.0f,1259.0f,1359.0f,1194.0f,1535.0f,1247.0f,1185.0f,1029.0f,889.0f,1182.0f,955.0f,1179.0f,1147.0f,1048.0f, - 1216.0f, 1087.0f, 1239.0f,1361.0f,1392.0f,1260.0f,1247.0f,1563.0f,1167.0f,1052.0f,942.0f,1214.0f,1045.0f,1134.0f,1264.0f,1126.0f, - 1125.0f, 966.0f, 1079.0f,1333.0f,1287.0f,1101.0f,1185.0f,1167.0f,1368.0f,990.0f,967.0f,1121.0f,971.0f,1086.0f,1130.0f,980.0f, - 999.0f, 902.0f, 1020.0f,1056.0f,1076.0f,929.0f,1029.0f,1052.0f,990.0f,1108.0f,823.0f,989.0f,759.0f,1041.0f,1003.0f,870.0f - }; - - bool passed = true; - - perform_gemm_test(matrixA, matrixB, expected_result, M, N, K); - - test_model model; - load_model(model, matrixA, matrixB, M, N, K, true); - - ggml_gallocr_t allocr = NULL; - - { - allocr = ggml_gallocr_new(ggml_backend_get_default_buffer_type(model.backend)); - - //create the worst case graph for memory usage estimation - struct ggml_cgraph * gf = build_graph(model); - - // compute the required memory - ggml_gallocr_reserve(allocr, gf); - size_t mem_size = ggml_gallocr_get_buffer_size(allocr, 0); - fprintf(stderr, "%s: compute buffer size: %.2f MB\n", __func__, mem_size/1024.0f/1024.0f); - } - - struct ggml_tensor * result = compute(model, allocr); - - std::vector out_data(ggml_nelements(result)); - - ggml_backend_tensor_get(result, out_data.data(), 0, ggml_nbytes(result)); - - printf("\nPerforming ggml_mul_mat test:\n"); - - passed = true; - for(int i = 0; i < M * N; i++) { - if(out_data[i] != expected_result[i]) { - passed = false; - break; - } - } - - for (int i = 0; i < M; i++) { - for (int j = 0; j < N; j++) { - printf("%.1f ", out_data[i * N + j]); - } - printf("\n"); - } - - printf("ggml_mul_mat (%d): %s\n", (int) ggml_nelements(result), passed && (ggml_nelements(result) == M * N) ? "\033[32mPASSED\033[0m" : "\033[31mFAILED\033[0m"); - - // free memory - ggml_free(model.ctx); - - ggml_backend_buffer_free(model.buffer); - ggml_backend_free(model.backend); - ggml_gallocr_free(allocr); - return 0; -} diff --git a/tests/test-mul-mat0.c b/tests/test-mul-mat0.c deleted file mode 100644 index d254a161..00000000 --- a/tests/test-mul-mat0.c +++ /dev/null @@ -1,336 +0,0 @@ -#define _CRT_SECURE_NO_DEPRECATE // Disables ridiculous "unsafe" warnigns on Windows -#include "ggml.h" -#include "ggml-cpu.h" - -#include -#include -#include -#include -#include - -#if defined(_MSC_VER) -#pragma warning(disable: 4244 4267) // possible loss of data -#endif - -#define MAX_NARGS 2 - -float frand(void) { - return (float)rand()/(float)RAND_MAX; -} - -int irand(int n) { - return rand()%n; -} - -void get_random_dims(int64_t * dims, int ndims) { - dims[0] = dims[1] = dims[2] = dims[3] = 1; - - for (int i = 0; i < ndims; i++) { - dims[i] = 1 + irand(4); - } -} - -struct ggml_tensor * get_random_tensor( - struct ggml_context * ctx0, - int ndims, - int64_t ne[], - float fmin, - float fmax) { - struct ggml_tensor * result = ggml_new_tensor(ctx0, GGML_TYPE_F32, ndims, ne); - - switch (ndims) { - case 1: - for (int i0 = 0; i0 < ne[0]; i0++) { - ((float *)result->data)[i0] = frand()*(fmax - fmin) + fmin; - } - break; - case 2: - for (int i1 = 0; i1 < ne[1]; i1++) { - for (int i0 = 0; i0 < ne[0]; i0++) { - ((float *)result->data)[i1*ne[0] + i0] = frand()*(fmax - fmin) + fmin; - } - } - break; - case 3: - for (int i2 = 0; i2 < ne[2]; i2++) { - for (int i1 = 0; i1 < ne[1]; i1++) { - for (int i0 = 0; i0 < ne[0]; i0++) { - ((float *)result->data)[i2*ne[1]*ne[0] + i1*ne[0] + i0] = frand()*(fmax - fmin) + fmin; - } - } - } - break; - case 4: - for (int i3 = 0; i3 < ne[3]; i3++) { - for (int i2 = 0; i2 < ne[2]; i2++) { - for (int i1 = 0; i1 < ne[1]; i1++) { - for (int i0 = 0; i0 < ne[0]; i0++) { - ((float *)result->data)[i3*ne[2]*ne[1]*ne[0] + i2*ne[1]*ne[0] + i1*ne[0] + i0] = frand()*(fmax - fmin) + fmin; - } - } - } - } - break; - default: - assert(false); - }; - - return result; -} - -float get_element(const struct ggml_tensor * t, int idx) { - return ((float *)t->data)[idx]; -} - -void set_element(struct ggml_tensor * t, int idx, float value) { - ((float *)t->data)[idx] = value; -} - -bool check_gradient( - const char * op_name, - struct ggml_context * ctx0, - struct ggml_tensor * x[], - struct ggml_tensor * f, - int ndims, - int nargs, - float eps, - float max_error_abs, - float max_error_rel) { - const int n_threads = 1; - ggml_set_loss(f); - - struct ggml_cgraph * gf = ggml_new_graph_custom(ctx0, GGML_DEFAULT_GRAPH_SIZE, true); - ggml_build_forward_expand(gf, f); - struct ggml_cgraph * gb = ggml_graph_dup(ctx0, gf, false); - ggml_build_backward_expand(ctx0, gb, false); - - ggml_graph_compute_with_ctx(ctx0, gf, n_threads); - ggml_graph_reset(gb); - ggml_graph_compute_with_ctx(ctx0, gb, n_threads); - - ggml_graph_dump_dot(gf, NULL, "test-grad0-forward.dot"); - ggml_graph_dump_dot(gb, gf, "test-grad0-backward.dot"); - - for (int i = 0; i < nargs; ++i) { - const int64_t nelements = ggml_nelements(x[i]); - for (int64_t k = 0; k < nelements; ++k) { - // compute gradient using finite differences - const float x0 = get_element(x[i], k); - - set_element(x[i], k, x0 + eps); - ggml_graph_compute_with_ctx(ctx0, gf, n_threads); - - const float f0 = ggml_get_f32_1d(f, 0); - - set_element(x[i], k, x0 - eps); - ggml_graph_compute_with_ctx(ctx0, gf, n_threads); - - const float f1 = ggml_get_f32_1d(f, 0); - - const float g0 = (f0 - f1)/(2.0f*eps); - - set_element(x[i], k, x0); - - // compute gradient using backward graph - ggml_graph_reset(gb); - ggml_graph_compute_with_ctx(ctx0, gb, n_threads); - - const float g1 = get_element(ggml_graph_get_grad(gb, x[i]), k); - - const float error_abs = fabsf(g0 - g1); - const float error_rel = g0 != 0 ? fabsf(g0 - g1)/fabs(g0) : 0; - - if (error_abs > max_error_abs || error_rel > max_error_rel) { - printf("%s: ndims=%d, i=%d, k=%" PRId64 ", g0=%f, g1=%f, error_abs=%f, error_rel=%f\n", op_name, ndims, i, k, g0, g1, error_abs, error_rel); - assert(false); - } - } - } - - return true; -} - - -float mat_get(const struct ggml_tensor * t, int i0, int i1, int i2, int i3) { - const size_t nb0 = t->nb[0]; - const size_t nb1 = t->nb[1]; - const size_t nb2 = t->nb[2]; - const size_t nb3 = t->nb[3]; - - return - *((float*) ((char*)t->data + i0*nb0 + i1*nb1 + i2*nb2 + i3*nb3)); -} - -bool check_mat_mul( - const struct ggml_tensor * y, - const struct ggml_tensor * x0, - const struct ggml_tensor * x1) { - const int64_t n00 = x0->ne[0]; - const int64_t n10 = x0->ne[1]; - const int64_t n20 = x0->ne[2]; - const int64_t n30 = x0->ne[3]; - - const int64_t n01 = x1->ne[0]; - const int64_t n11 = x1->ne[1]; - const int64_t n21 = x1->ne[2]; - const int64_t n31 = x1->ne[3]; - - const int64_t n02 = y->ne[0]; - const int64_t n12 = y->ne[1]; - const int64_t n22 = y->ne[2]; - const int64_t n32 = y->ne[3]; - - printf("x0: [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "]\n", n00, n10, n20, n30); - for (int j = 0; j < n10; ++j) { - for (int i = 0; i < n00; ++i) { - printf("%6.3f ", mat_get(x0, i, j, 0, 0)); - } - printf("\n"); - } - printf("\n"); - - printf("x1: [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "]\n", n01, n11, n21, n31); - for (int j = 0; j < n11; ++j) { - for (int i = 0; i < n01; ++i) { - printf("%6.3f ", mat_get(x1, i, j, 0, 0)); - } - printf("\n"); - } - printf("\n"); - - printf("y: [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "]\n", n02, n12, n22, n32); - for (int j = 0; j < n12; ++j) { - for (int i = 0; i < n02; ++i) { - printf("%6.3f ", mat_get(y, i, j, 0, 0)); - } - printf("\n"); - } - - for (int i3 = 0; i3 < n32; ++i3) { - for (int i2 = 0; i2 < n22; ++i2) { - for (int i1 = 0; i1 < n12; ++i1) { - for (int i0 = 0; i0 < n02; ++i0) { - float sum = 0.0f; - for (int k = 0; k < n00; ++k) { - sum += mat_get(x0, k, i0, i2, i3) * mat_get(x1, k, i1, i2, i3); - } - if (fabsf(sum - mat_get(y, i0, i1, i2, i3)) > 1e-5) { - printf("error: i0=%d, i1=%d, i2=%d, i3=%d, sum=%f, y=%f\n", - i0, i1, i2, i3, sum, mat_get(y, i0, i1, i2, i3)); - assert(false); - return false; - } - } - } - } - } - - return true; -} - -int main(int argc, const char ** argv) { - struct ggml_init_params params = { - .mem_size = 128*1024*1024, - .mem_buffer = NULL, - .no_alloc = false, - }; - - int64_t ne[4]; - - // original loop: 500 - int niter = 500; - const char *env = getenv("GGML_NLOOP"); - if (env != NULL) { - niter = atoi(env); - } - if (argc > 1) { - niter = atoi(argv[1]); - } - - int n_threads = 1; - - for (int iter = 0; iter < niter; ++iter) { - printf("test-mul-mat0: iter:%d/%d\n", iter, niter); - struct ggml_context * ctx0 = ggml_init(params); - - get_random_dims(ne, 4); - - struct ggml_tensor * x[MAX_NARGS]; - - // mul_mat - { - const int nargs = 1; - - for (int ndims = 2; ndims <= 4; ++ndims) { - x[0] = get_random_tensor(ctx0, ndims, ne, -1.0f, 1.0f); - ne[1] = rand()%4 + 1; - x[1] = get_random_tensor(ctx0, ndims, ne, -1.0f, 1.0f); - - ggml_set_param(x[0]); - - struct ggml_tensor * m = ggml_mul_mat(ctx0, x[1], x[0]); - struct ggml_tensor * f = ggml_sum(ctx0, m); - - printf("testing: mul_mat, [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "] = [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "] * [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "]\n", - m->ne[0], m->ne[1], m->ne[2], m->ne[3], - x[1]->ne[0], x[1]->ne[1], x[1]->ne[2], x[1]->ne[3], - x[0]->ne[0], x[0]->ne[1], x[0]->ne[2], x[0]->ne[3]); - - assert(m->ne[0] == x[1]->ne[1]); - assert(m->ne[1] == x[0]->ne[1]); - assert(m->ne[2] == x[0]->ne[2]); - assert(m->ne[3] == x[0]->ne[3]); - - if (ndims <= 2) { - check_gradient("mul_mat", ctx0, x, f, ndims, nargs, 1e-3f, 1e-3f, INFINITY); - } else { - struct ggml_cgraph * gf = ggml_new_graph(ctx0); - ggml_build_forward_expand(gf, m); - ggml_graph_compute_with_ctx(ctx0, gf, n_threads); - } - - check_mat_mul(m, x[1], x[0]); - } - } - - // mul_mat (transposed) - { - const int nargs = 1; - - for (int ndims = 2; ndims <= 4; ++ndims) { - x[0] = get_random_tensor(ctx0, ndims, ne, -1.0f, 1.0f); - ne[1] = ne[0]; - ne[0] = rand()%4 + 1; - x[1] = ggml_cont(ctx0, ggml_transpose(ctx0, get_random_tensor(ctx0, ndims, ne, -1.0f, 1.0f))); - - ggml_set_param(x[0]); - - struct ggml_tensor * m = ggml_mul_mat(ctx0, x[1], x[0]); - struct ggml_tensor * f = ggml_sum(ctx0, m); - - printf("testing: mul_mat, [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "] = [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "] * [%" PRId64 ", %" PRId64 ", %" PRId64 ", %" PRId64 "]\n", - m->ne[0], m->ne[1], m->ne[2], m->ne[3], - x[1]->ne[0], x[1]->ne[1], x[1]->ne[2], x[1]->ne[3], - x[0]->ne[0], x[0]->ne[1], x[0]->ne[2], x[0]->ne[3]); - - assert(m->ne[0] == x[1]->ne[1]); - assert(m->ne[1] == x[0]->ne[1]); - assert(m->ne[2] == x[0]->ne[2]); - assert(m->ne[3] == x[0]->ne[3]); - - if (ndims <= 2) { - check_gradient("mul_mat", ctx0, x, f, ndims, nargs, 1e-3f, 1e-3f, INFINITY); - } else { - struct ggml_cgraph * gf = ggml_new_graph(ctx0); - ggml_build_forward_expand(gf, m); - ggml_graph_compute_with_ctx(ctx0, gf, n_threads); - } - - check_mat_mul(m, x[1], x[0]); - } - } - ggml_free(ctx0); - } - - return 0; -} diff --git a/tests/test-mul-mat1.c b/tests/test-mul-mat1.c deleted file mode 100644 index b725a587..00000000 --- a/tests/test-mul-mat1.c +++ /dev/null @@ -1,312 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -#include - -const int M = 1280; -const int N = 1536; -const int K = 1280; - -uint64_t get_time_us(void) { - struct timeval tv; - gettimeofday(&tv, NULL); - return tv.tv_sec * 1000000 + tv.tv_usec; -} - -// -// naive implementation -// - -void mul_mat_f32_0( - const float * restrict src0, // M x K - const float * restrict src1, // N x K (transposed) - float * dst, - int m, int n, int k) { - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - float sum = 0; - for (int l = 0; l < k; l++) { - sum += src0[i*k + l] * src1[j*k + l]; - } - dst[i*n + j] = sum; - } - } -} - -void mul_mat_f16_0( - const __fp16 * src0, - const __fp16 * src1, - float * dst, - int m, int n, int k) { - const int k32 = k & ~31; - - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - float sumf = 0.0; - - float16x8_t sum0 = vdupq_n_f16(0.0f); - float16x8_t sum1 = vdupq_n_f16(0.0f); - float16x8_t sum2 = vdupq_n_f16(0.0f); - float16x8_t sum3 = vdupq_n_f16(0.0f); - - float16x8_t x0, x1, x2, x3; - float16x8_t y0, y1, y2, y3; - - const __fp16 * restrict p0 = src0 + i*k; - const __fp16 * restrict p1 = src1 + j*k; - - for (int l = 0; l < k32; l += 32) { - x0 = vld1q_f16(p0 + l + 0 ); - x1 = vld1q_f16(p0 + l + 8 ); - x2 = vld1q_f16(p0 + l + 16); - x3 = vld1q_f16(p0 + l + 24); - - y0 = vld1q_f16(p1 + l + 0 ); - y1 = vld1q_f16(p1 + l + 8 ); - y2 = vld1q_f16(p1 + l + 16); - y3 = vld1q_f16(p1 + l + 24); - - sum0 = vfmaq_f16(sum0, x0, y0); - sum1 = vfmaq_f16(sum1, x1, y1); - sum2 = vfmaq_f16(sum2, x2, y2); - sum3 = vfmaq_f16(sum3, x3, y3); - } - - // reduce sum0..sum3 to sum0 - sum0 = vaddq_f16(sum0, sum1); - sum2 = vaddq_f16(sum2, sum3); - sum0 = vaddq_f16(sum0, sum2); - - // load sum0 into 2 float32x4_t - float32x4_t sum0f32 = vcvt_f32_f16(vget_low_f16(sum0)); - float32x4_t sum1f32 = vcvt_f32_f16(vget_high_f16(sum0)); - - // reduce sum0f32 and sum1f32 to sumf - sum0f32 = vaddq_f32(sum0f32, sum1f32); - - float32x2_t sumf32 = vadd_f32(vget_low_f32(sum0f32), vget_high_f32(sum0f32)); - sumf = vget_lane_f32(sumf32, 0) + vget_lane_f32(sumf32, 1); - - //sumf = sum0[0] + sum0[1] + sum0[2] + sum0[3] + sum0[4] + sum0[5] + sum0[6] + sum0[7]; - - for (int l = k32; l < k32; l++) { - sumf += p0[l]*p1[l]; - } - - dst[i*n + j] = sumf; - } - } -} - -// blocking with block size 32 -void mul_mat_f16_1( - const __fp16 * src0, - const __fp16 * src1, - float * dst, - int m, int n, int k) { - - const int k32 = k & ~31; - const int bs = 32; - - memset(dst, 0, m*n*sizeof(float)); - - for (int i = 0; i < m; i += bs) { - for (int j = 0; j < n; j += bs) { - for (int l = 0; l < k; l += bs) { - for (int ii = i; ii < i + bs; ii++) { - const __fp16 * restrict p0 = src0 + ii*k; - - float16x8_t x0, x1, x2, x3; - - x0 = vld1q_f16(p0 + l + 0 ); - x1 = vld1q_f16(p0 + l + 8 ); - x2 = vld1q_f16(p0 + l + 16); - x3 = vld1q_f16(p0 + l + 24); - - for (int jj = j; jj < j + bs; jj++) { - float sumf = 0.0; - - float16x8_t sum0 = vdupq_n_f16(0.0f); - float16x8_t sum1 = vdupq_n_f16(0.0f); - float16x8_t sum2 = vdupq_n_f16(0.0f); - float16x8_t sum3 = vdupq_n_f16(0.0f); - - float16x8_t y0, y1, y2, y3; - - const __fp16 * restrict p1 = src1 + jj*k; - - y0 = vld1q_f16(p1 + l + 0 ); - y1 = vld1q_f16(p1 + l + 8 ); - y2 = vld1q_f16(p1 + l + 16); - y3 = vld1q_f16(p1 + l + 24); - - sum0 = vfmaq_f16(sum0, x0, y0); - sum1 = vfmaq_f16(sum1, x1, y1); - sum2 = vfmaq_f16(sum2, x2, y2); - sum3 = vfmaq_f16(sum3, x3, y3); - - // reduce sum0..sum3 to sum0 - sum0 = vaddq_f16(sum0, sum1); - sum2 = vaddq_f16(sum2, sum3); - sum0 = vaddq_f16(sum0, sum2); - - // load sum0 into 2 float32x4_t - float32x4_t sum0f32 = vcvt_f32_f16(vget_low_f16(sum0)); - float32x4_t sum1f32 = vcvt_f32_f16(vget_high_f16(sum0)); - - // reduce sum0f32 and sum1f32 to sumf - sum0f32 = vaddq_f32(sum0f32, sum1f32); - - float32x2_t sumf32 = vadd_f32(vget_low_f32(sum0f32), vget_high_f32(sum0f32)); - sumf = vget_lane_f32(sumf32, 0) + vget_lane_f32(sumf32, 1); - - //sumf = sum0[0] + sum0[1] + sum0[2] + sum0[3] + sum0[4] + sum0[5] + sum0[6] + sum0[7]; - - dst[ii*n + jj] += sumf; - } - } - } - } - } - -} - -void mul_mat_f8_0( - const uint8_t * src0, - const uint8_t * src1, - float * dst, - int m, int n, int k) { - const int k32 = k & ~31; - - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - float sumf = 0.0; - - const uint8_t * restrict p0 = src0 + i*k; - const uint8_t * restrict p1 = src1 + j*k; - - for (int l = 0; l < k32; l += 32) { - uint8x16_t x0 = vld1q_u8(p0 + l + 0 ); - uint8x16_t x1 = vld1q_u8(p0 + l + 16); - - uint8x16_t y0 = vld1q_u8(p1 + l + 0 ); - uint8x16_t y1 = vld1q_u8(p1 + l + 16); - - x0 = vmulq_u8(x0, y0); - x1 = vmulq_u8(x1, y1); - - sumf += vaddvq_u8(x0) + vaddvq_u8(x1); - } - - dst[i*n + j] = sumf; - } - } -} - -int main(int argc, const char ** argv) { - float * src0 = malloc(sizeof(float)*M*K); - float * src1 = malloc(sizeof(float)*N*K); - float * dst = malloc(sizeof(float)*M*N); - - for (int i = 0; i < M*K; i++) { - src0[i] = rand() / (float)RAND_MAX; - } - - for (int i = 0; i < N*K; i++) { - src1[i] = rand() / (float)RAND_MAX; - } - - // convert src0 and src1 to __fp16 - __fp16 * src0_fp16 = (__fp16 *)(malloc(sizeof(__fp16)*M*K)); - __fp16 * src1_fp16 = (__fp16 *)(malloc(sizeof(__fp16)*N*K)); - - uint8_t * src0_fp8 = (uint8_t *)(malloc(sizeof(__fp16)*M*K)); - uint8_t * src1_fp8 = (uint8_t *)(malloc(sizeof(__fp16)*N*K)); - - { - const uint64_t t_start = get_time_us(); - - for (int i = 0; i < M*K; i++) { - src0_fp16[i] = src0[i]; - //printf("%f %f\n", src0[i], src0_fp16[i]); - //assert(!isnan(src0_fp16[i])); - } - - for (int i = 0; i < N*K; i++) { - src1_fp16[i] = src1[i]; - } - - const uint64_t t_end = get_time_us(); - printf("convert time: %f ms\n", (t_end - t_start) / 1000.0); - } - - for (int i = 0; i < 16; ++i) { - printf("%f %f\n", src0[i], src0_fp16[i]); - } - - int method = 0; - if (argc > 1) { - method = atoi(argv[1]); - } - - const int nIter = 1; - - const clock_t start = clock(); - const uint64_t start_us = get_time_us(); - - double iM = 1.0/M; - double sum = 0.0f; - for (int i = 0; i < nIter; i++) { - if (method == 0) { - mul_mat_f32_0(src0, src1, dst, M, N, K); - } - - if (method == 1) { - mul_mat_f16_0(src0_fp16, src1_fp16, dst, M, N, K); - } - - if (method == 2) { - mul_mat_f16_1(src0_fp16, src1_fp16, dst, M, N, K); - } - - if (method == 3) { - mul_mat_f8_0(src0_fp8, src1_fp8, dst, M, N, K); - } - - if (method == 4) { - // Use BLAS sgemm from Accelerate framework - cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, M, N, K, 1.0f, src0, K, src1, K, 0.0f, dst, N); - } - } - - for (int i = 0; i < N; i++) { - sum += dst[i]*iM; - } - - { - const clock_t end = clock(); - const uint64_t end_us = get_time_us(); - printf("%s: elapsed ticks: %ld\n", __func__, end - start); - printf("%s: elapsed us: %llu / %f ms\n", __func__, end_us - start_us, (end_us - start_us) / 1000.0 / nIter); - } - - printf("%f\n", sum); - - free(src0); - free(src1); - free(dst); - - free(src0_fp16); - free(src1_fp16); - - return 0; -} diff --git a/tests/test-mul-mat2.c b/tests/test-mul-mat2.c deleted file mode 100644 index 89af2863..00000000 --- a/tests/test-mul-mat2.c +++ /dev/null @@ -1,2585 +0,0 @@ -// quantized matrix multiplication - -#include "ggml.h" - -#include -#include -#include -#include -#include -#include -#include -#include - -#if defined(__ARM_NEON) -#include "arm_neon.h" -#elif defined(__AVX__) || defined(__AVX2__) -#include "immintrin.h" -#endif - -#ifndef MIN -#define MAX(a, b) ((a) > (b) ? (a) : (b)) -#define MIN(a, b) ((a) < (b) ? (a) : (b)) -#endif - -#if defined(_MSC_VER) -#pragma warning(disable: 4244 4267) // possible loss of data -#include -#define __builtin_popcountll __popcnt64 -#endif - -const int M = 1280; -const int N = 1536; -const int K = 1280; - -//const int M = 64; -//const int N = 64; -//const int K = 64; - -#define QK 64 -#define QB 4 - -//#define GGML_GQ_USE_FP16_SCALE - -#if defined(GGML_GQ_USE_FP16_SCALE) -#define gq_scale_t ggml_fp16_t -#define GGML_FP32_TO_GQ(x) ggml_fp32_to_fp16(x) -#define GGML_GQ_TO_FP32(x) ggml_fp16_to_fp32(x) -#else -#define gq_scale_t float -#define GGML_FP32_TO_GQ(x) (x) -#define GGML_GQ_TO_FP32(x) (x) -#endif - -#define gq_t_bits 64 -#define gq_quant_t uint64_t - -float frand(void) { - return (float) rand() / (float) RAND_MAX; -} - -#if defined(__AVX2__) -// horizontally reduce 8 32-bit integers -static inline uint32_t _mm256_hadd_epi32_gg(__m256i v) { - __m128i v0 = _mm256_extractf128_si256(v, 0); - __m128i v1 = _mm256_extractf128_si256(v, 1); - - v0 = _mm_add_epi32(v0, v1); - - v1 = _mm_shuffle_epi32(v0, 0x0e); - v0 = _mm_add_epi32(v0, v1); - - v1 = _mm_shuffle_epi32(v0, 0x01); - v0 = _mm_add_epi32(v0, v1); - - return _mm_cvtsi128_si32(v0); -} - -//static inline float _mm256_hadd_epi32_gg(__m256i v) { -// const __m256 v0 = _mm256_cvtepi32_ps(v); -// const __m128 t0 = _mm_add_ps(_mm256_castps256_ps128(v0), _mm256_extractf128_ps(v0, 1)); -// const __m128 t1 = _mm_hadd_ps(t0, t0); -// -// return _mm_cvtss_f32(_mm_hadd_ps(t1, t1)); -//} - -// horizontally reduce 32 8-bit integers -static inline int32_t _mm256_hadd_epi8_gg(__m256i v0) { - __m256i v1 = _mm256_maddubs_epi16(v0, _mm256_set1_epi8(1)); - __m256i v2 = _mm256_madd_epi16 (v1, _mm256_set1_epi16(1)); - - return _mm256_hadd_epi32_gg(v2); -} - -static inline float _mm256_hadd_ps_gg(__m256 v) { - const __m128 t0 = _mm_add_ps(_mm256_castps256_ps128(v), _mm256_extractf128_ps(v, 1)); - const __m128 t1 = _mm_hadd_ps(t0, t0); - - return _mm_cvtss_f32(_mm_hadd_ps(t1, t1)); -} -#endif - -// -// naive implementation -// - -void mul_mat_f32_naive( - const float * restrict src0, // M x K - const float * restrict src1, // N x K (transposed) - float * dst, - int m, int n, int k) { - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - float sum = 0; - for (int l = 0; l < k; l++) { - sum += src0[i*k + l] * src1[j*k + l]; - } - dst[i*n + j] = sum; - } - } -} - -// -// method 1 -// - -static inline int quantize_1_blocks_per_row(int k) { - return k/QK; -} - -static inline int quantize_1_quants_per_block(void) { - return QK/gq_t_bits; -} - -static inline int quantize_1_row_size(int k) { - const int nb = quantize_1_blocks_per_row(k); - const int nq = quantize_1_quants_per_block(); - - return nb*(2*sizeof(gq_scale_t) + nq*QB*sizeof(gq_quant_t)); -} - -void quantize_1(const float * src, void * dst, int n, int k) { - char * p0 = dst; - - gq_quant_t pp[QB]; - - for (int j = 0; j < n; j++) { - for (int i = 0; i < k/QK; i++) { - float min = FLT_MAX; - float max = -FLT_MAX; - - // find min/max -#ifdef __ARM_NEON - { - float32x4_t minv = vdupq_n_f32(FLT_MAX); - float32x4_t maxv = vdupq_n_f32(-FLT_MAX); - - for (int l = 0; l < QK; l += 4) { - float32x4_t v = vld1q_f32(src + j*k + i*QK + l); - minv = vminq_f32(minv, v); - maxv = vmaxq_f32(maxv, v); - } - - float32x2_t minv32 = vpmin_f32(vget_low_f32(minv), vget_high_f32(minv)); - float32x2_t maxv32 = vpmax_f32(vget_low_f32(maxv), vget_high_f32(maxv)); - - min = MIN(vget_lane_f32(minv32, 0), vget_lane_f32(minv32, 1)); - max = MAX(vget_lane_f32(maxv32, 0), vget_lane_f32(maxv32, 1)); - - //printf("SIMD min/max: %f %f\n", min, max); - } -#else - { - for (int l = 0; l < QK; l++) { - const float v = src[j*k + i*QK + l]; - if (v < min) min = v; - if (v > max) max = v; - } - - //printf("NORM min/max: %f %f\n", min, max); - } -#endif - - const float d = (max - min) / ((1 << QB) - 1); - const float id = d ? 1.0/d : 0.0; - - memcpy(p0, &min, sizeof(float)); p0 += sizeof(float); - memcpy(p0, &d, sizeof(float)); p0 += sizeof(float); - - //printf("min/max/d/id: %f %f %f %f\n", min, max, d, id); - - for (int s = 0; s < QK/gq_t_bits; ++s) { - memset(pp, 0, sizeof(pp)); - - for (int l = 0; l < gq_t_bits; l++) { - const float v = src[j*k + i*QK + s*gq_t_bits + l]; - const uint8_t q = (v - min)*id; - - for (int b = 0; b < QB; b++) { - pp[b] |= q & (1 << b) ? (1ULL << l) : 0; - } - } - - for (int b = 0; b < QB; b++) { - memcpy(p0, &pp[b], sizeof(gq_quant_t)); p0 += sizeof(gq_quant_t); - } - } - } - } -} - -void mul_mat_gq_1( - const void * src0, - const void * src1, - float * dst, - int m, int n, int k) { - const int kp = k & ~(gq_t_bits - 1); - - const char * restrict p0 = src0; - const char * restrict p1 = src1; - - float s0[QB + 1]; - float s1[QB + 1]; - - gq_quant_t m0[QB + 1]; - gq_quant_t m1[QB + 1]; - - for (int ir0 = 0; ir0 < m; ir0++) { - for (int ir1 = 0; ir1 < n; ir1++) { - float sumf = 0.0; - - const char * restrict pp0 = p0 + ir0*((2*sizeof(float) + (QK/gq_t_bits)*QB*sizeof(gq_quant_t))*(k/QK)); - const char * restrict pp1 = p1 + ir1*((2*sizeof(float) + (QK/gq_t_bits)*QB*sizeof(gq_quant_t))*(k/QK)); - - for (int i = 0; i < kp/QK; i++) { - float min0, d0; - memcpy(&min0, pp0, sizeof(float)); pp0 += sizeof(float); - memcpy(&d0, pp0, sizeof(float)); pp0 += sizeof(float); - - float min1, d1; - memcpy(&min1, pp1, sizeof(float)); pp1 += sizeof(float); - memcpy(&d1, pp1, sizeof(float)); pp1 += sizeof(float); - - //printf("min0/d0 = %f %f | min1/d1 = %f %f\n", min0, d0, min1, d1); - -#if 1 - // >>> General case for any QB - - s0[0] = min0; - s1[0] = min1; - - for (int b = 0; b < QB; b++) { - s0[b + 1] = d0*(1 << b); - s1[b + 1] = d1*(1 << b); - } - - m0[0] = 0-1ULL; - m1[0] = 0-1ULL; - - for (int s = 0; s < QK/gq_t_bits; ++s) { - for (int b = 0; b < QB; b++) { - memcpy(&m0[b + 1], pp0, sizeof(gq_quant_t)); pp0 += sizeof(gq_quant_t); - memcpy(&m1[b + 1], pp1, sizeof(gq_quant_t)); pp1 += sizeof(gq_quant_t); - } - - for (int q0 = 0; q0 < QB + 1; q0++) { - for (int q1 = 0; q1 < QB + 1; q1++) { - sumf += s0[q0]*s1[q1]*__builtin_popcountll(m0[q0] & m1[q1]); - } - } - } -#else -#endif - } - - dst[ir0*n + ir1] = sumf; - } - } -} - -// -// method 2 -// n-bit quantization (2nd attempt) -// - -static inline int quantize_2_blocks_per_row(int k) { - return k/QK; -} - -static inline int quantize_2_quants_per_block(void) { - return QK/gq_t_bits; -} - -static inline int quantize_2_row_size(int k) { - const int nb = quantize_2_blocks_per_row(k); - const int nq = quantize_2_quants_per_block(); - - return nb*(2*sizeof(gq_scale_t) + nq*QB*sizeof(gq_quant_t)); -} - -void quantize_2_row(const float * restrict src, void * restrict dst, int k) { - assert(k % QK == 0); - - const int nb = quantize_2_blocks_per_row(k); - const int nq = quantize_2_quants_per_block(); - - gq_scale_t * restrict pm = (gq_scale_t *) (dst); - gq_scale_t * restrict pd = (gq_scale_t *) (pm + nb); - gq_quant_t * restrict pb = (gq_quant_t *) (pd + nb); - - gq_quant_t pp[QB]; - - static const int32_t sh[32] = { - 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, - 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, - }; - - for (int i = 0; i < nb; i++) { - float min = FLT_MAX; - float max = -FLT_MAX; - -#ifdef __ARM_NEON - { - float32x4_t minv = vdupq_n_f32(FLT_MAX); - float32x4_t maxv = vdupq_n_f32(-FLT_MAX); - - for (int l = 0; l < QK; l += 4) { - float32x4_t v = vld1q_f32(src + i*QK + l); - minv = vminq_f32(minv, v); - maxv = vmaxq_f32(maxv, v); - } - - float32x2_t minv32 = vpmin_f32(vget_low_f32(minv), vget_high_f32(minv)); - float32x2_t maxv32 = vpmax_f32(vget_low_f32(maxv), vget_high_f32(maxv)); - - min = MIN(vget_lane_f32(minv32, 0), vget_lane_f32(minv32, 1)); - max = MAX(vget_lane_f32(maxv32, 0), vget_lane_f32(maxv32, 1)); - } -#else - { - for (int l = 0; l < QK; l++) { - const float v = src[i*QK + l]; - if (v < min) min = v; - if (v > max) max = v; - } - } -#endif - - const float d = (max - min) / ((1 << QB) - 1); - const float id = d ? 1.0/d : 0.0; - - pm[i] = GGML_FP32_TO_GQ(min); - pd[i] = GGML_FP32_TO_GQ(d); - - for (int s = 0; s < nq; ++s) { - memset(pp, 0, sizeof(pp)); - -#if 1 - for (int l = 0; l < gq_t_bits; l++) { - const float v = src[i*QK + s*gq_t_bits + l]; - const uint8_t q = (v - min)*id + frand(); - - for (int b = 0; b < QB; b++) { - pp[b] |= q & (1 << b) ? (1ULL << l) : 0; - } - } -#elif defined(__ARM_NEON) -#if 1 - { - uint32_t ppt[2*4*QB]; - - float32x4_t minv = vdupq_n_f32(min); - float32x4_t idv = vdupq_n_f32(id); - - assert(gq_t_bits % 16 == 0); - - uint32x4_t p0[QB] = { vdupq_n_u32(0) }; - uint32x4_t p1[QB] = { vdupq_n_u32(0) }; - - for (int l = 0; l < gq_t_bits; l += 16) { - float32x4_t v0 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 0); - float32x4_t v1 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 4); - float32x4_t v2 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 8); - float32x4_t v3 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 12); - - v0 = vsubq_f32(v0, minv); - v1 = vsubq_f32(v1, minv); - v2 = vsubq_f32(v2, minv); - v3 = vsubq_f32(v3, minv); - - v0 = vmulq_f32(v0, idv); - v1 = vmulq_f32(v1, idv); - v2 = vmulq_f32(v2, idv); - v3 = vmulq_f32(v3, idv); - -#if 1 - v0[0] += frand(); v0[1] += frand(); v0[2] += frand(); v0[3] += frand(); - v1[0] += frand(); v1[1] += frand(); v1[2] += frand(); v1[3] += frand(); - v2[0] += frand(); v2[1] += frand(); v2[2] += frand(); v2[3] += frand(); - v3[0] += frand(); v3[1] += frand(); v3[2] += frand(); v3[3] += frand(); -#endif - - uint32x4_t q0 = vcvtq_u32_f32(v0); - uint32x4_t q1 = vcvtq_u32_f32(v1); - uint32x4_t q2 = vcvtq_u32_f32(v2); - uint32x4_t q3 = vcvtq_u32_f32(v3); - - for (int b = 0; b < QB; ++b) { - uint32x4_t m = vdupq_n_u32(1 << b); - uint32x4_t r = vdupq_n_u32(-b); - - if (l < 32) { - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q0, m), r), vld1q_s32(sh + l + 0))); - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q1, m), r), vld1q_s32(sh + l + 4))); - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q2, m), r), vld1q_s32(sh + l + 8))); - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q3, m), r), vld1q_s32(sh + l + 12))); - } else { - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q0, m), r), vld1q_s32(sh + l - 32))); - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q1, m), r), vld1q_s32(sh + l - 28))); - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q2, m), r), vld1q_s32(sh + l - 24))); - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q3, m), r), vld1q_s32(sh + l - 20))); - } - } - } - -#if QB == 4 - vst1q_u32((uint32_t *) ppt + 0, p0[0]); - vst1q_u32((uint32_t *) ppt + 4, p1[0]); - vst1q_u32((uint32_t *) ppt + 8, p0[1]); - vst1q_u32((uint32_t *) ppt + 12, p1[1]); - vst1q_u32((uint32_t *) ppt + 16, p0[2]); - vst1q_u32((uint32_t *) ppt + 20, p1[2]); - vst1q_u32((uint32_t *) ppt + 24, p0[3]); - vst1q_u32((uint32_t *) ppt + 28, p1[3]); - - pp[0] = (ppt[0] | ppt[1] | ppt[2] | ppt[3] ) | ((uint64_t) (ppt[4] | ppt[5] | ppt[6] | ppt[7]) ) << 32; - pp[1] = (ppt[8] | ppt[9] | ppt[10] | ppt[11]) | ((uint64_t) (ppt[12] | ppt[13] | ppt[14] | ppt[15])) << 32; - pp[2] = (ppt[16] | ppt[17] | ppt[18] | ppt[19]) | ((uint64_t) (ppt[20] | ppt[21] | ppt[22] | ppt[23])) << 32; - pp[3] = (ppt[24] | ppt[25] | ppt[26] | ppt[27]) | ((uint64_t) (ppt[28] | ppt[29] | ppt[30] | ppt[31])) << 32; -#else - for (int b = 0; b < QB; ++b) { - vst1q_u32((uint32_t *) ppt + 0, p0[b]); - vst1q_u32((uint32_t *) ppt + 4, p1[b]); - - pp[b] = (ppt[0] | ppt[1] | ppt[2] | ppt[3]) | ((uint64_t) (ppt[4] | ppt[5] | ppt[6] | ppt[7])) << 32; - } -#endif - } -#else - // less optimal SIMD - { - float32x4_t minv = vdupq_n_f32(min); - float32x4_t idv = vdupq_n_f32(id); - - assert(gq_t_bits == 64); - uint8_t qq[gq_t_bits]; - - for (int l = 0; l < gq_t_bits; l += 16) { - float32x4_t v0 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 0); - float32x4_t v1 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 4); - float32x4_t v2 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 8); - float32x4_t v3 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 12); - - v0 = vsubq_f32(v0, minv); - v1 = vsubq_f32(v1, minv); - v2 = vsubq_f32(v2, minv); - v3 = vsubq_f32(v3, minv); - - v0 = vmulq_f32(v0, idv); - v1 = vmulq_f32(v1, idv); - v2 = vmulq_f32(v2, idv); - v3 = vmulq_f32(v3, idv); - -#if 0 - v0[0] += frand(); v0[1] += frand(); v0[2] += frand(); v0[3] += frand(); - v1[0] += frand(); v1[1] += frand(); v1[2] += frand(); v1[3] += frand(); - v2[0] += frand(); v2[1] += frand(); v2[2] += frand(); v2[3] += frand(); - v3[0] += frand(); v3[1] += frand(); v3[2] += frand(); v3[3] += frand(); -#endif - - uint32x4_t q0 = vcvtq_u32_f32(v0); - uint32x4_t q1 = vcvtq_u32_f32(v1); - uint32x4_t q2 = vcvtq_u32_f32(v2); - uint32x4_t q3 = vcvtq_u32_f32(v3); - - // store in qq as uint8_t - vst1_u8(qq + l + 0, vmovn_u16(vcombine_u16(vmovn_u32(q0), vmovn_u32(q1)))); - vst1_u8(qq + l + 8, vmovn_u16(vcombine_u16(vmovn_u32(q2), vmovn_u32(q3)))); - } - - for (int l = 0; l < gq_t_bits; l++) { - for (int b = 0; b < QB; b++) { - const uint64_t ql = qq[l]; - /*pp[b] |= qq[l] & (1 << b) ? (1ULL << l) : 0;*/ - pp[b] |= ((ql & (1 << b)) >> b) << l; - } - } - } -#endif -#endif - memcpy(pb + i*nq*QB + s*QB, pp, sizeof(pp)); - } - } -} - -// reimplementation of quantize_2 using quantize_2_row -void quantize_2(const float * restrict src, char * restrict dst, int n, int k) { - assert(k % QK == 0); - - for (int j = 0; j < n; j++) { - quantize_2_row(src + j*k, dst, k); - dst = (char *) dst + quantize_2_row_size(k); - } -} - -void vec_dot_gq_2(const int n, float * restrict s, const void * restrict x, const void * restrict y) { - const int nb = quantize_2_blocks_per_row(n); - const int nq = quantize_2_quants_per_block(); - - const gq_scale_t * restrict pm0 = (const gq_scale_t *) x; - const gq_scale_t * restrict pm1 = (const gq_scale_t *) y; - - const gq_scale_t * restrict pd0 = pm0 + nb; - const gq_scale_t * restrict pd1 = pm1 + nb; - - const gq_quant_t * restrict pb0 = (const gq_quant_t *) (pd0 + nb); - const gq_quant_t * restrict pb1 = (const gq_quant_t *) (pd1 + nb); - - float sumf = 0.0; - -#if 1 - for (int i = 0; i < nb; i++) { - const float m0 = GGML_GQ_TO_FP32(pm0[i]); - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - - const float m1 = GGML_GQ_TO_FP32(pm1[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - -#if QB == 4 - int isum01 = 0; - int isum10 = 0; - int isum11 = 0; - - for (int s = 0; s < nq; ++s) { - const gq_quant_t * restrict mm0 = pb0 + i*nq*QB + s*QB; - const gq_quant_t * restrict mm1 = pb1 + i*nq*QB + s*QB; - -#define bpcnt(x) __builtin_popcountll(x) - isum01 += (1 << 0)*(bpcnt(mm1[0])); - isum01 += (1 << 1)*(bpcnt(mm1[1])); - isum01 += (1 << 2)*(bpcnt(mm1[2])); - isum01 += (1 << 3)*(bpcnt(mm1[3])); - - isum10 += (1 << 0)*(bpcnt(mm0[0])); - isum10 += (1 << 1)*(bpcnt(mm0[1])); - isum10 += (1 << 2)*(bpcnt(mm0[2])); - isum10 += (1 << 3)*(bpcnt(mm0[3])); - - isum11 += (1 << 0)*(bpcnt(mm0[0] & mm1[0])); - isum11 += (1 << 1)*(bpcnt(mm0[0] & mm1[1]) + bpcnt(mm0[1] & mm1[0])); - isum11 += (1 << 2)*(bpcnt(mm0[0] & mm1[2]) + bpcnt(mm0[1] & mm1[1]) + bpcnt(mm0[2] & mm1[0])); - isum11 += (1 << 3)*(bpcnt(mm0[0] & mm1[3]) + bpcnt(mm0[1] & mm1[2]) + bpcnt(mm0[2] & mm1[1]) + bpcnt(mm0[3] & mm1[0])); - isum11 += (1 << 4)*(bpcnt(mm0[1] & mm1[3]) + bpcnt(mm0[2] & mm1[2]) + bpcnt(mm0[3] & mm1[1])); - isum11 += (1 << 5)*(bpcnt(mm0[2] & mm1[3]) + bpcnt(mm0[3] & mm1[2])); - isum11 += (1 << 6)*(bpcnt(mm0[3] & mm1[3])); -#undef bpcnt - } - - sumf += nq*gq_t_bits*(m0*m1) + isum01*(m0*d1) + isum10*(m1*d0) + isum11*(d0*d1); -#elif QB == 3 - int isum01 = 0; - int isum10 = 0; - int isum11 = 0; - - for (int s = 0; s < nq; ++s) { - const gq_quant_t * restrict mm0 = pb0 + i*nq*QB + s*QB; - const gq_quant_t * restrict mm1 = pb1 + i*nq*QB + s*QB; - -#if gq_t_bits == 32 -#define bpcnt(x) __builtin_popcount(x) -#else -#define bpcnt(x) __builtin_popcountll(x) -#endif - isum01 += (1 << 0)*(bpcnt(mm1[0])); - isum01 += (1 << 1)*(bpcnt(mm1[1])); - isum01 += (1 << 2)*(bpcnt(mm1[2])); - - isum10 += (1 << 0)*(bpcnt(mm0[0])); - isum10 += (1 << 1)*(bpcnt(mm0[1])); - isum10 += (1 << 2)*(bpcnt(mm0[2])); - - isum11 += (1 << 0)*(bpcnt(mm0[0] & mm1[0])); - isum11 += (1 << 1)*(bpcnt(mm0[0] & mm1[1]) + bpcnt(mm0[1] & mm1[0])); - isum11 += (1 << 2)*(bpcnt(mm0[0] & mm1[2]) + bpcnt(mm0[1] & mm1[1]) + bpcnt(mm0[2] & mm1[0])); - isum11 += (1 << 3)*(bpcnt(mm0[1] & mm1[2]) + bpcnt(mm0[2] & mm1[1])); - isum11 += (1 << 4)*(bpcnt(mm0[2] & mm1[2])); -#undef bpcnt - } - - sumf += nq*gq_t_bits*(m0*m1) + isum01*(m0*d1) + isum10*(m1*d0) + isum11*(d0*d1); -#elif QB == 2 - int isum01 = 0; - int isum10 = 0; - int isum11 = 0; - - for (int s = 0; s < nq; ++s) { - const gq_quant_t * restrict mm0 = pb0 + i*nq*QB + s*QB; - const gq_quant_t * restrict mm1 = pb1 + i*nq*QB + s*QB; - -#if gq_t_bits == 32 -#define bpcnt(x) __builtin_popcount(x) -#else -#define bpcnt(x) __builtin_popcountll(x) -#endif - isum01 += (1 << 0)*(bpcnt(mm1[0])); - isum01 += (1 << 1)*(bpcnt(mm1[1])); - - isum10 += (1 << 0)*(bpcnt(mm0[0])); - isum10 += (1 << 1)*(bpcnt(mm0[1])); - - isum11 += (1 << 0)*(bpcnt(mm0[0] & mm1[0])); - isum11 += (1 << 1)*(bpcnt(mm0[0] & mm1[1]) + bpcnt(mm0[1] & mm1[0])); - isum11 += (1 << 2)*(bpcnt(mm0[1] & mm1[1])); -#undef bpcnt - } - - sumf += nq*gq_t_bits*(m0*m1) + isum01*(m0*d1) + isum10*(m1*d0) + isum11*(d0*d1); -#else - float s0[QB + 1]; - float s1[QB + 1]; - - s0[0] = m0; - s1[0] = m1; - - for (int b = 0; b < QB; b++) { - s0[b + 1] = d0*(1 << b); - s1[b + 1] = d1*(1 << b); - } - - for (int s = 0; s < nq; ++s) { - for (int q0 = 0; q0 < QB + 1; q0++) { - const gq_quant_t mm0 = q0 ? pb0[i*nq*QB + s*QB + q0 - 1] : -1ULL; - for (int q1 = 0; q1 < QB + 1; q1++) { - const gq_quant_t mm1 = q1 ? pb1[i*nq*QB + s*QB + q1 - 1] : -1ULL; - sumf += s0[q0]*s1[q1]*__builtin_popcountll(mm0 & mm1); - } - } - } -#endif - } -#else -#error "not implemented" -#endif - - *s = sumf; -} - -// use vec_dot_gq_2 to compute the dot product of two rows -void mul_mat_gq_2( - const void * src0, - const void * src1, // transposed - float * dst, - int m, int n, int k) { - assert(k % QK == 0); - - for (int ir0 = 0; ir0 < m; ir0++) { - for (int ir1 = 0; ir1 < n; ir1++) { - vec_dot_gq_2(k, dst + ir1, src0, src1); - src1 = (const char *) src1 + quantize_2_row_size(k); - } - src0 = (const char *) src0 + quantize_2_row_size(k); - src1 = (const char *) src1 - n*quantize_2_row_size(k); - - dst = (float *) dst + n; - } -} - -// -// method 3 -// (does not work) -// - -static inline int quantize_3_blocks_per_row(int k) { - return k/QK; -} - -static inline int quantize_3_quants_per_block(void) { - return QK/gq_t_bits; -} - -static inline int quantize_3_row_size(int k) { - const int nb = quantize_3_blocks_per_row(k); - const int nq = quantize_3_quants_per_block(); - - return nb*(sizeof(gq_scale_t) + nq*QB*sizeof(gq_quant_t)); -} - -void quantize_3_row(const float * restrict src, void * restrict dst, int k) { - assert(k % QK == 0); - - const int nb = quantize_3_blocks_per_row(k); - const int nq = quantize_3_quants_per_block(); - - gq_scale_t * restrict pd = (gq_scale_t *) (dst); - gq_quant_t * restrict pb = (gq_quant_t *) (pd + nb); - - gq_quant_t pp[QB]; - - static const int32_t sh[32] = { - 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, - 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, - }; - - for (int i = 0; i < nb; i++) { - float amax = 0.0f; // abs max - -#ifdef __ARM_NEON - { - // min / max - //float32x4_t minv = vdupq_n_f32(FLT_MAX); - //float32x4_t maxv = vdupq_n_f32(-FLT_MAX); - - //for (int l = 0; l < QK; l += 4) { - // float32x4_t v = vld1q_f32(src + i*QK + l); - // minv = vminq_f32(minv, v); - // maxv = vmaxq_f32(maxv, v); - //} - - //float32x2_t minv32 = vpmin_f32(vget_low_f32(minv), vget_high_f32(minv)); - //float32x2_t maxv32 = vpmax_f32(vget_low_f32(maxv), vget_high_f32(maxv)); - - //min = MIN(vget_lane_f32(minv32, 0), vget_lane_f32(minv32, 1)); - //max = MAX(vget_lane_f32(maxv32, 0), vget_lane_f32(maxv32, 1)); - - // abs max - float32x4_t amaxv = vdupq_n_f32(0.0f); - - for (int l = 0; l < QK; l += 4) { - float32x4_t v = vld1q_f32(src + i*QK + l); - amaxv = vmaxq_f32(amaxv, vabsq_f32(v)); - } - - float32x2_t amaxv32 = vpmax_f32(vget_low_f32(amaxv), vget_high_f32(amaxv)); - - amax = MAX(vget_lane_f32(amaxv32, 0), vget_lane_f32(amaxv32, 1)); - } -#else - { - for (int l = 0; l < QK; l++) { - const float v = src[i*QK + l]; - amax = MAX(amax, fabsf(v)); - } - } -#endif - - const float d = amax / ((1 << (QB - 1)) - 1); - const float id = d ? 1.0/d : 0.0; - - pd[i] = GGML_FP32_TO_GQ(d); - - for (int s = 0; s < nq; ++s) { - memset(pp, 0, sizeof(pp)); - -#if 0 - for (int l = 0; l < gq_t_bits; l++) { - const float v = src[i*QK + s*gq_t_bits + l]; - const uint8_t q = v*id + frand(); - - for (int b = 0; b < QB; b++) { - pp[b] |= q & (1 << b) ? (1ULL << l) : 0; - } - } -#elif defined(__ARM_NEON) - { - uint32_t ppt[2*4*QB]; - - float32x4_t idv = vdupq_n_f32(id); - - assert(gq_t_bits == 64); - - uint32x4_t p0[QB] = { vdupq_n_u32(0) }; - uint32x4_t p1[QB] = { vdupq_n_u32(0) }; - - for (int l = 0; l < gq_t_bits; l += 16) { - float32x4_t v0 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 0); - float32x4_t v1 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 4); - float32x4_t v2 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 8); - float32x4_t v3 = vld1q_f32(src + i*QK + s*gq_t_bits + l + 12); - - v0 = vmulq_f32(v0, idv); - v1 = vmulq_f32(v1, idv); - v2 = vmulq_f32(v2, idv); - v3 = vmulq_f32(v3, idv); - -#if 1 - v0[0] += frand(); v0[1] += frand(); v0[2] += frand(); v0[3] += frand(); - v1[0] += frand(); v1[1] += frand(); v1[2] += frand(); v1[3] += frand(); - v2[0] += frand(); v2[1] += frand(); v2[2] += frand(); v2[3] += frand(); - v3[0] += frand(); v3[1] += frand(); v3[2] += frand(); v3[3] += frand(); -#endif - - uint32x4_t q0 = vcvtq_u32_f32(v0); - uint32x4_t q1 = vcvtq_u32_f32(v1); - uint32x4_t q2 = vcvtq_u32_f32(v2); - uint32x4_t q3 = vcvtq_u32_f32(v3); - - for (int b = 0; b < QB; ++b) { - uint32x4_t m = vdupq_n_u32(1 << b); - int32x4_t r = vdupq_n_s32(-b); - - if (l < 32) { - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q0, m), r), vld1q_s32(sh + l + 0))); - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q1, m), r), vld1q_s32(sh + l + 4))); - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q2, m), r), vld1q_s32(sh + l + 8))); - p0[b] = vorrq_u32(p0[b], vshlq_u32(vshlq_u32(vandq_u32(q3, m), r), vld1q_s32(sh + l + 12))); - } else { - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q0, m), r), vld1q_s32(sh + l - 32))); - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q1, m), r), vld1q_s32(sh + l - 28))); - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q2, m), r), vld1q_s32(sh + l - 24))); - p1[b] = vorrq_u32(p1[b], vshlq_u32(vshlq_u32(vandq_u32(q3, m), r), vld1q_s32(sh + l - 20))); - } - } - } - -#if QB == 4 - vst1q_u32((uint32_t *) ppt + 0, p0[0]); - vst1q_u32((uint32_t *) ppt + 4, p1[0]); - vst1q_u32((uint32_t *) ppt + 8, p0[1]); - vst1q_u32((uint32_t *) ppt + 12, p1[1]); - vst1q_u32((uint32_t *) ppt + 16, p0[2]); - vst1q_u32((uint32_t *) ppt + 20, p1[2]); - vst1q_u32((uint32_t *) ppt + 24, p0[3]); - vst1q_u32((uint32_t *) ppt + 28, p1[3]); - - pp[0] = (ppt[0] | ppt[1] | ppt[2] | ppt[3] ) | ((uint64_t) (ppt[4] | ppt[5] | ppt[6] | ppt[7]) ) << 32; - pp[1] = (ppt[8] | ppt[9] | ppt[10] | ppt[11]) | ((uint64_t) (ppt[12] | ppt[13] | ppt[14] | ppt[15])) << 32; - pp[2] = (ppt[16] | ppt[17] | ppt[18] | ppt[19]) | ((uint64_t) (ppt[20] | ppt[21] | ppt[22] | ppt[23])) << 32; - pp[3] = (ppt[24] | ppt[25] | ppt[26] | ppt[27]) | ((uint64_t) (ppt[28] | ppt[29] | ppt[30] | ppt[31])) << 32; -#else - for (int q = 0; q < QB; ++q) { - vst1q_u32((uint32_t *) ppt + 0, p0[q]); - vst1q_u32((uint32_t *) ppt + 4, p1[q]); - - pp[q] = (ppt[0] | ppt[1] | ppt[2] | ppt[3]) | ((uint64_t) (ppt[4] | ppt[5] | ppt[6] | ppt[7])) << 32; - } -#endif - } -#endif - memcpy(pb + i*nq*QB + s*QB, pp, sizeof(pp)); - } - } -} - -// reimplementation of quantize_3 using quantize_3_row -void quantize_3(const float * restrict src, char * restrict dst, int n, int k) { - assert(k % QK == 0); - - for (int j = 0; j < n; j++) { - quantize_3_row(src + j*k, dst, k); - dst = (char *) dst + quantize_3_row_size(k); - } -} - -void vec_dot_gq_3(const int n, float * restrict s, const void * restrict x, const void * restrict y) { - float sumf = 0.0f; - - const int nb = quantize_3_blocks_per_row(n); - const int nq = quantize_3_quants_per_block(); - - const gq_scale_t * restrict pd0 = (const gq_scale_t *) x; - const gq_scale_t * restrict pd1 = (const gq_scale_t *) y; - - const gq_quant_t * restrict pb0 = (const gq_quant_t *) (pd0 + nb); - const gq_quant_t * restrict pb1 = (const gq_quant_t *) (pd1 + nb); - -#if 1 - for (int i = 0; i < nb; i++) { - int isum = 0; - -#if QB == 4 - for (int s = 0; s < nq; ++s) { - const gq_quant_t * restrict m0 = pb0 + i*nq*QB + s*QB; - const gq_quant_t * restrict m1 = pb1 + i*nq*QB + s*QB; - - isum += (1 << 0)*(__builtin_popcountll(m0[0] & m1[0])); - isum += (1 << 1)*(__builtin_popcountll(m0[0] & m1[1]) + __builtin_popcountll(m0[1] & m1[0])); - isum += (1 << 2)*(__builtin_popcountll(m0[0] & m1[2]) + __builtin_popcountll(m0[1] & m1[1]) + __builtin_popcountll(m0[2] & m1[0])); - isum += (1 << 3)*(__builtin_popcountll(m0[0] & m1[3]) + __builtin_popcountll(m0[1] & m1[2]) + __builtin_popcountll(m0[2] & m1[1]) + __builtin_popcountll(m0[3] & m1[0])); - isum += (1 << 4)*(__builtin_popcountll(m0[1] & m1[3]) + __builtin_popcountll(m0[2] & m1[2]) + __builtin_popcountll(m0[3] & m1[1])); - isum += (1 << 5)*(__builtin_popcountll(m0[2] & m1[3]) + __builtin_popcountll(m0[3] & m1[2])); - isum += (1 << 6)*(__builtin_popcountll(m0[3] & m1[3])); - } -#else - for (int s = 0; s < nq; ++s) { - for (int q0 = 0; q0 < QB; q0++) { - const gq_quant_t mm0 = pb0[i*nq*QB + s*QB + q0]; - for (int q1 = 0; q1 < QB; q1++) { - const gq_quant_t mm1 = pb1[i*nq*QB + s*QB + q1]; - isum += (1 << (q0 + q1))*(__builtin_popcountll(mm0 & mm1)); - } - } - } -#endif - - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - sumf += d0*d1*isum; - } -#else -#ifdef __ARM_NEON - // gq_quant_t == uint64_t - for (int i = 0; i < nb; i += 4) { - int isum[4] = {0, 0, 0, 0}; - - for (int k = 0; k < 4; ++k) { - for (int s = 0; s < nq; ++s) { - const gq_quant_t * restrict m0 = pb0 + (i+k)*nq*QB + s*QB; - const gq_quant_t * restrict m1 = pb1 + (i+k)*nq*QB + s*QB; - -#if QB == 4 -#define bpcnt(x) __builtin_popcountll(x) - //isum[k] += (1ULL << 0)*(bpcnt(m0[0] & m1[0])) + - // (1ULL << 1)*(bpcnt(m0[0] & m1[1]) + bpcnt(m0[1] & m1[0])) + - // (1ULL << 2)*(bpcnt(m0[0] & m1[2]) + bpcnt(m0[1] & m1[1]) + bpcnt(m0[2] & m1[0])) + - // (1ULL << 3)*(bpcnt(m0[0] & m1[3]) + bpcnt(m0[1] & m1[2]) + bpcnt(m0[2] & m1[1]) + bpcnt(m0[3] & m1[0])) + - // (1ULL << 4)*(bpcnt(m0[1] & m1[3]) + bpcnt(m0[2] & m1[2]) + bpcnt(m0[3] & m1[1])) + - // (1ULL << 5)*(bpcnt(m0[2] & m1[3]) + bpcnt(m0[3] & m1[2])) + - // (1ULL << 6)*(bpcnt(m0[3] & m1[3])); -#undef bpcnt - - const uint8x8_t m00 = vld1_u8((const uint8_t *) (m0 + 0)); - const uint8x8_t m01 = vld1_u8((const uint8_t *) (m0 + 1)); - const uint8x8_t m02 = vld1_u8((const uint8_t *) (m0 + 2)); - const uint8x8_t m03 = vld1_u8((const uint8_t *) (m0 + 3)); - - const uint8x8_t m10 = vld1_u8((const uint8_t *) (m1 + 0)); - const uint8x8_t m11 = vld1_u8((const uint8_t *) (m1 + 1)); - const uint8x8_t m12 = vld1_u8((const uint8_t *) (m1 + 2)); - const uint8x8_t m13 = vld1_u8((const uint8_t *) (m1 + 3)); - - const uint8x8_t m00m10 = vand_u8(m00, m10); - - const uint8x8_t m00m11 = vand_u8(m00, m11); - const uint8x8_t m01m10 = vand_u8(m01, m10); - - const uint8x8_t m00m12 = vand_u8(m00, m12); - const uint8x8_t m01m11 = vand_u8(m01, m11); - const uint8x8_t m02m10 = vand_u8(m02, m10); - - const uint8x8_t m00m13 = vand_u8(m00, m13); - const uint8x8_t m01m12 = vand_u8(m01, m12); - const uint8x8_t m02m11 = vand_u8(m02, m11); - const uint8x8_t m03m10 = vand_u8(m03, m10); - - const uint8x8_t m01m13 = vand_u8(m01, m13); - const uint8x8_t m02m12 = vand_u8(m02, m12); - const uint8x8_t m03m11 = vand_u8(m03, m11); - - const uint8x8_t m02m13 = vand_u8(m02, m13); - const uint8x8_t m03m12 = vand_u8(m03, m12); - - const uint8x8_t m03m13 = vand_u8(m03, m13); - -#define bpcnt(x) vaddv_u8(vcnt_u8(x)) - isum[k] += (1ULL << 0)*(bpcnt(m00m10)) + - (1ULL << 1)*(bpcnt(m00m11) + bpcnt(m01m10)) + - (1ULL << 2)*(bpcnt(m00m12) + bpcnt(m01m11) + bpcnt(m02m10)) + - (1ULL << 3)*(bpcnt(m00m13) + bpcnt(m01m12) + bpcnt(m02m11) + bpcnt(m03m10)) + - (1ULL << 4)*(bpcnt(m01m13) + bpcnt(m02m12) + bpcnt(m03m11)) + - (1ULL << 5)*(bpcnt(m02m13) + bpcnt(m03m12)) + - (1ULL << 6)*(bpcnt(m03m13)); -#undef bpcnt -#else - for (int q0 = 0; q0 < QB; q0++) { - const gq_quant_t mm0 = m0[q0]; - for (int q1 = 0; q1 < QB; q1++) { - const gq_quant_t mm1 = m1[q1]; - isum[k] += (1ULL << (q0 + q1))*(__builtin_popcountll(mm0 & mm1)); - } - } -#endif - } - } - - int32x4_t isumv = vld1q_s32(isum); - - float32x4_t d0v = vld1q_f32(pd0 + i); - float32x4_t d1v = vld1q_f32(pd1 + i); - - float32x4_t sumfv = vmulq_f32(d0v, d1v); - - sumfv = vmulq_f32(sumfv, vcvtq_f32_s32(isumv)); - sumf += vaddvq_f32(sumfv); - } -#else -#error "not implemented" -#endif - -#endif - *s = sumf; -} - -// use vec_dot_gq_3 to compute the dot product of two rows -void mul_mat_gq_3( - const void * src0, - const void * src1, // transposed - float * dst, - int m, int n, int k) { - assert(k % QK == 0); - - const int nb = quantize_3_blocks_per_row(k); - const int nq = quantize_3_quants_per_block(); - - for (int ir0 = 0; ir0 < m; ir0++) { - for (int ir1 = 0; ir1 < n; ir1++) { - vec_dot_gq_3(k, dst + ir1, src0, src1); - src1 = (const char *) src1 + quantize_3_row_size(k); - } - src0 = (const char *) src0 + quantize_3_row_size(k); - src1 = (const char *) src1 - n*quantize_3_row_size(k); - - dst = (float *) dst + n; - } -} - -// -// method 4 -// 4-bit quantization -// - -static inline int quantize_4_blocks_per_row(int k) { - return k/QK; -} - -static inline int quantize_4_row_size(int k) { - const int nb = quantize_4_blocks_per_row(k); - - return nb*(2*sizeof(gq_scale_t) + QK/2); -} - -void quantize_4_row(const float * restrict src, void * restrict dst, int k) { - assert(k % QK == 0); - assert(QB == 4); - - const int nb = quantize_4_blocks_per_row(k); - - gq_scale_t * restrict pm = (gq_scale_t *) (dst); - gq_scale_t * restrict pd = (gq_scale_t *) (pm + nb); - uint8_t * restrict pb = (uint8_t *) (pd + nb); - - uint8_t pp[QK/2]; - - for (int i = 0; i < nb; i++) { - memset(pp, 0, sizeof(pp)); - - float min = FLT_MAX; - float max = -FLT_MAX; - -#if defined(__AVX2__) - { - assert(QK == 64); - enum { QK8 = QK/8 }; - - __m256 srcv[QK8]; - __m256 minv[QK8]; - __m256 maxv[QK8]; - - for (int l = 0; l < QK8; l++) { - srcv[l] = _mm256_loadu_ps(src + i*QK + 8*l); - } - - for (int l = 0; l < QK8/2; l++) { - minv[2*l] = _mm256_min_ps(srcv[2*l], srcv[2*l+1]); - maxv[2*l] = _mm256_max_ps(srcv[2*l], srcv[2*l+1]); - } - - for (int l = 0; l < QK8/4; l++) { - minv[4*l] = _mm256_min_ps(minv[4*l], minv[4*l+2]); - maxv[4*l] = _mm256_max_ps(maxv[4*l], maxv[4*l+2]); - } - - for (int l = 0; l < QK8/8; l++) { - minv[8*l] = _mm256_min_ps(minv[8*l], minv[8*l+4]); - maxv[8*l] = _mm256_max_ps(maxv[8*l], maxv[8*l+4]); - } - - //min = MIN(minv[0][0], MIN(minv[0][1], MIN(minv[0][2], MIN(minv[0][3], MIN(minv[0][4], MIN(minv[0][5], MIN(minv[0][6], minv[0][7]))))))); - //max = MAX(maxv[0][0], MAX(maxv[0][1], MAX(maxv[0][2], MAX(maxv[0][3], MAX(maxv[0][4], MAX(maxv[0][5], MAX(maxv[0][6], maxv[0][7]))))))); - - const __m256 minv0_0 = _mm256_permute2f128_ps(minv[0], minv[0], 3); - const __m256 minv0_1 = _mm256_min_ps(minv[0], minv0_0); - const __m256 minv0_2 = _mm256_permute_ps(minv0_1, 0x4e); - const __m256 minv0_3 = _mm256_min_ps(minv0_1, minv0_2); - const __m256 minv0_4 = _mm256_permute_ps(minv0_3, 0xb1); - const __m256 minv0_5 = _mm256_min_ps(minv0_3, minv0_4); - - const __m256 maxv0_0 = _mm256_permute2f128_ps(maxv[0], maxv[0], 3); - const __m256 maxv0_1 = _mm256_max_ps(maxv[0], maxv0_0); - const __m256 maxv0_2 = _mm256_permute_ps(maxv0_1, 0x4e); - const __m256 maxv0_3 = _mm256_max_ps(maxv0_1, maxv0_2); - const __m256 maxv0_4 = _mm256_permute_ps(maxv0_3, 0xb1); - const __m256 maxv0_5 = _mm256_max_ps(maxv0_3, maxv0_4); - - min = _mm256_cvtss_f32(minv0_5); - max = _mm256_cvtss_f32(maxv0_5); - - const float d = (max - min) / ((1 << QB) - 2); - const float id = d ? 1.0/d : 0.0; - - pm[i] = GGML_FP32_TO_GQ(min); - pd[i] = GGML_FP32_TO_GQ(d); - - const __m256 idv = _mm256_set1_ps(id); - - for (int l = 0; l < QK/8; l++) { - __m256 v = _mm256_mul_ps(_mm256_sub_ps(srcv[l], _mm256_set1_ps(min)), idv); -#if 0 - v[0] += frand(); v[1] += frand(); v[2] += frand(); v[3] += frand(); - v[4] += frand(); v[5] += frand(); v[6] += frand(); v[7] += frand(); -#endif - - // convert to uint8 - __m256i vi = _mm256_cvtps_epi32(v); - - uint32_t vi_0 = _mm256_extract_epi32(vi, 0); - uint32_t vi_1 = _mm256_extract_epi32(vi, 1); - uint32_t vi_2 = _mm256_extract_epi32(vi, 2); - uint32_t vi_3 = _mm256_extract_epi32(vi, 3); - - uint32_t vi_4 = _mm256_extract_epi32(vi, 4); - uint32_t vi_5 = _mm256_extract_epi32(vi, 5); - uint32_t vi_6 = _mm256_extract_epi32(vi, 6); - uint32_t vi_7 = _mm256_extract_epi32(vi, 7); - - // convert to 4-bit, 2 consecutive packed into 1 byte - pp[4*l + 0] = vi_0 | (vi_1 << 4); - pp[4*l + 1] = vi_2 | (vi_3 << 4); - pp[4*l + 2] = vi_4 | (vi_5 << 4); - pp[4*l + 3] = vi_6 | (vi_7 << 4); - - //printf("vi: %7d %7d %7d %7d %7d %7d %7d %7d\n", vi_0, vi_1, vi_2, vi_3, vi_4, vi_5, vi_6, vi_7); - //printf("v : %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]); - } - - memcpy(pb + i*QK/2, pp, sizeof(pp)); - } -#elif defined(__ARM_NEON) && 0 - { - // TODO - } -#else - { - for (int l = 0; l < QK; l++) { - const float v = src[i*QK + l]; - if (v < min) min = v; - if (v > max) max = v; - } - - const float d = (max - min) / ((1 << QB) - 1); - const float id = d ? 1.0/d : 0.0; - - pm[i] = GGML_FP32_TO_GQ(min); - pd[i] = GGML_FP32_TO_GQ(d); - - for (int l = 0; l < QK; l++) { - const float v = (src[i*QK + l] - min) * id; - const uint8_t vi = (uint8_t) (v + frand()); - pp[l/2] |= (vi & 0xf) << (4*(l & 1)); - } - - memcpy(pb + i*QK/2, pp, sizeof(pp)); - } -#endif - //printf("min %f max %f\n", min, max); - } -} - -// reimplementation of quantize_4 using quantize_4_row -void quantize_4(const float * restrict src, char * restrict dst, int n, int k) { - assert(k % QK == 0); - - for (int j = 0; j < n; j++) { - quantize_4_row(src + j*k, dst, k); - dst = (char *) dst + quantize_4_row_size(k); - } -} - -void vec_dot_gq_4(const int n, float * restrict s, const void * restrict x, const void * restrict y) { - const int nb = quantize_4_blocks_per_row(n); - - const gq_scale_t * restrict pm0 = (const gq_scale_t *) x; - const gq_scale_t * restrict pm1 = (const gq_scale_t *) y; - - const gq_scale_t * restrict pd0 = pm0 + nb; - const gq_scale_t * restrict pd1 = pm1 + nb; - - const uint8_t * restrict pb0 = (const uint8_t *) (pd0 + nb); - const uint8_t * restrict pb1 = (const uint8_t *) (pd1 + nb); - - float sumf = 0.0; - -#if 0 - // scalar - for (int i = 0; i < nb; i++) { - const float m0 = GGML_GQ_TO_FP32(pm0[i]); - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - - const float m1 = GGML_GQ_TO_FP32(pm1[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - for (int j = 0; j < QK/2; j++) { - const uint8_t v0 = p0[j]; - const uint8_t v1 = p1[j]; - - const float f0 = d0*(v0 & 0xf) + m0; - const float f1 = d0*(v0 >> 4) + m0; - - const float f2 = d1*(v1 & 0xf) + m1; - const float f3 = d1*(v1 >> 4) + m1; - - sumf += f0*f2 + f1*f3; - } - } -#else -#if defined(__AVX2__) -#if QK == 64 && 0 - __m256 sumv0 = _mm256_setzero_ps(); - __m256 sumv1 = _mm256_setzero_ps(); - - for (int i = 0; i < nb; i++) { - const float m0 = GGML_GQ_TO_FP32(pm0[i]); - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - - const float m1 = GGML_GQ_TO_FP32(pm1[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - const __m256 m0v = _mm256_set1_ps(m0); - const __m256 d0v = _mm256_set1_ps(d0); - - const __m256 m1v = _mm256_set1_ps(m1); - const __m256 d1v = _mm256_set1_ps(d1); - - const __m256i m4b = _mm256_set1_epi8(0xf); - - __m256i v0 = _mm256_loadu_si256((__m256i *) p0); - - //_mm_prefetch((const char *) (p0 + 32), _MM_HINT_T0); - //_mm_prefetch((const char *) (p1 + 32), _MM_HINT_T0); - //_mm_prefetch((const char *) (pm0 + i + 1), _MM_HINT_T0); - //_mm_prefetch((const char *) (pm1 + i + 1), _MM_HINT_T0); - //_mm_prefetch((const char *) (pd0 + i + 1), _MM_HINT_T0); - //_mm_prefetch((const char *) (pd1 + i + 1), _MM_HINT_T0); - - __m256i v00 = _mm256_and_si256(v0, _mm256_set1_epi32(0x000000FF)); - __m256i v01 = _mm256_srli_epi32(_mm256_and_si256(v0, _mm256_set1_epi32(0x0000FFFF)), 8); - __m256i v02 = _mm256_srli_epi32(_mm256_and_si256(v0, _mm256_set1_epi32(0x00FFFFFF)), 16); - __m256i v03 = _mm256_srli_epi32(v0, 24); - - ////////////////////// - - //{ - // uint32_t vi_0 = _mm256_extract_epi32(v00, 0); - // uint32_t vi_1 = _mm256_extract_epi32(v00, 1); - // uint32_t vi_2 = _mm256_extract_epi32(v00, 2); - // uint32_t vi_3 = _mm256_extract_epi32(v00, 3); - // uint32_t vi_4 = _mm256_extract_epi32(v00, 4); - // uint32_t vi_5 = _mm256_extract_epi32(v00, 5); - // uint32_t vi_6 = _mm256_extract_epi32(v00, 6); - // uint32_t vi_7 = _mm256_extract_epi32(v00, 7); - // printf("v0: %7d %7d %7d %7d %7d %7d %7d %7d\n", vi_0, vi_1, vi_2, vi_3, vi_4, vi_5, vi_6, vi_7); - // printf("p0: %7d %7d %7d %7d %7d %7d %7d %7d\n", p0[0], p0[4], p0[8], p0[12], p0[16], p0[20], p0[24], p0[28]); - // printf("p1: %7d %7d %7d %7d %7d %7d %7d %7d\n", p0[1], p0[5], p0[9], p0[13], p0[17], p0[21], p0[25], p0[29]); - // printf("p2: %7d %7d %7d %7d %7d %7d %7d %7d\n", p0[2], p0[6], p0[10], p0[14], p0[18], p0[22], p0[26], p0[30]); - // printf("p3: %7d %7d %7d %7d %7d %7d %7d %7d\n", p0[3], p0[7], p0[11], p0[15], p0[19], p0[23], p0[27], p0[31]); - //} - - // compute 32 x 4-bit values (low and high) - __m256i v00l = _mm256_and_si256(v00, m4b); - __m256i v01l = _mm256_and_si256(v01, m4b); - __m256i v02l = _mm256_and_si256(v02, m4b); - __m256i v03l = _mm256_and_si256(v03, m4b); - - __m256i v00h = _mm256_srli_epi32(v00, 4); - __m256i v01h = _mm256_srli_epi32(v01, 4); - __m256i v02h = _mm256_srli_epi32(v02, 4); - __m256i v03h = _mm256_srli_epi32(v03, 4); - - //{ - // uint32_t vi_0 = _mm256_extract_epi32(v00l, 0); - // uint32_t vi_1 = _mm256_extract_epi32(v00l, 1); - // uint32_t vi_2 = _mm256_extract_epi32(v00l, 2); - // uint32_t vi_3 = _mm256_extract_epi32(v00l, 3); - // uint32_t vi_4 = _mm256_extract_epi32(v00l, 4); - // uint32_t vi_5 = _mm256_extract_epi32(v00l, 5); - // uint32_t vi_6 = _mm256_extract_epi32(v00l, 6); - // uint32_t vi_7 = _mm256_extract_epi32(v00l, 7); - - // printf("v0l: %7d %7d %7d %7d %7d %7d %7d %7d\n", vi_0, vi_1, vi_2, vi_3, vi_4, vi_5, vi_6, vi_7); - - // vi_0 = _mm256_extract_epi32(v00h, 0); - // vi_1 = _mm256_extract_epi32(v00h, 1); - // vi_2 = _mm256_extract_epi32(v00h, 2); - // vi_3 = _mm256_extract_epi32(v00h, 3); - // vi_4 = _mm256_extract_epi32(v00h, 4); - // vi_5 = _mm256_extract_epi32(v00h, 5); - // vi_6 = _mm256_extract_epi32(v00h, 6); - // vi_7 = _mm256_extract_epi32(v00h, 7); - - // printf("v0h: %7d %7d %7d %7d %7d %7d %7d %7d\n", vi_0, vi_1, vi_2, vi_3, vi_4, vi_5, vi_6, vi_7); - //} - - // convert to float - __m256 vf00l = _mm256_cvtepi32_ps(v00l); - __m256 vf01l = _mm256_cvtepi32_ps(v01l); - __m256 vf02l = _mm256_cvtepi32_ps(v02l); - __m256 vf03l = _mm256_cvtepi32_ps(v03l); - - __m256 vf00h = _mm256_cvtepi32_ps(v00h); - __m256 vf01h = _mm256_cvtepi32_ps(v01h); - __m256 vf02h = _mm256_cvtepi32_ps(v02h); - __m256 vf03h = _mm256_cvtepi32_ps(v03h); - - //{ - // printf("vf00l: %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", vf00l[0], vf00l[1], vf00l[2], vf00l[3], vf00l[4], vf00l[5], vf00l[6], vf00l[7]); - // printf("vf01l: %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", vf01l[0], vf01l[1], vf01l[2], vf01l[3], vf01l[4], vf01l[5], vf01l[6], vf01l[7]); - // printf("vf02l: %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", vf02l[0], vf02l[1], vf02l[2], vf02l[3], vf02l[4], vf02l[5], vf02l[6], vf02l[7]); - // printf("vf03l: %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", vf03l[0], vf03l[1], vf03l[2], vf03l[3], vf03l[4], vf03l[5], vf03l[6], vf03l[7]); - //} - - // multiply by scale and add offset - vf00l = _mm256_fmadd_ps(vf00l, d0v, m0v); - vf01l = _mm256_fmadd_ps(vf01l, d0v, m0v); - vf02l = _mm256_fmadd_ps(vf02l, d0v, m0v); - vf03l = _mm256_fmadd_ps(vf03l, d0v, m0v); - - vf00h = _mm256_fmadd_ps(vf00h, d0v, m0v); - vf01h = _mm256_fmadd_ps(vf01h, d0v, m0v); - vf02h = _mm256_fmadd_ps(vf02h, d0v, m0v); - vf03h = _mm256_fmadd_ps(vf03h, d0v, m0v); - - __m256i v1 = _mm256_loadu_si256((__m256i *) p1); - - __m256i v10 = _mm256_and_si256(v1, _mm256_set1_epi32(0x000000FF)); - __m256i v11 = _mm256_srli_epi32(_mm256_and_si256(v1, _mm256_set1_epi32(0x0000FFFF)), 8); - __m256i v12 = _mm256_srli_epi32(_mm256_and_si256(v1, _mm256_set1_epi32(0x00FFFFFF)), 16); - __m256i v13 = _mm256_srli_epi32(v1, 24); - - __m256i v10l = _mm256_and_si256(v10, m4b); - __m256i v11l = _mm256_and_si256(v11, m4b); - __m256i v12l = _mm256_and_si256(v12, m4b); - __m256i v13l = _mm256_and_si256(v13, m4b); - - __m256i v10h = _mm256_srli_epi32(v10, 4); - __m256i v11h = _mm256_srli_epi32(v11, 4); - __m256i v12h = _mm256_srli_epi32(v12, 4); - __m256i v13h = _mm256_srli_epi32(v13, 4); - - __m256 vf10l = _mm256_cvtepi32_ps(v10l); - __m256 vf11l = _mm256_cvtepi32_ps(v11l); - __m256 vf12l = _mm256_cvtepi32_ps(v12l); - __m256 vf13l = _mm256_cvtepi32_ps(v13l); - - __m256 vf10h = _mm256_cvtepi32_ps(v10h); - __m256 vf11h = _mm256_cvtepi32_ps(v11h); - __m256 vf12h = _mm256_cvtepi32_ps(v12h); - __m256 vf13h = _mm256_cvtepi32_ps(v13h); - - vf10l = _mm256_fmadd_ps(vf10l, d1v, m1v); - vf11l = _mm256_fmadd_ps(vf11l, d1v, m1v); - vf12l = _mm256_fmadd_ps(vf12l, d1v, m1v); - vf13l = _mm256_fmadd_ps(vf13l, d1v, m1v); - - vf10h = _mm256_fmadd_ps(vf10h, d1v, m1v); - vf11h = _mm256_fmadd_ps(vf11h, d1v, m1v); - vf12h = _mm256_fmadd_ps(vf12h, d1v, m1v); - vf13h = _mm256_fmadd_ps(vf13h, d1v, m1v); - - // compute dot product - sumv0 = _mm256_fmadd_ps(vf00l, vf10l, sumv0); - sumv0 = _mm256_fmadd_ps(vf01l, vf11l, sumv0); - sumv0 = _mm256_fmadd_ps(vf02l, vf12l, sumv0); - sumv0 = _mm256_fmadd_ps(vf03l, vf13l, sumv0); - - sumv1 = _mm256_fmadd_ps(vf00h, vf10h, sumv1); - sumv1 = _mm256_fmadd_ps(vf01h, vf11h, sumv1); - sumv1 = _mm256_fmadd_ps(vf02h, vf12h, sumv1); - sumv1 = _mm256_fmadd_ps(vf03h, vf13h, sumv1); - } - - // accumulate (horizontal sum) - const __m256 vdot = _mm256_add_ps(sumv0, sumv1); - const __m128 t0 = _mm_add_ps(_mm256_castps256_ps128(vdot), _mm256_extractf128_ps(vdot, 1)); - const __m128 t1 = _mm_hadd_ps(t0, t0); - - sumf += _mm_cvtss_f32(_mm_hadd_ps(t1, t1)); -#elif QK == 64 && 0 - float sum00 = 0.0f; - float sum01 = 0.0f; - float sum10 = 0.0f; - float sum11 = 0.0f; - - const __m256i m4b = _mm256_set1_epi8(0xf); - - for (int i = 0; i < nb; i++) { - const float m0 = GGML_GQ_TO_FP32(pm0[i]); - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - - const float m1 = GGML_GQ_TO_FP32(pm1[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - // 64 x 4 - const __m256i v0 = _mm256_loadu_si256((__m256i *) p0); - const __m256i v1 = _mm256_loadu_si256((__m256i *) p1); - - // 32 x 8 - const __m256i v0l = _mm256_and_si256(v0, m4b); - const __m256i v1l = _mm256_and_si256(v1, m4b); - - const __m256i v0h = _mm256_and_si256(_mm256_srli_epi16(v0, 4), m4b); - const __m256i v1h = _mm256_and_si256(_mm256_srli_epi16(v1, 4), m4b); - - const __m256i pl = _mm256_maddubs_epi16(v0l, v1l); - const __m256i ph = _mm256_maddubs_epi16(v0h, v1h); - - const __m256i p16 = _mm256_add_epi16(ph, pl); - const __m256i p = _mm256_madd_epi16(_mm256_set1_epi16(1), p16); - - sum00 += m0*m1; - sum01 += m1*d0*(_mm256_hadd_epi8_gg(_mm256_add_epi8(v0l, v0h))); - sum10 += m0*d1*(_mm256_hadd_epi8_gg(_mm256_add_epi8(v1l, v1h))); - sum11 += d0*d1*(_mm256_hadd_epi32_gg(p)); - } - - sumf = 64.0*sum00 + sum01 + sum10 + sum11; -#elif QK == 64 && 1 // this is the best when using min + d - float sum00 = 0.0f; - - __m256 sum01 = _mm256_setzero_ps(); - __m256 sum10 = _mm256_setzero_ps(); - __m256 sum11 = _mm256_setzero_ps(); - - for (int i = 0; i < nb; i++) { - const float m0 = GGML_GQ_TO_FP32(pm0[i]); - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - - const float m1 = GGML_GQ_TO_FP32(pm1[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - const __m256 m0v = _mm256_set1_ps(m0); - const __m256 d0v = _mm256_set1_ps(d0); - - const __m256 m1v = _mm256_set1_ps(m1); - const __m256 d1v = _mm256_set1_ps(d1); - - const __m256 m1d0v = _mm256_mul_ps(m1v, d0v); - const __m256 m0d1v = _mm256_mul_ps(m0v, d1v); - const __m256 d0d1v = _mm256_mul_ps(d0v, d1v); - - const __m256i m4b = _mm256_set1_epi8(0xf); - - // 64 x 4 - const __m256i v0 = _mm256_loadu_si256((__m256i *) p0); - const __m256i v1 = _mm256_loadu_si256((__m256i *) p1); - - // 32 x 8 - const __m256i v0l = _mm256_and_si256(v0, m4b); - const __m256i v1l = _mm256_and_si256(v1, m4b); - - const __m256i v0h = _mm256_and_si256(_mm256_srli_epi16(v0, 4), m4b); - const __m256i v1h = _mm256_and_si256(_mm256_srli_epi16(v1, 4), m4b); - - const __m256i v0a = _mm256_add_epi8(v0l, v0h); - const __m256i v1a = _mm256_add_epi8(v1l, v1h); - - const __m128i v0al = _mm256_extracti128_si256(v0a, 0); - const __m128i v0ah = _mm256_extracti128_si256(v0a, 1); - - const __m128i v1al = _mm256_extracti128_si256(v1a, 0); - const __m128i v1ah = _mm256_extracti128_si256(v1a, 1); - - const __m128i v0as = _mm_add_epi8(v0al, v0ah); - const __m128i v1as = _mm_add_epi8(v1al, v1ah); - - const __m256i v0as_0 = _mm256_cvtepu8_epi32(v0as); - const __m256i v0as_1 = _mm256_cvtepu8_epi32(_mm_srli_si128(v0as, 8)); - - const __m256i v1as_0 = _mm256_cvtepu8_epi32(v1as); - const __m256i v1as_1 = _mm256_cvtepu8_epi32(_mm_srli_si128(v1as, 8)); - - const __m256i v0ass = _mm256_add_epi32(v0as_0, v0as_1); - const __m256i v1ass = _mm256_add_epi32(v1as_0, v1as_1); - - const __m256 v0f = _mm256_cvtepi32_ps(v0ass); - const __m256 v1f = _mm256_cvtepi32_ps(v1ass); - - const __m256i pl = _mm256_maddubs_epi16(v0l, v1l); - const __m256i ph = _mm256_maddubs_epi16(v0h, v1h); - - const __m256i p16 = _mm256_add_epi16(ph, pl); - const __m256i p = _mm256_madd_epi16(_mm256_set1_epi16(1), p16); - - sum00 += m0*m1; - sum01 = _mm256_fmadd_ps(m1d0v, v0f, sum01); - sum10 = _mm256_fmadd_ps(m0d1v, v1f, sum10); - sum11 = _mm256_fmadd_ps(d0d1v, _mm256_cvtepi32_ps(p), sum11); - } - - sumf = 64.0*sum00 + _mm256_hadd_ps_gg(sum01) + _mm256_hadd_ps_gg(sum10) + _mm256_hadd_ps_gg(sum11); -#endif -#elif defined (__ARM_NEON) - float sum00 = 0.0f; - float sum01 = 0.0f; - float sum10 = 0.0f; - float sum11 = 0.0f; - - for (int i = 0; i < nb; i++) { - const float m0 = GGML_GQ_TO_FP32(pm0[i]); - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - - const float m1 = GGML_GQ_TO_FP32(pm1[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - const uint8x16_t m4b = vdupq_n_u8(0xf); - - const uint8x16_t v0_0 = vld1q_u8(p0); - const uint8x16_t v0_1 = vld1q_u8(p0 + 16); - const uint8x16_t v1_0 = vld1q_u8(p1); - const uint8x16_t v1_1 = vld1q_u8(p1 + 16); - - // and with 0xf - const uint8x16_t v0_0l = vandq_u8(v0_0, m4b); - const uint8x16_t v0_1l = vandq_u8(v0_1, m4b); - const uint8x16_t v1_0l = vandq_u8(v1_0, m4b); - const uint8x16_t v1_1l = vandq_u8(v1_1, m4b); - - const uint8x16_t v0_0h = vshrq_n_u8(v0_0, 4); - const uint8x16_t v0_1h = vshrq_n_u8(v0_1, 4); - const uint8x16_t v1_0h = vshrq_n_u8(v1_0, 4); - const uint8x16_t v1_1h = vshrq_n_u8(v1_1, 4); - - // dot product into uint16x8_t - const uint16x8_t pl0l = vmull_u8(vget_low_u8 (v0_0l), vget_low_u8 (v1_0l)); - const uint16x8_t pl0h = vmull_u8(vget_high_u8(v0_0l), vget_high_u8(v1_0l)); - const uint16x8_t pl1l = vmull_u8(vget_low_u8 (v0_1l), vget_low_u8 (v1_1l)); - const uint16x8_t pl1h = vmull_u8(vget_high_u8(v0_1l), vget_high_u8(v1_1l)); - - const uint16x8_t ph0l = vmull_u8(vget_low_u8 (v0_0h), vget_low_u8 (v1_0h)); - const uint16x8_t ph0h = vmull_u8(vget_high_u8(v0_0h), vget_high_u8(v1_0h)); - const uint16x8_t ph1l = vmull_u8(vget_low_u8 (v0_1h), vget_low_u8 (v1_1h)); - const uint16x8_t ph1h = vmull_u8(vget_high_u8(v0_1h), vget_high_u8(v1_1h)); - - const uint16x8_t pl0 = vaddq_u16(pl0l, pl0h); - const uint16x8_t pl1 = vaddq_u16(pl1l, pl1h); - const uint16x8_t ph0 = vaddq_u16(ph0l, ph0h); - const uint16x8_t ph1 = vaddq_u16(ph1l, ph1h); - - const uint16x8_t pl = vaddq_u16(pl0, pl1); - const uint16x8_t ph = vaddq_u16(ph0, ph1); - - sum00 += m0*m1; - sum01 += m1*d0*(vaddvq_u8(v0_0l) + vaddvq_u8(v0_0h) + vaddvq_u8(v0_1l) + vaddvq_u8(v0_1h)); - sum10 += m0*d1*(vaddvq_u8(v1_0l) + vaddvq_u8(v1_0h) + vaddvq_u8(v1_1l) + vaddvq_u8(v1_1h)); - //sum11 += d0*d1*( - // vaddvq_u16(vaddq_u16(vaddq_u16(pl0l, pl0h), vaddq_u16(pl1l, pl1h))) + - // vaddvq_u16(vaddq_u16(vaddq_u16(ph0l, ph0h), vaddq_u16(ph1l, ph1h)))); - sum11 += d0*d1*vaddvq_u16(vaddq_u16(pl, ph)); - } - - sumf = 64.0*sum00 + sum01 + sum10 + sum11; -#endif -#endif - - *s = sumf; -} - -// use vec_dot_gq_4 to compute the dot product of two rows -void mul_mat_gq_4( - const void * src0, - const void * src1, // transposed - float * dst, - int m, int n, int k) { - assert(k % QK == 0); - - const int nb = quantize_4_blocks_per_row(k); - - for (int ir0 = 0; ir0 < m; ir0++) { - for (int ir1 = 0; ir1 < n; ir1++) { - vec_dot_gq_4(k, dst + ir1, src0, src1); - src1 = (const char *) src1 + quantize_4_row_size(k); - } - src0 = (const char *) src0 + quantize_4_row_size(k); - src1 = (const char *) src1 - n*quantize_4_row_size(k); - - dst = (float *) dst + n; - } -} - -// -// method 5 -// 4-bit quantization (without min, only delta) -// - -static inline int quantize_5_blocks_per_row(int k) { - return k/QK; -} - -static inline int quantize_5_row_size(int k) { - const int nb = quantize_5_blocks_per_row(k); - - return nb*(sizeof(gq_scale_t) + QK/2); -} - -void quantize_5_row(const float * restrict src, void * restrict dst, int k) { - assert(k % QK == 0); - assert(QB == 4); - - const int nb = quantize_5_blocks_per_row(k); - - gq_scale_t * restrict pd = (gq_scale_t *) (dst); - uint8_t * restrict pb = (uint8_t *) (pd + nb); - - uint8_t pp[QK/2]; - - for (int i = 0; i < nb; i++) { - memset(pp, 0, sizeof(pp)); - - float amax = 0.0f; // absolute max - -#if defined(__AVX2__) - { - assert(QK == 64); - enum { QK8 = QK/8 }; - - __m256 srcv [QK8]; - __m256 asrcv[QK8]; - __m256 amaxv[QK8]; - - for (int l = 0; l < QK8; l++) { - srcv[l] = _mm256_loadu_ps(src + i*QK + 8*l); - } - - for (int l = 0; l < QK8; l++) { - asrcv[l] = _mm256_and_ps(srcv[l], _mm256_castsi256_ps(_mm256_set1_epi32(0x7fffffff))); - } - - - for (int l = 0; l < QK8/2; l++) { - amaxv[2*l] = _mm256_max_ps(asrcv[2*l], asrcv[2*l+1]); - } - - for (int l = 0; l < QK8/4; l++) { - amaxv[4*l] = _mm256_max_ps(amaxv[4*l], amaxv[4*l+2]); - } - - for (int l = 0; l < QK8/8; l++) { - amaxv[8*l] = _mm256_max_ps(amaxv[8*l], amaxv[8*l+4]); - } - - //amax = MAX(amaxv[0][0], MAX(amaxv[0][1], MAX(amaxv[0][2], MAX(amaxv[0][3], MAX(amaxv[0][4], MAX(amaxv[0][5], MAX(amaxv[0][6], amaxv[0][7]))))))); - - const __m256 amaxv0_0 = _mm256_permute2f128_ps(amaxv[0], amaxv[0], 3); - const __m256 amaxv0_1 = _mm256_max_ps(amaxv[0], amaxv0_0); - const __m256 amaxv0_2 = _mm256_permute_ps(amaxv0_1, 0x4e); - const __m256 amaxv0_3 = _mm256_max_ps(amaxv0_1, amaxv0_2); - const __m256 amaxv0_4 = _mm256_permute_ps(amaxv0_3, 0xb1); - const __m256 amaxv0_5 = _mm256_max_ps(amaxv0_3, amaxv0_4); - - amax = _mm256_cvtss_f32(amaxv0_5); - - //printf("amax = %f\n", amax); - - const float d = amax / ((1 << (QB - 1)) - 1); - const float id = d ? 1.0/d : 0.0; - - pd[i] = GGML_FP32_TO_GQ(d); - - const __m256 idv = _mm256_set1_ps(id); - - for (int l = 0; l < QK/8; l++) { - __m256 v = _mm256_mul_ps(srcv[l], idv); -#if 0 - v[0] += frand(); v[1] += frand(); v[2] += frand(); v[3] += frand(); - v[4] += frand(); v[5] += frand(); v[6] += frand(); v[7] += frand(); -#endif - - // convert to int8 - __m256i vi = _mm256_cvtps_epi32(v); - vi = _mm256_add_epi32(vi, _mm256_set1_epi32(8)); - - int32_t vi_0 = _mm256_extract_epi32(vi, 0); - int32_t vi_1 = _mm256_extract_epi32(vi, 1); - int32_t vi_2 = _mm256_extract_epi32(vi, 2); - int32_t vi_3 = _mm256_extract_epi32(vi, 3); - - int32_t vi_4 = _mm256_extract_epi32(vi, 4); - int32_t vi_5 = _mm256_extract_epi32(vi, 5); - int32_t vi_6 = _mm256_extract_epi32(vi, 6); - int32_t vi_7 = _mm256_extract_epi32(vi, 7); - - // convert to 4-bit, 2 consecutive packed into 1 byte - pp[4*l + 0] = vi_0 | (vi_1 << 4); - pp[4*l + 1] = vi_2 | (vi_3 << 4); - pp[4*l + 2] = vi_4 | (vi_5 << 4); - pp[4*l + 3] = vi_6 | (vi_7 << 4); - - //printf("vi: %7d %7d %7d %7d %7d %7d %7d %7d\n", vi_0, vi_1, vi_2, vi_3, vi_4, vi_5, vi_6, vi_7); - ////printf("v : %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]); - - assert(vi_0 >= 0 && vi_0 < 16); - assert(vi_1 >= 0 && vi_1 < 16); - assert(vi_2 >= 0 && vi_2 < 16); - assert(vi_3 >= 0 && vi_3 < 16); - - assert(vi_4 >= 0 && vi_4 < 16); - assert(vi_5 >= 0 && vi_5 < 16); - assert(vi_6 >= 0 && vi_6 < 16); - assert(vi_7 >= 0 && vi_7 < 16); - } - - memcpy(pb + i*QK/2, pp, sizeof(pp)); - } -#elif defined(__ARM_NEON) && 0 - { - // TODO - } -#else - { - for (int l = 0; l < QK; l++) { - const float v = src[i*QK + l]; - amax = MAX(amax, fabsf(v)); - } - - const float d = amax / ((1 << (QB - 1)) - 1); - const float id = d ? 1.0/d : 0.0; - - pd[i] = GGML_FP32_TO_GQ(d); - - for (int l = 0; l < QK; l++) { - const float v = src[i*QK + l]*id; - const int8_t vi = ((int8_t) (round(v))) + 8; - assert(vi >= 0 && vi < 16); - pp[l/2] |= (vi & 0xf) << (4*(l & 1)); - } - - memcpy(pb + i*QK/2, pp, sizeof(pp)); - } -#endif - //printf("min %f max %f\n", min, max); - } -} - -// reimplementation of quantize_5 using quantize_5_row -void quantize_5(const float * restrict src, char * restrict dst, int n, int k) { - assert(k % QK == 0); - - for (int j = 0; j < n; j++) { - quantize_5_row(src + j*k, dst, k); - dst = (char *) dst + quantize_5_row_size(k); - } -} - -void vec_dot_gq_5(const int n, float * restrict s, const void * restrict x, const void * restrict y) { - const int nb = quantize_5_blocks_per_row(n); - - const gq_scale_t * restrict pd0 = (const gq_scale_t *) x; - const gq_scale_t * restrict pd1 = (const gq_scale_t *) y; - - const uint8_t * restrict pb0 = (const uint8_t *) (pd0 + nb); - const uint8_t * restrict pb1 = (const uint8_t *) (pd1 + nb); - - float sumf = 0.0; - -#if 0 - // scalar - for (int i = 0; i < nb; i++) { - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - for (int j = 0; j < QK/2; j++) { - const uint8_t v0 = p0[j]; - const uint8_t v1 = p1[j]; - - const float f0 = d0*((int8_t) (v0 & 0xf) - 8); - const float f1 = d0*((int8_t) (v0 >> 4) - 8); - - const float f2 = d1*((int8_t) (v1 & 0xf) - 8); - const float f3 = d1*((int8_t) (v1 >> 4) - 8); - - sumf += f0*f2 + f1*f3; - } - } -#else -#if defined(__AVX2__) -#if QK == 64 && 1 - __m256 sum11 = _mm256_setzero_ps(); - - for (int i = 0; i < nb; i++) { - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - const __m256 d0v = _mm256_set1_ps(d0); - const __m256 d1v = _mm256_set1_ps(d1); - - const __m256 d0d1v = _mm256_mul_ps(d0v, d1v); - - const __m256i m4b = _mm256_set1_epi8(0xf); - - // 64 x 4 - const __m256i v0 = _mm256_loadu_si256((__m256i *) p0); - const __m256i v1 = _mm256_loadu_si256((__m256i *) p1); - - // 32 x 8 - __m256i v0l = _mm256_and_si256(v0, m4b); - __m256i v1l = _mm256_and_si256(v1, m4b); - - __m256i v0h = _mm256_and_si256(_mm256_srli_epi16(v0, 4), m4b); - __m256i v1h = _mm256_and_si256(_mm256_srli_epi16(v1, 4), m4b); - - // sub 8 - v0l = _mm256_sub_epi8(v0l, _mm256_set1_epi8(8)); - v0h = _mm256_sub_epi8(v0h, _mm256_set1_epi8(8)); - - v1l = _mm256_sub_epi8(v1l, _mm256_set1_epi8(8)); - v1h = _mm256_sub_epi8(v1h, _mm256_set1_epi8(8)); - - // abs - const __m256i v0la = _mm256_sign_epi8(v0l, v0l); - const __m256i v0ha = _mm256_sign_epi8(v0h, v0h); - - // sign - const __m256i v1ls = _mm256_sign_epi8(v1l, v0l); - const __m256i v1hs = _mm256_sign_epi8(v1h, v0h); - - const __m256i pl = _mm256_maddubs_epi16(v0la, v1ls); - const __m256i ph = _mm256_maddubs_epi16(v0ha, v1hs); - - const __m256i p16 = _mm256_add_epi16(ph, pl); - const __m256i p = _mm256_madd_epi16(_mm256_set1_epi16(1), p16); - - sum11 = _mm256_fmadd_ps(d0d1v, _mm256_cvtepi32_ps(p), sum11); - } - - sumf = _mm256_hadd_ps_gg(sum11); -#endif -#elif defined (__ARM_NEON) - float sum11 = 0.0f; - - //float32x4_t sum_0 = vdupq_n_f32(0.0f); - //float32x4_t sum_1 = vdupq_n_f32(0.0f); - - //float16x8_t sum_0 = vdupq_n_f16(0.0f); - //float16x8_t sum_1 = vdupq_n_f16(0.0f); - - for (int i = 0; i < nb; i++) { - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - //float32x4_t d0d1v = vdupq_n_f32(d0*d1); - //float16x8_t d0d1v = vdupq_n_f16(d0*d1); - - const uint8_t * restrict p0 = pb0 + i*QK/2; - const uint8_t * restrict p1 = pb1 + i*QK/2; - - const uint8x16_t m4b = vdupq_n_u8(0xf); - const int8x16_t s8b = vdupq_n_s8(0x8); - - const uint8x16_t v0_0 = vld1q_u8(p0); - const uint8x16_t v0_1 = vld1q_u8(p0 + 16); - const uint8x16_t v1_0 = vld1q_u8(p1); - const uint8x16_t v1_1 = vld1q_u8(p1 + 16); - - // 4-bit -> 8-bit - const int8x16_t v0_0l = vreinterpretq_s8_u8(vandq_u8(v0_0, m4b)); - const int8x16_t v0_1l = vreinterpretq_s8_u8(vandq_u8(v0_1, m4b)); - const int8x16_t v1_0l = vreinterpretq_s8_u8(vandq_u8(v1_0, m4b)); - const int8x16_t v1_1l = vreinterpretq_s8_u8(vandq_u8(v1_1, m4b)); - - const int8x16_t v0_0h = vreinterpretq_s8_u8(vshrq_n_u8(v0_0, 4)); - const int8x16_t v0_1h = vreinterpretq_s8_u8(vshrq_n_u8(v0_1, 4)); - const int8x16_t v1_0h = vreinterpretq_s8_u8(vshrq_n_u8(v1_0, 4)); - const int8x16_t v1_1h = vreinterpretq_s8_u8(vshrq_n_u8(v1_1, 4)); - - // sub 8 - const int8x16_t v0_0ls = vsubq_s8(v0_0l, s8b); - const int8x16_t v0_1ls = vsubq_s8(v0_1l, s8b); - const int8x16_t v1_0ls = vsubq_s8(v1_0l, s8b); - const int8x16_t v1_1ls = vsubq_s8(v1_1l, s8b); - - const int8x16_t v0_0hs = vsubq_s8(v0_0h, s8b); - const int8x16_t v0_1hs = vsubq_s8(v0_1h, s8b); - const int8x16_t v1_0hs = vsubq_s8(v1_0h, s8b); - const int8x16_t v1_1hs = vsubq_s8(v1_1h, s8b); - - // dot product into int16x8_t - const int16x8_t pl0l = vmull_s8(vget_low_s8 (v0_0ls), vget_low_s8 (v1_0ls)); - const int16x8_t pl0h = vmull_s8(vget_high_s8(v0_0ls), vget_high_s8(v1_0ls)); - const int16x8_t pl1l = vmull_s8(vget_low_s8 (v0_1ls), vget_low_s8 (v1_1ls)); - const int16x8_t pl1h = vmull_s8(vget_high_s8(v0_1ls), vget_high_s8(v1_1ls)); - - const int16x8_t ph0l = vmull_s8(vget_low_s8 (v0_0hs), vget_low_s8 (v1_0hs)); - const int16x8_t ph0h = vmull_s8(vget_high_s8(v0_0hs), vget_high_s8(v1_0hs)); - const int16x8_t ph1l = vmull_s8(vget_low_s8 (v0_1hs), vget_low_s8 (v1_1hs)); - const int16x8_t ph1h = vmull_s8(vget_high_s8(v0_1hs), vget_high_s8(v1_1hs)); - - const int16x8_t pl0 = vaddq_s16(pl0l, pl0h); - const int16x8_t pl1 = vaddq_s16(pl1l, pl1h); - const int16x8_t ph0 = vaddq_s16(ph0l, ph0h); - const int16x8_t ph1 = vaddq_s16(ph1l, ph1h); - - const int16x8_t pl = vaddq_s16(pl0, pl1); - const int16x8_t ph = vaddq_s16(ph0, ph1); - - //const int8x16_t pl0 = vmulq_s8(v0_0ls, v1_0ls); - //const int8x16_t pl1 = vmulq_s8(v0_1ls, v1_1ls); - //const int8x16_t ph0 = vmulq_s8(v0_0hs, v1_0hs); - //const int8x16_t ph1 = vmulq_s8(v0_1hs, v1_1hs); - - //const int16x8_t pll = vaddl_s8(vget_low_s8(pl0), vget_low_s8(pl1)); - //const int16x8_t plh = vaddl_s8(vget_high_s8(pl0), vget_high_s8(pl1)); - //const int16x8_t phl = vaddl_s8(vget_low_s8(ph0), vget_low_s8(ph1)); - //const int16x8_t phh = vaddl_s8(vget_high_s8(ph0), vget_high_s8(ph1)); - - //const int16x8_t pl = vaddq_s16(pll, plh); - //const int16x8_t ph = vaddq_s16(phl, phh); - - const int16x8_t p = vaddq_s16(pl, ph); - - // convert to float - //const float32x4_t pf0 = vcvtq_f32_s32(vmovl_s16(vget_low_s16 (p))); - //const float32x4_t pf1 = vcvtq_f32_s32(vmovl_s16(vget_high_s16(p))); - - // scalar - sum11 += d0*d1*vaddvq_s16(p); - //sum11 += d0*d1*(vaddvq_s16(pl) + vaddvq_s16(ph)); - //sum11 += d0*d1*vaddvq_s16(vaddq_s16(pl, ph)); - //sum11 += d0*d1*(vaddvq_s8(pl0) + vaddvq_s8(pl1) + vaddvq_s8(ph0) + vaddvq_s8(ph1)); - //sum11 += d0*d1*(vaddvq_s16(pll) + vaddvq_s16(plh) + vaddvq_s16(phl) + vaddvq_s16(phh)); - - //sum_0 = vfmaq_f16(sum_0, d0d1v, vcvtq_f16_s16(p)); - //sum_0 = vfmaq_f16(sum_0, d0d1v, vcvtq_f16_s16(pl)); - //sum_1 = vfmaq_f16(sum_1, d0d1v, vcvtq_f16_s16(ph)); - - // vectorize - //sum_0 = vmlaq_f32(sum_0, d0d1v, pf0); - //sum_1 = vmlaq_f32(sum_1, d0d1v, pf1); - } - - sumf = sum11; - //sumf = vaddvq_f32(sum_0) + vaddvq_f32(sum_1); - //sumf = sum_0[0] + sum_0[1] + sum_0[2] + sum_0[3] + sum_0[4] + sum_0[5] + sum_0[6] + sum_0[7]; - //sum_0 = vaddq_f16(sum_0, sum_1); - //sumf = sum_0[0] + sum_0[1] + sum_0[2] + sum_0[3] + sum_0[4] + sum_0[5] + sum_0[6] + sum_0[7]; -#endif -#endif - - *s = sumf; -} - -// use vec_dot_gq_5 to compute the dot product of two rows -void mul_mat_gq_5( - const void * src0, - const void * src1, // transposed - float * dst, - int m, int n, int k) { - assert(k % QK == 0); - - const int nb = quantize_5_blocks_per_row(k); - - for (int ir0 = 0; ir0 < m; ir0++) { - for (int ir1 = 0; ir1 < n; ir1++) { - vec_dot_gq_5(k, dst + ir1, src0, src1); - src1 = (const char *) src1 + quantize_5_row_size(k); - } - src0 = (const char *) src0 + quantize_5_row_size(k); - src1 = (const char *) src1 - n*quantize_5_row_size(k); - - dst = (float *) dst + n; - } -} - -// -// method 6 -// same as 5 but with 32 element blocks -// - -static inline int quantize_6_blocks_per_row(int k) { - return k/32; -} - -static inline int quantize_6_row_size(int k) { - const int nb = quantize_6_blocks_per_row(k); - - return nb*(sizeof(gq_scale_t) + 16); -} - -void quantize_6_row(const float * restrict src, void * restrict dst, int k) { - assert(k % 32 == 0); - assert(QB == 4); - - const int nb = quantize_6_blocks_per_row(k); - - gq_scale_t * restrict pd = (gq_scale_t *) (dst); - uint8_t * restrict pb = (uint8_t *) (pd + nb); - - uint8_t pp[16]; - - for (int i = 0; i < nb; i++) { - memset(pp, 0, sizeof(pp)); - - float amax = 0.0f; // absolute max - -#if defined(__AVX2__) - { - enum { QK8 = 4 }; - - __m256 srcv [QK8]; - __m256 asrcv[QK8]; - __m256 amaxv[QK8]; - - for (int l = 0; l < QK8; l++) { - srcv[l] = _mm256_loadu_ps(src + i*32 + 8*l); - } - - for (int l = 0; l < QK8; l++) { - asrcv[l] = _mm256_and_ps(srcv[l], _mm256_castsi256_ps(_mm256_set1_epi32(0x7fffffff))); - } - - for (int l = 0; l < QK8/2; l++) { - amaxv[2*l] = _mm256_max_ps(asrcv[2*l], asrcv[2*l+1]); - } - - for (int l = 0; l < QK8/4; l++) { - amaxv[4*l] = _mm256_max_ps(amaxv[4*l], amaxv[4*l+2]); - } - - const __m256 amaxv0_0 = _mm256_permute2f128_ps(amaxv[0], amaxv[0], 3); - const __m256 amaxv0_1 = _mm256_max_ps(amaxv[0], amaxv0_0); - const __m256 amaxv0_2 = _mm256_permute_ps(amaxv0_1, 0x4e); - const __m256 amaxv0_3 = _mm256_max_ps(amaxv0_1, amaxv0_2); - const __m256 amaxv0_4 = _mm256_permute_ps(amaxv0_3, 0xb1); - const __m256 amaxv0_5 = _mm256_max_ps(amaxv0_3, amaxv0_4); - - amax = _mm256_cvtss_f32(amaxv0_5); - - const float d = amax / ((1 << (QB - 1)) - 1); - const float id = d ? 1.0/d : 0.0; - - pd[i] = GGML_FP32_TO_GQ(d); - - const __m256 idv = _mm256_set1_ps(id); - - for (int l = 0; l < 4; l++) { - __m256 v = _mm256_mul_ps(srcv[l], idv); - - // convert to int8 - __m256i vi = _mm256_cvtps_epi32(v); - vi = _mm256_add_epi32(vi, _mm256_set1_epi32(8)); - - int32_t vi_0 = _mm256_extract_epi32(vi, 0); - int32_t vi_1 = _mm256_extract_epi32(vi, 1); - int32_t vi_2 = _mm256_extract_epi32(vi, 2); - int32_t vi_3 = _mm256_extract_epi32(vi, 3); - - int32_t vi_4 = _mm256_extract_epi32(vi, 4); - int32_t vi_5 = _mm256_extract_epi32(vi, 5); - int32_t vi_6 = _mm256_extract_epi32(vi, 6); - int32_t vi_7 = _mm256_extract_epi32(vi, 7); - - // convert to 4-bit, 2 consecutive packed into 1 byte - pp[4*l + 0] = vi_0 | (vi_1 << 4); - pp[4*l + 1] = vi_2 | (vi_3 << 4); - pp[4*l + 2] = vi_4 | (vi_5 << 4); - pp[4*l + 3] = vi_6 | (vi_7 << 4); - - assert(vi_0 >= 0 && vi_0 < 16); - assert(vi_1 >= 0 && vi_1 < 16); - assert(vi_2 >= 0 && vi_2 < 16); - assert(vi_3 >= 0 && vi_3 < 16); - - assert(vi_4 >= 0 && vi_4 < 16); - assert(vi_5 >= 0 && vi_5 < 16); - assert(vi_6 >= 0 && vi_6 < 16); - assert(vi_7 >= 0 && vi_7 < 16); - } - - memcpy(pb + i*16, pp, sizeof(pp)); - } -#elif defined(__ARM_NEON) - { - float32x4_t srcv [8]; - float32x4_t asrcv[8]; - float32x4_t amaxv[8]; - - for (int l = 0; l < 8; l++) srcv[l] = vld1q_f32(src + i*32 + 4*l); - for (int l = 0; l < 8; l++) asrcv[l] = vabsq_f32(srcv[l]); - - for (int l = 0; l < 4; l++) amaxv[2*l] = vmaxq_f32(asrcv[2*l], asrcv[2*l+1]); - for (int l = 0; l < 2; l++) amaxv[4*l] = vmaxq_f32(amaxv[4*l], amaxv[4*l+2]); - for (int l = 0; l < 1; l++) amaxv[8*l] = vmaxq_f32(amaxv[8*l], amaxv[8*l+4]); - - amax = MAX( - MAX(vgetq_lane_f32(amaxv[0], 0), vgetq_lane_f32(amaxv[0], 1)), - MAX(vgetq_lane_f32(amaxv[0], 2), vgetq_lane_f32(amaxv[0], 3))); - - const float d = amax / ((1 << 3) - 1); - const float id = d ? 1.0/d : 0.0; - - pd[i] = GGML_FP32_TO_GQ(d); - - for (int l = 0; l < 8; l++) { - const float32x4_t v = vmulq_n_f32(srcv[l], id); - const float32x4_t vf = vaddq_f32(v, vdupq_n_f32(8.5f)); - const int32x4_t vi = vcvtq_s32_f32(vf); - - pp[2*l + 0] = vgetq_lane_s32(vi, 0) | (vgetq_lane_s32(vi, 1) << 4); - pp[2*l + 1] = vgetq_lane_s32(vi, 2) | (vgetq_lane_s32(vi, 3) << 4); - } - - memcpy(pb + i*16, pp, sizeof(pp)); - } -#else - { - for (int l = 0; l < 32; l++) { - const float v = src[i*32 + l]; - amax = MAX(amax, fabsf(v)); - } - - const float d = amax / ((1 << (QB - 1)) - 1); - const float id = d ? 1.0/d : 0.0; - - pd[i] = GGML_FP32_TO_GQ(d); - - for (int l = 0; l < 32; l++) { - const float v = src[i*32 + l]*id; - const int8_t vi = ((int8_t) (round(v))) + 8; - assert(vi >= 0 && vi < 16); - pp[l/2] |= (vi & 0xf) << (4*(l & 1)); - } - - memcpy(pb + i*16, pp, sizeof(pp)); - } -#endif - //printf("amax = %f\n", amax); - } -} - -// reimplementation of quantize__6using quantize_6_row -void quantize_6(const float * restrict src, char * restrict dst, int n, int k) { - assert(k % 32 == 0); - - for (int j = 0; j < n; j++) { - quantize_6_row(src + j*k, dst, k); - dst = (char *) dst + quantize_6_row_size(k); - } -} - -void vec_dot_gq_6(const int n, float * restrict s, const void * restrict x, const void * restrict y) { - const int nb = quantize_6_blocks_per_row(n); - - const gq_scale_t * restrict pd0 = (const gq_scale_t *) x; - const gq_scale_t * restrict pd1 = (const gq_scale_t *) y; - - const uint8_t * restrict pb0 = (const uint8_t *) (pd0 + nb); - const uint8_t * restrict pb1 = (const uint8_t *) (pd1 + nb); - - float sumf = 0.0; - -#if 0 - // scalar - for (int i = 0; i < nb; i++) { - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - const uint8_t * restrict p0 = pb0 + i*16; - const uint8_t * restrict p1 = pb1 + i*16; - - for (int j = 0; j < 16; j++) { - const uint8_t v0 = p0[j]; - const uint8_t v1 = p1[j]; - - const float f0 = d0*((int8_t) (v0 & 0xf) - 8); - const float f1 = d0*((int8_t) (v0 >> 4) - 8); - - const float f2 = d1*((int8_t) (v1 & 0xf) - 8); - const float f3 = d1*((int8_t) (v1 >> 4) - 8); - - sumf += f0*f2 + f1*f3; - } - } -#else -#if defined(__AVX2__) - // TODO -#elif defined (__ARM_NEON) -#if 0 - float sum0 = 0.0f; - - for (int i = 0; i < nb; i++) { - const float d0 = GGML_GQ_TO_FP32(pd0[i]); - const float d1 = GGML_GQ_TO_FP32(pd1[i]); - - //float32x4_t d0d1v = vdupq_n_f32(d0*d1); - //float16x8_t d0d1v = vdupq_n_f16(d0*d1); - - const uint8_t * restrict p0 = pb0 + i*16; - const uint8_t * restrict p1 = pb1 + i*16; - - const uint8x16_t m4b = vdupq_n_u8(0xf); - const int8x16_t s8b = vdupq_n_s8(0x8); - - const uint8x16_t v0_0 = vld1q_u8(p0); - const uint8x16_t v1_0 = vld1q_u8(p1); - - // 4-bit -> 8-bit - const uint8x16_t v0_0l = vandq_u8(v0_0, m4b); - const uint8x16_t v1_0l = vandq_u8(v1_0, m4b); - - const uint8x16_t v0_0h = vshrq_n_u8(v0_0, 4); - const uint8x16_t v1_0h = vshrq_n_u8(v1_0, 4); - - // sub 8 - const int8x16_t v0_0ls = vsubq_s8(v0_0l, s8b); - const int8x16_t v1_0ls = vsubq_s8(v1_0l, s8b); - - const int8x16_t v0_0hs = vsubq_s8(v0_0h, s8b); - const int8x16_t v1_0hs = vsubq_s8(v1_0h, s8b); - - // dot product into int16x8_t - const int16x8_t pl0l = vmull_s8(vget_low_s8 (v0_0ls), vget_low_s8 (v1_0ls)); - const int16x8_t pl0h = vmull_s8(vget_high_s8(v0_0ls), vget_high_s8(v1_0ls)); - - const int16x8_t ph0l = vmull_s8(vget_low_s8 (v0_0hs), vget_low_s8 (v1_0hs)); - const int16x8_t ph0h = vmull_s8(vget_high_s8(v0_0hs), vget_high_s8(v1_0hs)); - - const int16x8_t pl = vaddq_s16(pl0l, pl0h); - const int16x8_t ph = vaddq_s16(ph0l, ph0h); - - const int16x8_t p = vaddq_s16(pl, ph); - - // scalar - sum0 += d0*d1*vaddvq_s16(p); - } - - sumf = sum0; -#elif 1 // this is a bit faster than the above - float sum0 = 0.0f; - float sum1 = 0.0f; - - for (int i = 0; i < nb; i += 2) { - const float d0_0 = GGML_GQ_TO_FP32(pd0[i + 0]); - const float d1_0 = GGML_GQ_TO_FP32(pd1[i + 0]); - const float d0_1 = GGML_GQ_TO_FP32(pd0[i + 1]); - const float d1_1 = GGML_GQ_TO_FP32(pd1[i + 1]); - - const uint8_t * restrict p0 = pb0 + i*16; - const uint8_t * restrict p1 = pb1 + i*16; - - const uint8x16_t m4b = vdupq_n_u8(0xf); - const int8x16_t s8b = vdupq_n_s8(0x8); - - const uint8x16_t v0_0 = vld1q_u8(p0); - const uint8x16_t v0_1 = vld1q_u8(p0 + 16); - const uint8x16_t v1_0 = vld1q_u8(p1); - const uint8x16_t v1_1 = vld1q_u8(p1 + 16); - - // 4-bit -> 8-bit - const int8x16_t v0_0l = vreinterpretq_s8_u8(vandq_u8(v0_0, m4b)); - const int8x16_t v1_0l = vreinterpretq_s8_u8(vandq_u8(v1_0, m4b)); - - const int8x16_t v0_0h = vreinterpretq_s8_u8(vshrq_n_u8(v0_0, 4)); - const int8x16_t v1_0h = vreinterpretq_s8_u8(vshrq_n_u8(v1_0, 4)); - - const int8x16_t v0_1l = vreinterpretq_s8_u8(vandq_u8(v0_1, m4b)); - const int8x16_t v1_1l = vreinterpretq_s8_u8(vandq_u8(v1_1, m4b)); - - const int8x16_t v0_1h = vreinterpretq_s8_u8(vshrq_n_u8(v0_1, 4)); - const int8x16_t v1_1h = vreinterpretq_s8_u8(vshrq_n_u8(v1_1, 4)); - - // sub 8 - const int8x16_t v0_0ls = vsubq_s8(v0_0l, s8b); - const int8x16_t v1_0ls = vsubq_s8(v1_0l, s8b); - - const int8x16_t v0_0hs = vsubq_s8(v0_0h, s8b); - const int8x16_t v1_0hs = vsubq_s8(v1_0h, s8b); - - const int8x16_t v0_1ls = vsubq_s8(v0_1l, s8b); - const int8x16_t v1_1ls = vsubq_s8(v1_1l, s8b); - - const int8x16_t v0_1hs = vsubq_s8(v0_1h, s8b); - const int8x16_t v1_1hs = vsubq_s8(v1_1h, s8b); - - // dot product into int16x8_t - const int16x8_t pl0l = vmull_s8(vget_low_s8 (v0_0ls), vget_low_s8 (v1_0ls)); - const int16x8_t pl0h = vmull_s8(vget_high_s8(v0_0ls), vget_high_s8(v1_0ls)); - - const int16x8_t ph0l = vmull_s8(vget_low_s8 (v0_0hs), vget_low_s8 (v1_0hs)); - const int16x8_t ph0h = vmull_s8(vget_high_s8(v0_0hs), vget_high_s8(v1_0hs)); - - const int16x8_t pl1l = vmull_s8(vget_low_s8 (v0_1ls), vget_low_s8 (v1_1ls)); - const int16x8_t pl1h = vmull_s8(vget_high_s8(v0_1ls), vget_high_s8(v1_1ls)); - - const int16x8_t ph1l = vmull_s8(vget_low_s8 (v0_1hs), vget_low_s8 (v1_1hs)); - const int16x8_t ph1h = vmull_s8(vget_high_s8(v0_1hs), vget_high_s8(v1_1hs)); - - const int16x8_t pl_0 = vaddq_s16(pl0l, pl0h); - const int16x8_t ph_0 = vaddq_s16(ph0l, ph0h); - - const int16x8_t pl_1 = vaddq_s16(pl1l, pl1h); - const int16x8_t ph_1 = vaddq_s16(ph1l, ph1h); - - const int16x8_t p_0 = vaddq_s16(pl_0, ph_0); - const int16x8_t p_1 = vaddq_s16(pl_1, ph_1); - - // scalar - sum0 += d0_0*d1_0*vaddvq_s16(p_0); - sum1 += d0_1*d1_1*vaddvq_s16(p_1); - } - - sumf = sum0 + sum1; -#endif -#endif -#endif - - *s = sumf; -} - -// use vec_dot_gq_6 to compute the dot product of two rows -void mul_mat_gq_6( - const void * src0, - const void * src1, // transposed - float * dst, - int m, int n, int k) { - assert(k % 32 == 0); - - for (int ir0 = 0; ir0 < m; ir0++) { - for (int ir1 = 0; ir1 < n; ir1++) { - vec_dot_gq_6(k, dst + ir1, src0, src1); - src1 = (const char *) src1 + quantize_6_row_size(k); - } - src0 = (const char *) src0 + quantize_6_row_size(k); - src1 = (const char *) src1 - n*quantize_6_row_size(k); - - dst = (float *) dst + n; - } -} - -int main(int argc, const char ** argv) { - assert(sizeof(gq_quant_t)*8 == gq_t_bits); - ggml_time_init(); - - // needed to initialize f16 tables - { - struct ggml_init_params params = { 0, NULL, false }; - struct ggml_context * ctx = ggml_init(params); - ggml_free(ctx); - } - - int method = 0; - if (argc > 1) { - method = atoi(argv[1]); - } - - float * src0 = malloc(sizeof(float)*M*K); - float * src1 = malloc(sizeof(float)*N*K); - float * dst = malloc(sizeof(float)*M*N); - - // allocate aligned memory - //float * src0 = (float *)aligned_alloc(32, sizeof(float)*M*K); - //float * src1 = (float *)aligned_alloc(32, sizeof(float)*N*K); - //float * dst = (float *)aligned_alloc(32, sizeof(float)*M*N); - - for (int i = 0; i < M*K; i++) { - src0[i] = 0.8 - rand() / (float)RAND_MAX; - /*src0[i] = rand() / (float)RAND_MAX;*/ - /*src0[i] = i % 2;*/ - } - - for (int i = 0; i < N*K; i++) { - src1[i] = 0.8 - rand() / (float)RAND_MAX; - /*src1[i] = rand() / (float)RAND_MAX;*/ - /*src1[i] = i % 3;*/ - } - - void * src0_gq = NULL; - void * src1_gq = NULL; - - size_t sizegq = 0; - - { - if (method == 1) { - src0_gq = calloc(1, quantize_1_row_size(K)*M); - src1_gq = calloc(1, quantize_1_row_size(K)*N); - - sizegq = quantize_1_row_size(K)*M + quantize_1_row_size(K)*N; - } - - if (method == 2) { - src0_gq = calloc(1, quantize_2_row_size(K)*M); - src1_gq = calloc(1, quantize_2_row_size(K)*N); - - sizegq = quantize_2_row_size(K)*M + quantize_2_row_size(K)*N; - } - - if (method == 3) { - src0_gq = calloc(1, quantize_3_row_size(K)*M); - src1_gq = calloc(1, quantize_3_row_size(K)*N); - - sizegq = quantize_3_row_size(K)*M + quantize_3_row_size(K)*N; - } - - if (method == 4) { - src0_gq = calloc(1, quantize_4_row_size(K)*M); - src1_gq = calloc(1, quantize_4_row_size(K)*N); - - sizegq = quantize_4_row_size(K)*M + quantize_4_row_size(K)*N; - } - - if (method == 5) { - src0_gq = calloc(1, quantize_5_row_size(K)*M); - src1_gq = calloc(1, quantize_5_row_size(K)*N); - - sizegq = quantize_5_row_size(K)*M + quantize_5_row_size(K)*N; - } - - if (method == 6) { - src0_gq = calloc(1, quantize_6_row_size(K)*M); - src1_gq = calloc(1, quantize_6_row_size(K)*N); - - sizegq = quantize_6_row_size(K)*M + quantize_6_row_size(K)*N; - } - } - - const size_t sizef16 = sizeof(ggml_fp16_t)*M*K + sizeof(ggml_fp16_t)*N*K; - - printf("compression: %f\n", (float)sizegq/sizef16); - - // convert fp32 -> gq - { - const int64_t t_start = ggml_time_us(); - - if (method == 1) { - quantize_1(src0, src0_gq, M, K); - quantize_1(src1, src1_gq, N, K); - } - - if (method == 2) { - quantize_2(src0, src0_gq, M, K); - quantize_2(src1, src1_gq, N, K); - } - - if (method == 3) { - quantize_3(src0, src0_gq, M, K); - quantize_3(src1, src1_gq, N, K); - } - - if (method == 4) { - quantize_4(src0, src0_gq, M, K); - quantize_4(src1, src1_gq, N, K); - } - - if (method == 5) { - quantize_5(src0, src0_gq, M, K); - quantize_5(src1, src1_gq, N, K); - } - - if (method == 6) { - quantize_6(src0, src0_gq, M, K); - quantize_6(src1, src1_gq, N, K); - } - - const int64_t t_end = ggml_time_us(); - printf("convert time: %f ms / method = %d\n", (t_end - t_start) / 1000.0, method); - } - - for (int i = 0; i < 16; ++i) { - printf("%f %f\n", src0[i], src1[i]); - } - - const int nIter = 1; - - const int64_t start = ggml_cycles(); - const int64_t start_us = ggml_time_us(); - - double iM = 1.0/M; - double sum = 0.0f; - for (int i = 0; i < nIter; i++) { - if (method == 0) { - mul_mat_f32_naive(src0, src1, dst, M, N, K); - } - - if (method == 1) { - mul_mat_gq_1(src0_gq, src1_gq, dst, M, N, K); - } - - if (method == 2) { - mul_mat_gq_2(src0_gq, src1_gq, dst, M, N, K); - } - - if (method == 3) { - mul_mat_gq_3(src0_gq, src1_gq, dst, M, N, K); - } - - if (method == 4) { - mul_mat_gq_4(src0_gq, src1_gq, dst, M, N, K); - } - - if (method == 5) { - mul_mat_gq_5(src0_gq, src1_gq, dst, M, N, K); - } - - if (method == 6) { - mul_mat_gq_6(src0_gq, src1_gq, dst, M, N, K); - } - } - - for (int i = 0; i < N; i++) { - sum += dst[i]*iM; - } - - { - const int64_t end = ggml_cycles(); - const int64_t end_us = ggml_time_us(); - printf("%s: elapsed ticks: %" PRIu64 "\n", __func__, end - start); - printf("%s: elapsed us: %d / %f ms\n", __func__, (int)(end_us - start_us), (end_us - start_us) / 1000.0 / nIter); - } - -#if 0 - // print src0 - printf("src0:\n"); - for (int i = 0; i < M; i++) { - for (int j = 0; j < K; j++) { - printf("%4.1f ", src0[i*K+j]); - } - printf("\n"); - } - - // print src1 - printf("src1:\n"); - for (int i = 0; i < N; i++) { - for (int j = 0; j < K; j++) { - printf("%4.1f ", src1[i*K+j]); - } - printf("\n"); - } - - printf("dst:\n"); - for (int i = 0; i < M; i++) { - for (int j = 0; j < N; j++) { - printf("%4.1f ", dst[i*N+j]); - } - printf("\n"); - } -#endif - - printf("%f\n", sum); - - free(src0); - free(src1); - free(dst); - - if (src0_gq) free(src0_gq); - if (src1_gq) free(src1_gq); - - return 0; -} diff --git a/tests/test-svd0.c b/tests/test-svd0.c deleted file mode 100644 index bae414ed..00000000 --- a/tests/test-svd0.c +++ /dev/null @@ -1,214 +0,0 @@ -// SVD dimensionality reduction - -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -float frand(void) { - return (float) rand() / (float) RAND_MAX; -} - -//int sgesvd_(char *__jobu, char *__jobvt, __CLPK_integer *__m, -// __CLPK_integer *__n, __CLPK_real *__a, __CLPK_integer *__lda, -// __CLPK_real *__s, __CLPK_real *__u, __CLPK_integer *__ldu, -// __CLPK_real *__vt, __CLPK_integer *__ldvt, __CLPK_real *__work, -// __CLPK_integer *__lwork, -// __CLPK_integer *__info) - -int main(int argc, const char ** argv) { - int m = 10; - int n = 5; - - float * A = malloc(n * m * sizeof(float)); - float * A0 = malloc(n * m * sizeof(float)); - - for (int i = 0; i < n; ++i) { - for (int j = 0; j < m; ++j) { - A[i * m + j] = (float) (10.0f*(i + 1) + 1.0f * frand()); - //A[i * m + j] = (float) (10.0f*(i%2 + 1) + 0.1f * frand()); - //if (i == 2) { - // A[i * m + j] += 20*frand(); - //} - if ((i == 1 || i == 3) && j > m/2) { - A[i * m + j] = -A[i * m + j]; - } - } - } - - // average vector - //float * M = malloc(m * sizeof(float)); - - //{ - // for (int j = 0; j < m; ++j) { - // M[j] = 0.0f; - // } - // for (int i = 0; i < n; ++i) { - // for (int j = 0; j < m; ++j) { - // M[j] += A[i * m + j]; - // } - // } - // for (int j = 0; j < m; ++j) { - // M[j] /= (float) n; - // } - //} - - //// subtract average vector - //for (int i = 0; i < n; ++i) { - // for (int j = 0; j < m; ++j) { - // A[i * m + j] -= M[j]; - // } - //} - - memcpy(A0, A, n * m * sizeof(float)); - - // print A - printf("A:\n"); - for (int i = 0; i < n; ++i) { - printf("col %d : ", i); - for (int j = 0; j < m; ++j) { - printf("%9.5f ", A[i * m + j]); - } - printf("\n"); - } - printf("\n"); - - // SVD - // A = U * S * V^T - - float * U = malloc(n * m * sizeof(float)); - float * S = malloc(n * sizeof(float)); - float * V = malloc(n * n * sizeof(float)); - - int lda = m; - int ldu = m; - int ldvt = n; - - float work_size; - int lwork = -1; - int info = 0; - - sgesvd_("S", "S", &m, &n, A, &lda, S, U, &ldu, V, &ldvt, &work_size, &lwork, &info); - - lwork = (int) work_size; - - printf("work_size = %f, info = %d, lwork = %d\n", work_size, info, lwork); - - float * work = malloc(lwork * sizeof(float)); - - sgesvd_("S", "S", &m, &n, A, &lda, S, U, &ldu, V, &ldvt, work, &lwork, &info); - - // print U - printf("U:\n"); - for (int i = 0; i < n; ++i) { - printf("col %d : ", i); - for (int j = 0; j < m; ++j) { - printf("%9.5f ", U[i * m + j]); - } - printf("\n"); - } - printf("\n"); - - // normalize S - { - double sum = 0.0; - for (int i = 0; i < n; ++i) { - sum += S[i]; - } - sum *= sqrt((double) m); - for (int i = 0; i < n; ++i) { - S[i] /= sum; - } - } - - // print S - printf("S:\n"); - for (int i = 0; i < n; ++i) { - printf("- %d = %9.5f\n", i, S[i]); - } - printf("\n"); - - // print V - printf("V:\n"); - for (int i = 0; i < n; ++i) { - printf("col %d : ", i); - for (int j = 0; j < n; ++j) { - printf("%9.5f ", V[i * n + j]); - } - printf("\n"); - } - printf("\n"); - - // print A - printf("A:\n"); - for (int i = 0; i < n; ++i) { - printf("col %d : ", i); - for (int j = 0; j < m; ++j) { - printf("%9.5f ", A[i * m + j]); - } - printf("\n"); - } - printf("\n"); - - // compute singular vectors in U - for (int i = 0; i < n; ++i) { - for (int j = 0; j < m; ++j) { - U[i * m + j] *= S[i]; - } - } - - // normalize U - for (int i = 0; i < n; ++i) { - double sum = 0.0; - for (int j = 0; j < m; ++j) { - sum += U[i * m + j] * U[i * m + j]; - } - sum = sqrt(sum); - for (int j = 0; j < m; ++j) { - U[i * m + j] /= sum*sqrt((double) m); - } - } - - // print U - printf("U:\n"); - for (int i = 0; i < n; ++i) { - printf("col %d : ", i); - for (int j = 0; j < m; ++j) { - printf("%9.5f ", U[i * m + j]); - } - printf("\n"); - } - printf("\n"); - - - // project A0 onto U - float * A1 = malloc(n * n * sizeof(float)); - - for (int i = 0; i < n; ++i) { - for (int j = 0; j < n; ++j) { - A1[i * n + j] = 0.0f; - for (int k = 0; k < m; ++k) { - A1[i * n + j] += A0[i * m + k] * U[j * m + k]; - } - } - } - - // print A1 - printf("A1:\n"); - for (int i = 0; i < n; ++i) { - printf("col %d : ", i); - for (int j = 0; j < n; ++j) { - printf("%9.5f ", A1[i * n + j]); - } - printf("\n"); - } - printf("\n"); - - return 0; -} diff --git a/tests/test-vec0.c b/tests/test-vec0.c deleted file mode 100644 index fc7b8ee9..00000000 --- a/tests/test-vec0.c +++ /dev/null @@ -1,133 +0,0 @@ -#include -#include -#include -#include - -const int N = 1 << 14; -const int M = 1 << 14; - -void mul_mat_vec_f32_0( - const float * src0, - const float * src1, - float * dst, - unsigned nrows, - unsigned ncols) { - for (unsigned i = 0; i < nrows; i++) { - float sum = 0.0f; - for (unsigned j = 0; j < ncols; j++) { - sum += src0[i*ncols + j]*src1[j]; - } - dst[i] = sum; - } -} -#if defined(_MSC_VER) -typedef float __declspec(align(32)) afloat; -#else -typedef float afloat __attribute__((__aligned__(32))); -#endif -void mul_mat_vec_f32_1( - const afloat *restrict src0, - const afloat *restrict src1, - afloat *restrict dst, - unsigned nrows, - unsigned ncols) { - for (unsigned i = 0; i < nrows; i++) { - const afloat * restrict row = src0 + i*ncols; - const afloat * restrict col = src1; - - float sum = 0.0f; - - for (unsigned j = 0; j < ncols; j++) { - sum += *row++ * *col++; - } - - dst[i] = sum; - - //float sum[8] = {0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f}; - - //for (unsigned j = 0; j < ncols; j += 8) { - // sum[0] += row[0]*col[0]; - // sum[1] += row[1]*col[1]; - // sum[2] += row[2]*col[2]; - // sum[3] += row[3]*col[3]; - // sum[4] += row[4]*col[4]; - // sum[5] += row[5]*col[5]; - // sum[6] += row[6]*col[6]; - // sum[7] += row[7]*col[7]; - - // row += 8; - // col += 8; - //} - - //dst[i] = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5] + sum[6] + sum[7]; - } -} - -void mul_mat_vec_f32_2( - const void * src0, - const void * src1, - void * dst, - unsigned nrows, - unsigned ncols) { - void * d = dst; - for (unsigned i = 0; i < nrows; i++) { - float sum = 0.0f; - - const char * row = (const char*)src0 + i*ncols*sizeof(float); - const char * col = (const char*)src1; - for (unsigned j = 0; j < ncols; j++) { - sum += (*(float *)row) * (*(float *)col); - row += sizeof(float); - col += sizeof(float); - } - *(float *)d = sum; - d = (char*)d + sizeof(float); - } -} - -#if defined(_MSC_VER) || defined(__MINGW32__) || defined(__MINGW64__) -void* aligned_alloc(size_t alignment, size_t size) { - return _aligned_malloc(size, alignment); -} -#endif - -int main(int argc, const char ** argv) { - //float * src0 = malloc(sizeof(float)*N*M); - //float * src1 = malloc(sizeof(float)*M); - //float * dst = malloc(sizeof(float)*N); - - afloat * src0 = (float *)(aligned_alloc(32, sizeof(float)*N*M)); - afloat * src1 = (float *)(aligned_alloc(32, sizeof(float)*M)); - afloat * dst = (float *)(aligned_alloc(32, sizeof(float)*N)); - - for (int i = 0; i < N*M; i++) { - src0[i] = (afloat)i; - } - - for (int i = 0; i < M; i++) { - src1[i] = (afloat)i; - } - - const int nIter = 10; - - const clock_t start = clock(); - - double sum = 0.0f; - for (int i = 0; i < nIter; i++) { - //mul_mat_vec_f32_0(src0, src1, dst, N, M); - mul_mat_vec_f32_1(src0, src1, dst, N, M); - //mul_mat_vec_f32_2(src0, src1, dst, N, M); - for (int i = 0; i < N; i++) { - sum += dst[i]; - } - } - - { - const clock_t end = clock(); - printf("%s: elapsed ticks: %ld\n", __func__, end - start); - } - - printf("%f\n", sum); - - return 0; -} diff --git a/tests/test-vec1.c b/tests/test-vec1.c deleted file mode 100644 index 567cb061..00000000 --- a/tests/test-vec1.c +++ /dev/null @@ -1,576 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include - -#include - -const int N = 1 << 14; -const int M = 768; - -// -// naive implementation -// - -void mul_mat_vec_f32_0( - const float * restrict src0, - const float * restrict src1, - float * dst, - int nrows, - int ncols) { - for (int i = 0; i < nrows; i++) { - float sum = 0.0f; - for (int j = 0; j < ncols; j++) { - sum += src0[i*ncols + j]*src1[j]; - } - dst[i] = sum; - } -} - -// -// SIMD with 8 32-bit floats -// - -float reduce_vector8_0(__m256 v) { - __m128 v1 = _mm256_extractf128_ps(v, 0); - __m128 v2 = _mm256_extractf128_ps(v, 1); - __m128 v3 = _mm_add_ps(v1, v2); - __m128 v4 = _mm_shuffle_ps(v3, v3, 0x4e); - __m128 v5 = _mm_add_ps(v3, v4); - __m128 v6 = _mm_shuffle_ps(v5, v5, 0x11); - __m128 v7 = _mm_add_ps(v5, v6); - return _mm_cvtss_f32(v7); -} - -// vectorized implementation using AVX -void mul_mat_vec_f32_1( - const float * restrict src0, - const float * restrict src1, - float * dst, - int nrows, - int ncols) { - - const int ncols8 = ncols & ~7; - - for (int i = 0; i < nrows; i++) { - __m256 sum = _mm256_setzero_ps(); - for (int j = 0; j < ncols8; j += 8) { - __m256 a = _mm256_loadu_ps(src0 + i*ncols + j); - __m256 b = _mm256_loadu_ps(src1 + j); - __m256 c = _mm256_mul_ps(a, b); - sum = _mm256_add_ps(sum, c); - } - dst[i] = reduce_vector8_0(sum); - - for (int j = ncols8; j < ncols; j++) { - dst[i] += src0[i*ncols + j]*src1[j]; - } - } -} - -void mul_mat_vec_f32_2( - const float * restrict src0, - const float * restrict src1, - float * dst, - int nrows, - int ncols) { - - const int ncols32 = ncols & ~31; - - for (int i = 0; i < nrows; i++) { - __m256 sum0 = _mm256_setzero_ps(); - __m256 sum1 = _mm256_setzero_ps(); - __m256 sum2 = _mm256_setzero_ps(); - __m256 sum3 = _mm256_setzero_ps(); - - const float * restrict src0_row = src0 + i*ncols; - for (int j = 0; j < ncols32; j += 32) { - __m256 a0 = _mm256_loadu_ps(src0_row + j + 0); - __m256 a1 = _mm256_loadu_ps(src0_row + j + 8); - __m256 a2 = _mm256_loadu_ps(src0_row + j + 16); - __m256 a3 = _mm256_loadu_ps(src0_row + j + 24); - __m256 b0 = _mm256_loadu_ps(src1 + j + 0); - __m256 b1 = _mm256_loadu_ps(src1 + j + 8); - __m256 b2 = _mm256_loadu_ps(src1 + j + 16); - __m256 b3 = _mm256_loadu_ps(src1 + j + 24); -#if defined(__FMA__) - sum0 = _mm256_fmadd_ps(a0, b0, sum0); - sum1 = _mm256_fmadd_ps(a1, b1, sum1); - sum2 = _mm256_fmadd_ps(a2, b2, sum2); - sum3 = _mm256_fmadd_ps(a3, b3, sum3); -#else - sum0 = _mm256_add_ps(_mm256_mul_ps(a0, b0), sum0); - sum1 = _mm256_add_ps(_mm256_mul_ps(a1, b1), sum1); - sum2 = _mm256_add_ps(_mm256_mul_ps(a2, b2), sum2); - sum3 = _mm256_add_ps(_mm256_mul_ps(a3, b3), sum3); -#endif - } - dst[i] = reduce_vector8_0(_mm256_add_ps(_mm256_add_ps(sum0, sum1), _mm256_add_ps(sum2, sum3))); - - for (int j = ncols32; j < ncols; j++) { - dst[i] += src0[i*ncols + j]*src1[j]; - } - } -} - -// -// SIMD with 8 16-bit floats -// - -static inline float fp32_from_bits(uint32_t w) { -#if defined(__OPENCL_VERSION__) - return as_float(w); -#elif defined(__CUDA_ARCH__) - return __uint_as_float((unsigned int) w); -#elif defined(__INTEL_COMPILER) - return _castu32_f32(w); -#elif defined(_MSC_VER) && (defined(_M_ARM) || defined(_M_ARM64)) - return _CopyFloatFromInt32((__int32) w); -#else - union { - uint32_t as_bits; - float as_value; - } fp32 = { w }; - return fp32.as_value; -#endif -} - -static inline uint32_t fp32_to_bits(float f) { -#if defined(__OPENCL_VERSION__) - return as_uint(f); -#elif defined(__CUDA_ARCH__) - return (uint32_t) __float_as_uint(f); -#elif defined(__INTEL_COMPILER) - return _castf32_u32(f); -#elif defined(_MSC_VER) && (defined(_M_ARM) || defined(_M_ARM64)) - return (uint32_t) _CopyInt32FromFloat(f); -#else - union { - float as_value; - uint32_t as_bits; - } fp32 = { f }; - return fp32.as_bits; -#endif -} - -/* - * Convert a 16-bit floating-point number in IEEE half-precision format, in bit representation, to - * a 32-bit floating-point number in IEEE single-precision format. - * - * @note The implementation relies on IEEE-like (no assumption about rounding mode and no operations on denormals) - * floating-point operations and bitcasts between integer and floating-point variables. - */ -static inline float fp16_ieee_to_fp32_value(uint16_t h) { - /* - * Extend the half-precision floating-point number to 32 bits and shift to the upper part of the 32-bit word: - * +---+-----+------------+-------------------+ - * | S |EEEEE|MM MMMM MMMM|0000 0000 0000 0000| - * +---+-----+------------+-------------------+ - * Bits 31 26-30 16-25 0-15 - * - * S - sign bit, E - bits of the biased exponent, M - bits of the mantissa, 0 - zero bits. - */ - const uint32_t w = (uint32_t) h << 16; - /* - * Extract the sign of the input number into the high bit of the 32-bit word: - * - * +---+----------------------------------+ - * | S |0000000 00000000 00000000 00000000| - * +---+----------------------------------+ - * Bits 31 0-31 - */ - const uint32_t sign = w & UINT32_C(0x80000000); - /* - * Extract mantissa and biased exponent of the input number into the high bits of the 32-bit word: - * - * +-----+------------+---------------------+ - * |EEEEE|MM MMMM MMMM|0 0000 0000 0000 0000| - * +-----+------------+---------------------+ - * Bits 27-31 17-26 0-16 - */ - const uint32_t two_w = w + w; - - /* - * Shift mantissa and exponent into bits 23-28 and bits 13-22 so they become mantissa and exponent - * of a single-precision floating-point number: - * - * S|Exponent | Mantissa - * +-+---+-----+------------+----------------+ - * |0|000|EEEEE|MM MMMM MMMM|0 0000 0000 0000| - * +-+---+-----+------------+----------------+ - * Bits | 23-31 | 0-22 - * - * Next, there are some adjustments to the exponent: - * - The exponent needs to be corrected by the difference in exponent bias between single-precision and half-precision - * formats (0x7F - 0xF = 0x70) - * - Inf and NaN values in the inputs should become Inf and NaN values after conversion to the single-precision number. - * Therefore, if the biased exponent of the half-precision input was 0x1F (max possible value), the biased exponent - * of the single-precision output must be 0xFF (max possible value). We do this correction in two steps: - * - First, we adjust the exponent by (0xFF - 0x1F) = 0xE0 (see exp_offset below) rather than by 0x70 suggested - * by the difference in the exponent bias (see above). - * - Then we multiply the single-precision result of exponent adjustment by 2**(-112) to reverse the effect of - * exponent adjustment by 0xE0 less the necessary exponent adjustment by 0x70 due to difference in exponent bias. - * The floating-point multiplication hardware would ensure than Inf and NaN would retain their value on at least - * partially IEEE754-compliant implementations. - * - * Note that the above operations do not handle denormal inputs (where biased exponent == 0). However, they also do not - * operate on denormal inputs, and do not produce denormal results. - */ - const uint32_t exp_offset = UINT32_C(0xE0) << 23; -#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L) || defined(__GNUC__) && !defined(__STRICT_ANSI__) - const float exp_scale = 0x1.0p-112f; -#else - const float exp_scale = fp32_from_bits(UINT32_C(0x7800000)); -#endif - const float normalized_value = fp32_from_bits((two_w >> 4) + exp_offset) * exp_scale; - - /* - * Convert denormalized half-precision inputs into single-precision results (always normalized). - * Zero inputs are also handled here. - * - * In a denormalized number the biased exponent is zero, and mantissa has on-zero bits. - * First, we shift mantissa into bits 0-9 of the 32-bit word. - * - * zeros | mantissa - * +---------------------------+------------+ - * |0000 0000 0000 0000 0000 00|MM MMMM MMMM| - * +---------------------------+------------+ - * Bits 10-31 0-9 - * - * Now, remember that denormalized half-precision numbers are represented as: - * FP16 = mantissa * 2**(-24). - * The trick is to construct a normalized single-precision number with the same mantissa and thehalf-precision input - * and with an exponent which would scale the corresponding mantissa bits to 2**(-24). - * A normalized single-precision floating-point number is represented as: - * FP32 = (1 + mantissa * 2**(-23)) * 2**(exponent - 127) - * Therefore, when the biased exponent is 126, a unit change in the mantissa of the input denormalized half-precision - * number causes a change of the constructud single-precision number by 2**(-24), i.e. the same ammount. - * - * The last step is to adjust the bias of the constructed single-precision number. When the input half-precision number - * is zero, the constructed single-precision number has the value of - * FP32 = 1 * 2**(126 - 127) = 2**(-1) = 0.5 - * Therefore, we need to subtract 0.5 from the constructed single-precision number to get the numerical equivalent of - * the input half-precision number. - */ - const uint32_t magic_mask = UINT32_C(126) << 23; - const float magic_bias = 0.5f; - const float denormalized_value = fp32_from_bits((two_w >> 17) | magic_mask) - magic_bias; - - /* - * - Choose either results of conversion of input as a normalized number, or as a denormalized number, depending on the - * input exponent. The variable two_w contains input exponent in bits 27-31, therefore if its smaller than 2**27, the - * input is either a denormal number, or zero. - * - Combine the result of conversion of exponent and mantissa with the sign of the input number. - */ - const uint32_t denormalized_cutoff = UINT32_C(1) << 27; - const uint32_t result = sign | - (two_w < denormalized_cutoff ? fp32_to_bits(denormalized_value) : fp32_to_bits(normalized_value)); - return fp32_from_bits(result); -} - -/* - * Convert a 32-bit floating-point number in IEEE single-precision format to a 16-bit floating-point number in - * IEEE half-precision format, in bit representation. - * - * @note The implementation relies on IEEE-like (no assumption about rounding mode and no operations on denormals) - * floating-point operations and bitcasts between integer and floating-point variables. - */ -static inline uint16_t fp16_ieee_from_fp32_value(float f) { -#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L) || defined(__GNUC__) && !defined(__STRICT_ANSI__) - const float scale_to_inf = 0x1.0p+112f; - const float scale_to_zero = 0x1.0p-110f; -#else - const float scale_to_inf = fp32_from_bits(UINT32_C(0x77800000)); - const float scale_to_zero = fp32_from_bits(UINT32_C(0x08800000)); -#endif - float base = (fabsf(f) * scale_to_inf) * scale_to_zero; - - const uint32_t w = fp32_to_bits(f); - const uint32_t shl1_w = w + w; - const uint32_t sign = w & UINT32_C(0x80000000); - uint32_t bias = shl1_w & UINT32_C(0xFF000000); - if (bias < UINT32_C(0x71000000)) { - bias = UINT32_C(0x71000000); - } - - base = fp32_from_bits((bias >> 1) + UINT32_C(0x07800000)) + base; - const uint32_t bits = fp32_to_bits(base); - const uint32_t exp_bits = (bits >> 13) & UINT32_C(0x00007C00); - const uint32_t mantissa_bits = bits & UINT32_C(0x00000FFF); - const uint32_t nonsign = exp_bits + mantissa_bits; - return (sign >> 16) | (shl1_w > UINT32_C(0xFF000000) ? UINT16_C(0x7E00) : nonsign); -} - -void mul_mat_vec_f16_0( - const uint16_t * src0, - const uint16_t * src1, - float * dst, - int nrows, - int ncols) { - - const int ncols8 = ncols & ~7; - - for (int i = 0; i < nrows; i++) { - __m256 sum = _mm256_setzero_ps(); - - const uint16_t * src0_row = src0 + i * ncols; - for (int j = 0; j < ncols8; j += 8) { - __m256 a = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j))); - __m256 b = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src1 + j))); -#if defined(__FMA__) - sum = _mm256_fmadd_ps(a, b, sum); -#else - sum = _mm256_add_ps(_mm256_mul_ps(a, b), sum); -#endif - } - dst[i] = reduce_vector8_0(sum); - - for (int j = ncols8; j < ncols; j++) { - dst[i] += fp16_ieee_to_fp32_value(src0_row[j]) * fp16_ieee_to_fp32_value(src1[j]); - } - } -} - -void mul_mat_vec_f16_1( - const uint16_t * src0, - const uint16_t * src1, - float * dst, - int nrows, - int ncols) { - - const int ncols16 = ncols & ~15; - - for (int i = 0; i < nrows; i++) { - __m256 sum0 = _mm256_setzero_ps(); - __m256 sum1 = _mm256_setzero_ps(); - - const uint16_t * src0_row = src0 + i * ncols; - for (int j = 0; j < ncols16; j += 16) { - __m256 a0 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 0))); - __m256 a1 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 8))); - __m256 b0 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src1 + j))); - __m256 b1 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src1 + j + 8))); -#if defined(__FMA__) - sum0 = _mm256_fmadd_ps(a0, b0, sum0); - sum1 = _mm256_fmadd_ps(a1, b1, sum1); -#else - sum0 = _mm256_add_ps(_mm256_mul_ps(a0, b0), sum0); - sum1 = _mm256_add_ps(_mm256_mul_ps(a1, b1), sum1); -#endif - } - dst[i] = reduce_vector8_0(sum0) + reduce_vector8_0(sum1); - - for (int j = ncols16; j < ncols; j++) { - dst[i] += fp16_ieee_to_fp32_value(src0_row[j]) * fp16_ieee_to_fp32_value(src1[j]); - } - } -} - -void mul_mat_vec_f16_2( - const uint16_t * src0, - const uint16_t * src1, - float * dst, - int nrows, - int ncols) { - - const int ncols32 = ncols & ~31; - - for (int i = 0; i < nrows; i++) { - __m256 sum0 = _mm256_setzero_ps(); - __m256 sum1 = _mm256_setzero_ps(); - __m256 sum2 = _mm256_setzero_ps(); - __m256 sum3 = _mm256_setzero_ps(); - - const uint16_t * src0_row = src0 + i * ncols; - for (int j = 0; j < ncols32; j += 32) { - __m256 a0 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 0))); - __m256 a1 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 8))); - __m256 a2 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 16))); - __m256 a3 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 24))); - __m256 b0 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src1 + j))); - __m256 b1 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src1 + j + 8))); - __m256 b2 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src1 + j + 16))); - __m256 b3 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src1 + j + 24))); -#if defined(__FMA__) - sum0 = _mm256_fmadd_ps(a0, b0, sum0); - sum1 = _mm256_fmadd_ps(a1, b1, sum1); - sum2 = _mm256_fmadd_ps(a2, b2, sum2); - sum3 = _mm256_fmadd_ps(a3, b3, sum3); -#else - sum0 = _mm256_add_ps(_mm256_mul_ps(a0, b0), sum0); - sum1 = _mm256_add_ps(_mm256_mul_ps(a1, b1), sum1); - sum2 = _mm256_add_ps(_mm256_mul_ps(a2, b2), sum2); - sum3 = _mm256_add_ps(_mm256_mul_ps(a3, b3), sum3); -#endif - } - dst[i] = reduce_vector8_0(sum0) + reduce_vector8_0(sum1) + reduce_vector8_0(sum2) + reduce_vector8_0(sum3); - - for (int j = ncols32; j < ncols; j++) { - dst[i] += fp16_ieee_to_fp32_value(src0_row[j]) * fp16_ieee_to_fp32_value(src1[j]); - } - } -} - -void mul_mat_vec_f16_3( - const uint16_t * src0, - const float * src1, - float * dst, - int nrows, - int ncols) { - - const int ncols32 = ncols & ~31; - - for (int i = 0; i < nrows; i++) { - __m256 sum0 = _mm256_setzero_ps(); - __m256 sum1 = _mm256_setzero_ps(); - __m256 sum2 = _mm256_setzero_ps(); - __m256 sum3 = _mm256_setzero_ps(); - - const uint16_t * src0_row = src0 + i * ncols; - for (int j = 0; j < ncols32; j += 32) { - __m256 a0 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 0))); - __m256 a1 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 8))); - __m256 a2 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 16))); - __m256 a3 = _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)(src0_row + j + 24))); - __m256 b0 = _mm256_loadu_ps(src1 + j); - __m256 b1 = _mm256_loadu_ps(src1 + j + 8); - __m256 b2 = _mm256_loadu_ps(src1 + j + 16); - __m256 b3 = _mm256_loadu_ps(src1 + j + 24); -#if defined(__FMA__) - sum0 = _mm256_fmadd_ps(a0, b0, sum0); - sum1 = _mm256_fmadd_ps(a1, b1, sum1); - sum2 = _mm256_fmadd_ps(a2, b2, sum2); - sum3 = _mm256_fmadd_ps(a3, b3, sum3); -#else - sum0 = _mm256_add_ps(_mm256_mul_ps(a0, b0), sum0); - sum1 = _mm256_add_ps(_mm256_mul_ps(a1, b1), sum1); - sum2 = _mm256_add_ps(_mm256_mul_ps(a2, b2), sum2); - sum3 = _mm256_add_ps(_mm256_mul_ps(a3, b3), sum3); -#endif - } - dst[i] = reduce_vector8_0(sum0) + reduce_vector8_0(sum1) + reduce_vector8_0(sum2) + reduce_vector8_0(sum3); - - for (int j = ncols32; j < ncols; j++) { - dst[i] += fp16_ieee_to_fp32_value(src0_row[j]) * fp16_ieee_to_fp32_value(src1[j]); - } - } -} - -uint64_t get_time_us(void) { - struct timeval tv; - gettimeofday(&tv, NULL); - return tv.tv_sec * 1000000 + tv.tv_usec; -} - -int main(int argc, const char ** argv) { - float * src0 = malloc(sizeof(float)*N*M); - float * src1 = malloc(sizeof(float)*M); - float * dst = malloc(sizeof(float)*N); - - //float * src0 = (float *)(aligned_alloc(64, sizeof(float)*N*M)); - //float * src1 = (float *)(aligned_alloc(64, sizeof(float)*M)); - //float * dst = (float *)(aligned_alloc(64, sizeof(float)*N)); - - for (int i = 0; i < N*M; i++) { - src0[i] = rand() / (float)RAND_MAX; - } - - for (int i = 0; i < M; i++) { - src1[i] = rand() / (float)RAND_MAX; - } - - // convert src0 and src1 to __fp16 - uint16_t * src0_fp16 = (uint16_t *)(malloc(sizeof(uint16_t)*N*M)); - uint16_t * src1_fp16 = (uint16_t *)(malloc(sizeof(uint16_t)*M)); - //uint16_t * src0_fp16 = (uint16_t *)(aligned_alloc(64, sizeof(uint16_t)*N*M)); - //uint16_t * src1_fp16 = (uint16_t *)(aligned_alloc(64, sizeof(uint16_t)*M)); - - { - const uint64_t t_start = get_time_us(); - - for (int i = 0; i < N*M; i++) { - src0_fp16[i] = fp16_ieee_from_fp32_value(src0[i]); - //printf("%f %f\n", src0[i], fp16_ieee_to_fp32_value(src0_fp16[i])); - //assert(!isnan(fp16_ieee_to_fp32_value(src0_fp16[i]))); - } - - for (int i = 0; i < M; i++) { - src1_fp16[i] = fp16_ieee_from_fp32_value(src1[i]); - } - - const uint64_t t_end = get_time_us(); - printf("convert time: %f ms\n", (t_end - t_start) / 1000.0); - } - - for (int i = 0; i < 16; ++i) { - printf("%f %f\n", src0[i], fp16_ieee_to_fp32_value(src0_fp16[i])); - } - - int method = 0; - if (argc > 1) { - method = atoi(argv[1]); - } - - const int nIter = 1000; - - const clock_t start = clock(); - const uint64_t start_us = get_time_us(); - - double iM = 1.0/M; - double sum = 0.0f; - for (int i = 0; i < nIter; i++) { - if (method == 0) { - mul_mat_vec_f32_0(src0, src1, dst, N, M); - } - - if (method == 1) { - mul_mat_vec_f32_1(src0, src1, dst, N, M); - } - - if (method == 2) { - mul_mat_vec_f32_2(src0, src1, dst, N, M); - } - - if (method == 3) { - mul_mat_vec_f16_0(src0_fp16, src1_fp16, dst, N, M); - } - - if (method == 4) { - mul_mat_vec_f16_1(src0_fp16, src1_fp16, dst, N, M); - } - - if (method == 5) { - mul_mat_vec_f16_2(src0_fp16, src1_fp16, dst, N, M); - } - - if (method == 6) { - mul_mat_vec_f16_3(src0_fp16, src1, dst, N, M); - } - } - - for (int i = 0; i < N; i++) { - sum += dst[i]*iM; - } - - { - const clock_t end = clock(); - const uint64_t end_us = get_time_us(); - printf("%s: elapsed ticks: %ld\n", __func__, end - start); - printf("%s: elapsed us: %ld\n", __func__, end_us - start_us); - } - - printf("%f\n", sum); - - free(src0); - free(src1); - free(dst); - - free(src0_fp16); - free(src1_fp16); - - return 0; -} diff --git a/tests/test-vec2.c b/tests/test-vec2.c deleted file mode 100644 index 4fa364ca..00000000 --- a/tests/test-vec2.c +++ /dev/null @@ -1,268 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include - -#include - -const int N = 1 << 12; -const int M = 1 << 12; - -// -// naive implementation -// - -void mul_mat_vec_f32_0( - const float * restrict src0, - const float * restrict src1, - float * dst, - int nrows, - int ncols) { - for (int i = 0; i < nrows; i++) { - float sum = 0.0f; - for (int j = 0; j < ncols; j++) { - sum += src0[i*ncols + j]*src1[j]; - } - dst[i] = sum; - } -} - -void mul_mat_vec_f16_0( - const __fp16 * src0, - const __fp16 * src1, - float * dst, - int nrows, - int ncols) { - - const int n64 = ncols & ~63; - - for (int r = 0; r < nrows; r++) { - float sumf = 0.0; - - float16x8_t sum0 = vdupq_n_f16(0.0f); - float16x8_t sum1 = vdupq_n_f16(0.0f); - float16x8_t sum2 = vdupq_n_f16(0.0f); - float16x8_t sum3 = vdupq_n_f16(0.0f); - float16x8_t sum4 = vdupq_n_f16(0.0f); - float16x8_t sum5 = vdupq_n_f16(0.0f); - float16x8_t sum6 = vdupq_n_f16(0.0f); - float16x8_t sum7 = vdupq_n_f16(0.0f); - - float16x8_t x0, x1, x2, x3, x4, x5, x6, x7; - float16x8_t y0, y1, y2, y3, y4, y5, y6, y7; - - const __fp16 * restrict p0 = src0 + r*ncols; - - for (int i = 0; i < n64; i += 64) { - x0 = vld1q_f16(p0 + i + 0 ); - x1 = vld1q_f16(p0 + i + 8 ); - x2 = vld1q_f16(p0 + i + 16); - x3 = vld1q_f16(p0 + i + 24); - x4 = vld1q_f16(p0 + i + 32); - x5 = vld1q_f16(p0 + i + 40); - x6 = vld1q_f16(p0 + i + 48); - x7 = vld1q_f16(p0 + i + 56); - - y0 = vld1q_f16(src1 + i + 0 ); - y1 = vld1q_f16(src1 + i + 8 ); - y2 = vld1q_f16(src1 + i + 16); - y3 = vld1q_f16(src1 + i + 24); - y4 = vld1q_f16(src1 + i + 32); - y5 = vld1q_f16(src1 + i + 40); - y6 = vld1q_f16(src1 + i + 48); - y7 = vld1q_f16(src1 + i + 56); - - sum0 = vfmaq_f16(sum0, x0, y0); - sum1 = vfmaq_f16(sum1, x1, y1); - sum2 = vfmaq_f16(sum2, x2, y2); - sum3 = vfmaq_f16(sum3, x3, y3); - sum4 = vfmaq_f16(sum4, x4, y4); - sum5 = vfmaq_f16(sum5, x5, y5); - sum6 = vfmaq_f16(sum6, x6, y6); - sum7 = vfmaq_f16(sum7, x7, y7); - } - - // TODO: F16 - better way to reduce this ? - float16x8_t sum = vaddq_f16(sum0, sum1); - - sum = vaddq_f16(sum, sum2); - sum = vaddq_f16(sum, sum3); - sum = vaddq_f16(sum, sum4); - sum = vaddq_f16(sum, sum5); - sum = vaddq_f16(sum, sum6); - sum = vaddq_f16(sum, sum7); - - sumf += sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5] + sum[6] + sum[7]; - - for (int j = n64; j < n64; j++) { - sumf += src0[r*ncols + j]*src1[j]; - } - - dst[r] = sumf; - } -} - -void mul_mat_vec_f16_1( - const __fp16 * src0, - const __fp16 * src1, - float * dst, - int nrows, - int ncols) { - - const int n32 = ncols & ~31; - - for (int r = 0; r < nrows; r++) { - float sumf = 0.0; - - float16x8_t sum0 = vdupq_n_f16(0.0f); - float16x8_t sum1 = vdupq_n_f16(0.0f); - float16x8_t sum2 = vdupq_n_f16(0.0f); - float16x8_t sum3 = vdupq_n_f16(0.0f); - - float16x8_t x0, x1, x2, x3; - float16x8_t y0, y1, y2, y3; - - const __fp16 * restrict p0 = src0 + r*ncols; - - for (int i = 0; i < n32; i += 32) { - x0 = vld1q_f16(p0 + i + 0 ); - x1 = vld1q_f16(p0 + i + 8 ); - x2 = vld1q_f16(p0 + i + 16); - x3 = vld1q_f16(p0 + i + 24); - - y0 = vld1q_f16(src1 + i + 0 ); - y1 = vld1q_f16(src1 + i + 8 ); - y2 = vld1q_f16(src1 + i + 16); - y3 = vld1q_f16(src1 + i + 24); - - sum0 = vfmaq_f16(sum0, x0, y0); - sum1 = vfmaq_f16(sum1, x1, y1); - sum2 = vfmaq_f16(sum2, x2, y2); - sum3 = vfmaq_f16(sum3, x3, y3); - } - - // reduce sum0..sum3 to sum0 - sum0 = vaddq_f16(sum0, sum1); - sum2 = vaddq_f16(sum2, sum3); - sum0 = vaddq_f16(sum0, sum2); - - // load sum0 into 2 float32x4_t - float32x4_t sum0f32 = vcvt_f32_f16(vget_low_f16(sum0)); - float32x4_t sum1f32 = vcvt_f32_f16(vget_high_f16(sum0)); - - // reduce sum0f32 and sum1f32 to sumf - sum0f32 = vaddq_f32(sum0f32, sum1f32); - - float32x2_t sumf32 = vadd_f32(vget_low_f32(sum0f32), vget_high_f32(sum0f32)); - sumf = vget_lane_f32(sumf32, 0) + vget_lane_f32(sumf32, 1); - - //sumf = sum0[0] + sum0[1] + sum0[2] + sum0[3] + sum0[4] + sum0[5] + sum0[6] + sum0[7]; - - for (int j = n32; j < n32; j++) { - sumf += src0[r*ncols + j]*src1[j]; - } - - dst[r] = sumf; - } -} - -uint64_t get_time_us(void) { - struct timeval tv; - gettimeofday(&tv, NULL); - return tv.tv_sec * 1000000 + tv.tv_usec; -} - -int main(int argc, const char ** argv) { - float * src0 = malloc(sizeof(float)*N*M); - float * src1 = malloc(sizeof(float)*M); - float * dst = malloc(sizeof(float)*N); - - //float * src0 = (float *)(aligned_alloc(64, sizeof(float)*N*M)); - //float * src1 = (float *)(aligned_alloc(64, sizeof(float)*M)); - //float * dst = (float *)(aligned_alloc(64, sizeof(float)*N)); - - for (int i = 0; i < N*M; i++) { - src0[i] = rand() / (float)RAND_MAX; - } - - for (int i = 0; i < M; i++) { - src1[i] = rand() / (float)RAND_MAX; - } - - // convert src0 and src1 to __fp16 - __fp16 * src0_fp16 = (__fp16 *)(malloc(sizeof(__fp16)*N*M)); - __fp16 * src1_fp16 = (__fp16 *)(malloc(sizeof(__fp16)*M)); - - { - const uint64_t t_start = get_time_us(); - - for (int i = 0; i < N*M; i++) { - src0_fp16[i] = src0[i]; - //printf("%f %f\n", src0[i], src0_fp16[i]); - //assert(!isnan(src0_fp16[i])); - } - - for (int i = 0; i < M; i++) { - src1_fp16[i] = src1[i]; - } - - const uint64_t t_end = get_time_us(); - printf("convert time: %f ms\n", (t_end - t_start) / 1000.0); - } - - for (int i = 0; i < 16; ++i) { - printf("%f %f\n", src0[i], src0_fp16[i]); - } - - int method = 0; - if (argc > 1) { - method = atoi(argv[1]); - } - - const int nIter = 1000; - - const clock_t start = clock(); - const uint64_t start_us = get_time_us(); - - double iM = 1.0/M; - double sum = 0.0f; - for (int i = 0; i < nIter; i++) { - if (method == 0) { - mul_mat_vec_f32_0(src0, src1, dst, N, M); - } - - if (method == 1) { - mul_mat_vec_f16_0(src0_fp16, src1_fp16, dst, N, M); - } - - if (method == 2) { - mul_mat_vec_f16_1(src0_fp16, src1_fp16, dst, N, M); - } - } - - for (int i = 0; i < N; i++) { - sum += dst[i]*iM; - } - - { - const clock_t end = clock(); - const uint64_t end_us = get_time_us(); - printf("%s: elapsed ticks: %ld\n", __func__, end - start); - printf("%s: elapsed us: %llu / %f ms\n", __func__, end_us - start_us, (end_us - start_us) / 1000.0 / nIter); - } - - printf("%f\n", sum); - - free(src0); - free(src1); - free(dst); - - free(src0_fp16); - free(src1_fp16); - - return 0; -} diff --git a/tests/test0.c b/tests/test0.c deleted file mode 100644 index 09154b46..00000000 --- a/tests/test0.c +++ /dev/null @@ -1,42 +0,0 @@ -#include "ggml.h" - -#include -#include - -int main(int argc, const char ** argv) { - struct ggml_init_params params = { - .mem_size = 128*1024*1024, - .mem_buffer = NULL, - .no_alloc = false, - }; - - struct ggml_context * ctx0 = ggml_init(params); - - struct ggml_tensor * t1 = ggml_new_tensor_1d(ctx0, GGML_TYPE_F32, 10); - struct ggml_tensor * t2 = ggml_new_tensor_2d(ctx0, GGML_TYPE_I16, 10, 20); - struct ggml_tensor * t3 = ggml_new_tensor_3d(ctx0, GGML_TYPE_I32, 10, 20, 30); - - GGML_ASSERT(ggml_n_dims(t1) == 1); - GGML_ASSERT(t1->ne[0] == 10); - GGML_ASSERT(t1->nb[1] == 10*sizeof(float)); - - GGML_ASSERT(ggml_n_dims(t2) == 2); - GGML_ASSERT(t2->ne[0] == 10); - GGML_ASSERT(t2->ne[1] == 20); - GGML_ASSERT(t2->nb[1] == 10*sizeof(int16_t)); - GGML_ASSERT(t2->nb[2] == 10*20*sizeof(int16_t)); - - GGML_ASSERT(ggml_n_dims(t3) == 3); - GGML_ASSERT(t3->ne[0] == 10); - GGML_ASSERT(t3->ne[1] == 20); - GGML_ASSERT(t3->ne[2] == 30); - GGML_ASSERT(t3->nb[1] == 10*sizeof(int32_t)); - GGML_ASSERT(t3->nb[2] == 10*20*sizeof(int32_t)); - GGML_ASSERT(t3->nb[3] == 10*20*30*sizeof(int32_t)); - - ggml_print_objects(ctx0); - - ggml_free(ctx0); - - return 0; -} diff --git a/tests/test0.zig b/tests/test0.zig deleted file mode 100644 index 5e625942..00000000 --- a/tests/test0.zig +++ /dev/null @@ -1,41 +0,0 @@ -const std = @import("std"); -const c = @cImport({ - @cInclude("ggml.h"); -}); - -pub fn main() !void { - const params = .{ - .mem_size = 128 * 1024 * 1024, - .mem_buffer = null, - .no_alloc = false, - }; - - const ctx0 = c.ggml_init(params); - defer c.ggml_free(ctx0); - - const t1 = c.ggml_new_tensor_1d(ctx0, c.GGML_TYPE_F32, 10); - const t2 = c.ggml_new_tensor_2d(ctx0, c.GGML_TYPE_I16, 10, 20); - const t3 = c.ggml_new_tensor_3d(ctx0, c.GGML_TYPE_I32, 10, 20, 30); - - try std.testing.expect(c.ggml_n_dims(t1) == 1); - try std.testing.expect(t1.*.ne[0] == 10); - try std.testing.expect(t1.*.nb[1] == 10 * @sizeOf(f32)); - - try std.testing.expect(c.ggml_n_dims(t2) == 2); - try std.testing.expect(t2.*.ne[0] == 10); - try std.testing.expect(t2.*.ne[1] == 20); - try std.testing.expect(t2.*.nb[1] == 10 * @sizeOf(i16)); - try std.testing.expect(t2.*.nb[2] == 10 * 20 * @sizeOf(i16)); - - try std.testing.expect(c.ggml_n_dims(t3) == 3); - try std.testing.expect(t3.*.ne[0] == 10); - try std.testing.expect(t3.*.ne[1] == 20); - try std.testing.expect(t3.*.ne[2] == 30); - try std.testing.expect(t3.*.nb[1] == 10 * @sizeOf(i32)); - try std.testing.expect(t3.*.nb[2] == 10 * 20 * @sizeOf(i32)); - try std.testing.expect(t3.*.nb[3] == 10 * 20 * 30 * @sizeOf(i32)); - - c.ggml_print_objects(ctx0); - - _ = try std.io.getStdIn().reader().readByte(); -}