ggml_mul_mat result is different with numpy matmul when use same data

[knightXun]

I dump the numpy matmul value into ggml_tensors and use the ggml_tensor's ggml_mul_mat. And I compare the result with numpy, they are totally different. below are my code, can anyone help me to find the reason?

import numpy as np

op1 = np.random.randint(
    low=0,        
    high=128,      
    size=(2048, 20480), 
    dtype=np.int8
)

op2 = np.random.rand(2048,2).astype(np.float32)

with open("op1.bin", "wb") as f:
    f.write(np.int32(op1.shape[0]).tobytes())
    f.write(np.int32(op1.shape[1]).tobytes())
    f.write(op1.tobytes())

with open("op2.bin", "wb") as f:
    f.write(np.int32(op2.shape[0]).tobytes())
    f.write(np.int32(op2.shape[1]).tobytes())
    f.write(op2.tobytes())

result = np.matmul(op1.T, op2)
result = result.astype(np.float32)
# result = result.T

print("op1 shape:", op1.shape)
for i in range(20):
    print(f"op1[{i}][0] = {op1[i][0]},, op1[{i}][1] = {op1[i][1]}")

# op2 = op2.T
print("op2 shape:", op2.shape)
for i in range(20):
    print(f"op2[{i}][0] = {op2[i][0]}, op2[{i}][1] = {op2[i][1]}")


print("Result shape:", result.shape)
for i in range(20):
    print(f"result[{i}][0] = {result[i][0]}, result[{i}][1] = {result[i][1]}")

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stdint.h>
#include <string.h>
#include <unistd.h> 

#include "ggml.h"
#include "ggml-cpu.h"
#include "ggml-backend.h"

void free_f32_array(float** array, int rows) {
    if (array) {
        for (int i = 0; i < rows; i++) {
            free(array[i]);
        }
        free(array);
    }
}

float**read_f32_array_from_file(const char* filename, int* rows, int* cols) {
    char cwd[1024];
    if (getcwd(cwd, sizeof(cwd)) != NULL) {
        printf("current: %s\n", cwd);
    } else {
        perror("fault to getcwd");
    }
    
    FILE* file = fopen(filename, "rb");
    if (!file) {
        perror("can't open file"); 
        printf("try to open file: %s\n", filename);
        return NULL;
    }
    
    if (fread(rows, sizeof(int), 1, file) != 1) {
        printf("read line\n");
        fclose(file);
        return NULL;
    }
    
    if (fread(cols, sizeof(int), 1, file) != 1) {
        printf("read line failed\n");
        fclose(file);
        return NULL;
    }
    
    float** array = (float**)malloc(*rows * sizeof(float*));
    if (!array) {
        printf("malloc failed\n");
        fclose(file);
        return NULL;
    }
    
    for (int i = 0; i < *rows; i++) {
        array[i] = (float*)malloc(*cols * sizeof(float));
        if (!array[i]) {
            printf("malloc failed\n");
            for (int j = 0; j < i; j++) {
                free(array[j]);
            }
            free(array);
            fclose(file);
            return NULL;
        }
    }
    
    for (int i = 0; i < *rows; i++) {
        if (fread(array[i], sizeof(float), *cols, file) != *cols) {
            printf("read failed: %d\n", i);
            free_f32_array(array, *rows);
            fclose(file);
            return NULL;
        }
    }
    
    fclose(file);
    return array;
}

