# 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 $<TARGET_FILE:${TEST_TARGET}> 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 $<TARGET_FILE:${TEST_TARGET}>)
-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 $<TARGET_FILE:${TEST_TARGET}>)
-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
add_test(NAME ${TEST_TARGET} COMMAND $<TARGET_FILE:${TEST_TARGET}>)
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 $<TARGET_FILE:${TEST_TARGET}>)
- set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw")
-
#
# test-pool
add_test(NAME ${TEST_TARGET} COMMAND $<TARGET_FILE:${TEST_TARGET}>)
set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw")
-
#
# test-conv2d
add_test(NAME ${TEST_TARGET} COMMAND $<TARGET_FILE:${TEST_TARGET}>)
set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw")
-
#
# test-conv2d-dw
add_test(NAME ${TEST_TARGET} COMMAND $<TARGET_FILE:${TEST_TARGET}>)
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 $<TARGET_FILE:${TEST_TARGET}>)
- set_property(TEST ${TEST_TARGET} PROPERTY ENVIRONMENT "LLVM_PROFILE_FILE=${TEST_TARGET}.profraw")
-
#
# test-cont
+++ /dev/null
-#include "ggml.h"
-#include "ggml-cpu.h"
-
-#include <stdint.h>
-#include <stdio.h>
-#include <assert.h>
-#include <stdlib.h>
-#include <string.h>
-#include <time.h>
-#include <math.h>
-
-#include <sys/time.h>
-
-#include <arm_neon.h>
-
-#include <Accelerate/Accelerate.h>
-
-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;
-}
+++ /dev/null
-#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 <cassert>
-#include <cmath>
-#include <cstdio>
-#include <cstring>
-#include <fstream>
-#include <map>
-#include <string>
-#include <vector>
-
-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<uint8_t> 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<float> 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<float> 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;
-}
+++ /dev/null
-#define _CRT_SECURE_NO_DEPRECATE // Disables ridiculous "unsafe" warnigns on Windows
-#include "ggml.h"
-#include "ggml-cpu.h"
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <assert.h>
-#include <inttypes.h>
-
-#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;
-}
+++ /dev/null
-#include <stdint.h>
-#include <stdio.h>
-#include <assert.h>
-#include <stdlib.h>
-#include <string.h>
-#include <time.h>
-#include <math.h>
-
-#include <sys/time.h>
-
-#include <arm_neon.h>
-
-#include <Accelerate/Accelerate.h>
-
-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;
-}
+++ /dev/null
-// quantized matrix multiplication
-
-#include "ggml.h"
-
-#include <float.h>
-#include <stdint.h>
-#include <stdio.h>
-#include <inttypes.h>
-#include <assert.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-
-#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 <intrin.h>
-#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;
-}
+++ /dev/null
-// SVD dimensionality reduction
-
-#include <float.h>
-#include <stdint.h>
-#include <stdio.h>
-#include <assert.h>
-#include <stdlib.h>
-#include <string.h>
-#include <time.h>
-#include <math.h>
-
-#include <sys/time.h>
-
-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;
-}
+++ /dev/null
-#include <stdio.h>
-#include <assert.h>
-#include <stdlib.h>
-#include <time.h>
-
-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;
-}
+++ /dev/null
-#include <stdint.h>
-#include <stdio.h>
-#include <assert.h>
-#include <stdlib.h>
-#include <time.h>
-#include <math.h>
-
-#include <sys/time.h>
-
-#include <immintrin.h>
-
-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;
-}
+++ /dev/null
-#include <stdint.h>
-#include <stdio.h>
-#include <assert.h>
-#include <stdlib.h>
-#include <time.h>
-#include <math.h>
-
-#include <sys/time.h>
-
-#include <arm_neon.h>
-
-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;
-}
+++ /dev/null
-#include "ggml.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-
-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;
-}
+++ /dev/null
-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();
-}
overview / pkg / ggml / sources / ggml / commitdiff