int8_t**read_q8_array_from_file(const char* filename, int* rows, int* cols) {
    char cwd[1024];
    if (getcwd(cwd, sizeof(cwd)) != NULL) {
        printf("current dir: %s\n", cwd);
    } else {
        perror("get current dir failed");
    }
    
    FILE* file = fopen(filename, "rb");
    if (!file) {
        perror("cant open file"); 
        printf("try to open file: %s\n", filename);
        return NULL;
    }
    
    fread(rows, sizeof(int), 1, file);
    fread(cols, sizeof(int), 1, file);
    
    int8_t** array = (int8_t**)malloc(*rows * sizeof(int8_t*));
    if (!array) {
        printf("malloc failed\n");
        fclose(file);
        return NULL;
    }
    
    for (int i = 0; i < *rows; i++) {
        array[i] = (int8_t*)malloc(*cols * sizeof(int8_t));
        if (!array[i]) {
            printf("malloc failed\n");
            for (int j = 0; j < i; j++) {
                free(array[j]);
            }
            free(array);
            fclose(file);
            return NULL;
        }
    }
    
    for (int i = 0; i < *rows; i++) {
        fread(array[i], sizeof(int8_t), *cols, file);
    }
    
    fclose(file);
    return array;
}


 void free_q8_array(int8_t** array, int rows) {
    if (array) {
        for (int i = 0; i < rows; i++) {
            free(array[i]);
        }
        free(array);
    }
}

int test_origin_mat_mul(int nthread, struct ggml_tensor *out,struct ggml_context * ctx)
{

    int n_threads = nthread;
    int n_rounds = 1;

    // Create graph
    struct ggml_cgraph *gf = ggml_new_graph(ctx);

    fprintf(stdout, "before ggml_build_forward_expand\n");

    ggml_build_forward_expand(gf, out);
    int n_nodes = ggml_graph_n_nodes(gf);

    fprintf(stdout, "before Create threadpool\n");
    // Create threadpool
    struct ggml_threadpool_params tpp = ggml_threadpool_params_default(n_threads);
    struct ggml_threadpool *threadpool = ggml_threadpool_new(&tpp);
    if (!threadpool)
    {
        fprintf(stderr, "threadpool create failed : n_threads %d\n", n_threads);
        exit(1);
    }

    // Create compute plan
    struct ggml_cplan cplan = ggml_graph_plan(gf, n_threads, threadpool);

    fprintf(stdout, "origin matmul work_size: %zu\n", cplan.work_size);

    cplan.work_data = (uint8_t *)ggml_new_buffer(ctx, cplan.work_size);

    fprintf(stdout, "graph-compute with n_threads: %d, n_nodes: %d, n_rounds: %d\n", n_threads, n_nodes, n_rounds);
    // ggml_graph_print(gf);

    // Warmup
    ggml_graph_compute(gf, &cplan);

    struct timespec t0, t1;

    int ret = clock_gettime(CLOCK_MONOTONIC, &t0);
    if (ret != 0) {
        fprintf(stderr, "clock_gettime failed for t0\n");
        return -1;
    }


    for (int i = 0; i < n_rounds; i++)
    {
        ggml_graph_compute(gf, &cplan);
    }

    // 获取结束时间
    ret = clock_gettime(CLOCK_MONOTONIC, &t1);
    if (ret != 0) {
        fprintf(stderr, "clock_gettime failed for t1\n");
        return -1;
    }

    int64_t nsec_total = (int64_t)(t1.tv_sec - t0.tv_sec) * 1000000000LL + 
                         (int64_t)(t1.tv_nsec - t0.tv_nsec);
    int64_t usec_total = nsec_total / 1000;  // 转换为微秒

    float usec_per_iter = (float)usec_total / n_rounds;
    float nsec_per_node = (float)nsec_total / (n_rounds * n_nodes);

    fprintf(stderr, "graph-compute took %lld usec \n", (long long)usec_total);
    fprintf(stderr, " %f usec per-iter\n", usec_per_iter);
    fprintf(stderr, " %f nsec per-node\n", nsec_per_node);
    ggml_threadpool_free(threadpool);
    // ggml_free(ctx);

    return 1;
}


int main(int argc, char *argv[])
{
    size_t mem_size = 1024 * 1024 * 1024;  // 1GB
    struct ggml_init_params params1 = {
        /* .mem_size   = */ mem_size,
        /* .mem_buffer = */ NULL,
        /* .no_alloc   = */ false,
    };

    struct ggml_context * ctx1 = ggml_init(params1);


    fprintf(stdout, "-----------------------------------------------\n");
    struct ggml_tensor *a = ggml_new_tensor_4d(
        ctx1, 
        GGML_TYPE_Q8_0, 
        2048, 20480, 1, 1
    );

    struct ggml_tensor *b = ggml_new_tensor_4d(ctx1, GGML_TYPE_F32, 2048, 2, 1, 1);
    struct ggml_tensor *op_tensor1 = ggml_mul_mat(ctx1, a, b);

    fprintf(stdout, "--------------------------------------------\n");
    
    bool src1_cont = ggml_is_contiguous(op_tensor1->src[0]);

    int rows;
    int cols = 0;
    int8_t** op1 = read_q8_array_from_file("/home/firefly/RongKit/tools/matmul/op1.bin", &rows, &cols);

    rows;
    cols = 0;
    float** op2 = read_f32_array_from_file("/home/firefly/RongKit/tools/matmul/op2.bin", &rows, &cols);


    for (int i = 0; i < 2048; i++)
    {
        for( int j = 0; j < 20480; j++)
        {
            int8_t * aptr = NULL;
            int8_t * aaptr = NULL;

            *(int8_t *) ((char *) op_tensor1->src[0]->data + i*op_tensor1->src[0]->nb[1] + j*op_tensor1->src[0]->nb[0]) = op1[i][j];

        }
    }
    fprintf(stdout, "--------------------------------------------\n");
    src1_cont = ggml_is_contiguous(op_tensor1->src[1]);

    for (int i = 0; i < 2048; i++)
    {
        for( int j = 0; j < 2; j++)
        {
            float *aptr = NULL;
            float *aaptr = NULL;
            *(float *) ((float *) op_tensor1->src[1]->data + i*op_tensor1->src[1]->ne[1] + j ) = op2[i][j];
        }
    }
    fprintf(stdout, "--------------------------------------------\n");

    src1_cont = ggml_is_contiguous(op_tensor1->src[0]);
    for (size_t i = 0; i < 20; i++)
    {
        for (size_t j = 0; j < 2; j++)
        {
            int8_t * aptr = NULL;

            aptr = (int8_t *)op_tensor1->src[0]->data +  i * op_tensor1->src[0]->nb[1] + j * op_tensor1->src[0]->nb[0];

            fprintf(stdout, "op_tensor1->src0[%zu][%zu] = %d\n", i, j, *aptr);
        }
    }
    fprintf(stdout, "--------------------------------------------\n");
    for (size_t i = 0; i < 20; i++)
    {
        for (size_t j = 0; j < 2; j++)
        {
            float *aptr = NULL;

            aptr = (float *) op_tensor1->src[1]->data + i*op_tensor1->src[1]->ne[1] + j;

            fprintf(stdout, "op_tensor1->src1[%zu][%zu] = %f\n", i, j, *aptr);
        }
    }


    test_origin_mat_mul(1, op_tensor1, ctx1);


    fprintf(stdout, "----------------------------- compare result ------------------------------------------\n");
    fprintf(stdout, "op_tensor result shape: %d,%d\n", op_tensor1->ne[0], op_tensor1->ne[1]);
    src1_cont = ggml_is_contiguous(op_tensor1);
    for (size_t i = 0; i < 20; i++)
    {
        for (size_t j = 0; j < 2; j++)
        {
            float *aptr = NULL;
            if (src1_cont)
            {
                aptr = (float *)op_tensor1->data +  i * op_tensor1->ne[1] + j;
            } else{
                aptr = (float *)op_tensor1->data +  i * op_tensor1->nb[1] + j * op_tensor1->nb[0];
            }

            fprintf(stdout, "op_tensor1[%zu][%zu] = %f\n", i, j, *aptr);
        }
    }

    fprintf(stdout, "----------------------------- end ------------------------------------------\n");
    fprintf(stdout, "op_tensor1 shape: [%d,%d,%d,%d]\n", op_tensor1->ne[0], op_tensor1->ne[1], op_tensor1->ne[2], op_tensor1->ne[3]);
    return 0;
}

[JohannesGaessler]

The data type GGML_TYPE_Q8_0 is not plain int8 values, it's a struct. The value at index i is effectively d * qs[i]

[knightXun]

@JohannesGaessler Thanks for this great hit! but I have a further question, I used the gguf writer to store f32 value. but the result is still different. Can you help me take a look ?

import numpy as np
import gguf

op1 = np.random.rand(2048, 6144).astype(np.float32)

op2 = np.random.rand(2048,2).astype(np.float32)

result = np.matmul(op1.T, op2)
result = result.astype(np.float32)
# result = result.T

writer = gguf.GGUFWriter("op1.gguf", "arm")
writer.add_string("description", "NumPy data converted to GGUF")
writer.add_uint32("version", 1)
writer.add_tensor("op1", op1.T)
writer.write_header_to_file()
writer.write_kv_data_to_file()
writer.write_tensors_to_file()
writer.close()

writer = gguf.GGUFWriter("op2.gguf", "arm")
writer.add_string("description", "NumPy data converted to GGUF")
writer.add_uint32("version", 1)
writer.add_tensor("op2", op2.T)
writer.write_header_to_file()
writer.write_kv_data_to_file()
writer.write_tensors_to_file()
writer.close()

print("op1 shape:", op1.shape)
for i in range(20):
    print(f"op1[{i}][0] = {op1[i][0]}, op1[{i}][1] = {op1[i][1]}")
print("-----------------------------------------------------------")

# op2 = op2.T
print("op2 shape:", op2.shape)
for i in range(20):
    print(f"op2[{i}][0] = {op2.T[0][i]}, op2[{i}][1] = {op2.T[1][i]}")


print("Result shape:", result.shape)
for i in range(20):
    print(f"result[{i}][0] = {result[i][0]}, result[{i}][1] = {result[i][1]}")

numpy result:

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stdint.h>
#include <string.h>
#include <unistd.h> 

#include "ggml.h"
#include "gguf.h"
#include "ggml-cpu.h"
#include "ggml-backend.h"


int test_origin_mat_mul(int nthread, struct ggml_tensor *out,struct ggml_context * ctx)
{

    int n_threads = nthread;
    int n_rounds = 1;

    // Create graph
    struct ggml_cgraph *gf = ggml_new_graph(ctx);

    fprintf(stdout, "before ggml_build_forward_expand\n");

    ggml_build_forward_expand(gf, out);
    int n_nodes = ggml_graph_n_nodes(gf);

    fprintf(stdout, "before Create threadpool\n");
    // Create threadpool
    struct ggml_threadpool_params tpp = ggml_threadpool_params_default(n_threads);
    struct ggml_threadpool *threadpool = ggml_threadpool_new(&tpp);
    if (!threadpool)
    {
        fprintf(stderr, "threadpool create failed : n_threads %d\n", n_threads);
        exit(1);
    }

    // Create compute plan
    struct ggml_cplan cplan = ggml_graph_plan(gf, n_threads, threadpool);

    fprintf(stdout, "origin matmul work_size: %zu\n", cplan.work_size);

    cplan.work_data = (uint8_t *)ggml_new_buffer(ctx, cplan.work_size);

    fprintf(stdout, "graph-compute with n_threads: %d, n_nodes: %d, n_rounds: %d\n", n_threads, n_nodes, n_rounds);
    // ggml_graph_print(gf);

    // Warmup
    ggml_graph_compute(gf, &cplan);

    struct timespec t0, t1;

    int ret = clock_gettime(CLOCK_MONOTONIC, &t0);
    if (ret != 0) {
        fprintf(stderr, "clock_gettime failed for t0\n");
        return -1;
    }


    for (int i = 0; i < n_rounds; i++)
    {
        ggml_graph_compute(gf, &cplan);
    }

    // 获取结束时间
    ret = clock_gettime(CLOCK_MONOTONIC, &t1);
    if (ret != 0) {
        fprintf(stderr, "clock_gettime failed for t1\n");
        return -1;
    }

    int64_t nsec_total = (int64_t)(t1.tv_sec - t0.tv_sec) * 1000000000LL + 
                         (int64_t)(t1.tv_nsec - t0.tv_nsec);
    int64_t usec_total = nsec_total / 1000;  // 转换为微秒

    float usec_per_iter = (float)usec_total / n_rounds;
    float nsec_per_node = (float)nsec_total / (n_rounds * n_nodes);

    fprintf(stderr, "graph-compute took %lld usec \n", (long long)usec_total);
    fprintf(stderr, " %f usec per-iter\n", usec_per_iter);
    fprintf(stderr, " %f nsec per-node\n", nsec_per_node);
    ggml_threadpool_free(threadpool);
    // ggml_free(ctx);

    return 1;
}


int main(int argc, char *argv[])
{

    size_t mem_size = 1024 * 1024 * 512;

    struct ggml_init_params params1 = {
        /* .mem_size   = */ mem_size,
        /* .mem_buffer = */ NULL,
        /* .no_alloc   = */ false,
    };

    struct ggml_context * ctx1 = ggml_init(params1);

    struct ggml_context * ctx2 = NULL;
    struct gguf_init_params gguf_params1 = {
        /*.no_alloc = */ false,
        /*.ctx      = */ &ctx2,
    };

    gguf_init_from_file("/home/firefly/RongKit/tools/matmul/op1.gguf", gguf_params1);

    if (!ctx2) {
        fprintf(stderr, "Failed to initialize gguf context from file\n");
        return -1;
    }
    
    struct ggml_tensor *op1 = ggml_get_tensor(ctx2, "op1");
    fprintf(stdout, "-----------------------------------------------\n");
    for (size_t i = 0; i < 20; i++)
    {
        for (size_t j = 0; j < 2; j++)
        {
            float val = ggml_get_f32_nd(op1, i, j, 0, 0);
            fprintf(stdout, "op1[%zu][%zu] = %f\n", i, j, val);
        }
    }
    fprintf(stdout, "--------------------------------------------\n");
    struct ggml_context * ctx3 = NULL;
    struct gguf_init_params gguf_params2 = {
        /*.no_alloc = */ false,
        /*.ctx      = */ &ctx3,
    };
    gguf_init_from_file("/home/firefly/RongKit/tools/matmul/op2.gguf", gguf_params2);

    struct ggml_tensor *op2 = ggml_get_tensor(ctx3, "op2");
    fprintf(stdout, "-----------------------------------------------\n");
    for (size_t i = 0; i < 20; i++)
    {
        for (size_t j = 0; j < 2; j++)
        {
            float val = ggml_get_f32_nd(op2, i, j, 0, 0);
            fprintf(stdout, "op2[%zu][%zu] = %f\n", i, j, val);
        }
    }
    fprintf(stdout, "--------------------------------------------\n");

    struct ggml_tensor *op_tensor1 = ggml_mul_mat(ctx1, op1, op2);


    fprintf(stdout, "-----------------------------------------------\n");
    
    test_origin_mat_mul(1, op_tensor1, ctx1);


    fprintf(stdout, "----------------------------- compare result ------------------------------------------\n");
    fprintf(stdout, "op_tensor result shape: %d,%d; type is F32 = %d\n", op_tensor1->ne[0], op_tensor1->ne[1], op_tensor1->type == GGML_TYPE_F32 );


    for (size_t i = 0; i < 20; i++)
    {
        for (size_t j = 0; j < 2; j++)
        {
            float val = ggml_get_f32_nd(op_tensor1, i, j, 0, 0);

            fprintf(stdout, "op_tensor1[%zu][%zu] = %f\n", i, j, val);
        }
    }

    fprintf(stdout, "----------------------------- end ------------------------------------------\n");
    return 0;
}

ggml result:

[knightXun]

@JohannesGaessler I found the reason, it's my mistake, thanks