CODEWITHSUNDEEP
c++matrixfloating-point-precision
Da Zheng
Da Zheng

Reputation: 121

My matrix multiplication doesn't have the same precision as BLAS

I wrote a matrix multiplication routine in C++ and used it in an eigensolver. Then I realized that my implementation has precision problems even though I already use double-precision float-points. It seems the result of my implementation doesn't have enough precision. I tested it with BLAS and found its result doesn't match with BLAS's. The simple code below illustrates the problem and has the following output.

0: 0
1: 0
2: -1.42109e-14
test: test.cpp:59: int main(): Assertion `C1[i] == C2[i]' failed.

I can simply call BLAS routines for matrix multiplication, but I want to know why this happens. I don't know if this is an architecture-related problem. I use Intel Xeon E5-4620.

BTW, I tried to define sum with "long double" and still got different results.

Thanks, Da

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <cblas.h>
#include <assert.h>

#include <vector>

void init_vec(std::vector<double> &vec)
{
    for (size_t i = 0; i < vec.size(); i++)
        vec[i] = ((double) random()) / ((double) random());
}

void mv(const std::vector<double> &A, const std::vector<double> &B,
        std::vector<double> &C)
{
    size_t num_rows = A.size() / B.size();
    size_t num_cols = B.size();
    for (size_t i = 0; i < num_rows; i++) {
        double sum = 0;
        for (size_t j = 0; j < num_cols; j++)
            sum += A[j * num_rows + i] * B[j];
        C[i] = sum;
    }
}

int main()
{
    std::vector<double> A(100 * 10);
    std::vector<double> B(10);
    std::vector<double> C1(100);
    std::vector<double> C2(100);
    init_vec(A);
    init_vec(B);

    mv(A, B, C1);
    cblas_dgemv(CblasColMajor, CblasNoTrans, 100, 10, 1, A.data(), 100, B.data(), 1, 0, C2.data(), 1);
    for (size_t i = 0; i < C1.size(); i++) {
        printf("%ld: %g\n", i, C1[i] - C2[i]);
        assert(C1[i] == C2[i]);
    }
}

Upvotes: 0

Views: 605

Answers (1)

Da Zheng
Da Zheng

Reputation: 121

I realized where the problem is now. BLAS indeed used 80-bit float-point for multiplication. The result of the code below matches with BLAS.

void mv(const std::vector<double> &A, const std::vector<double> &B,
       std::vector<double> &C)
{
    size_t num_rows = A.size() / B.size();
    size_t num_cols = B.size();
    for (size_t i = 0; i < num_rows; i++) {
        double sum = 0;
        for (size_t j = 0; j < num_cols; j++) {
            long double tmp1 = A[j * num_rows + i];
            long double tmp2 = B[j];
            sum += tmp2 * tmp1;
        }
        C[i] = sum;
    }
}

It seems 80-bit float points are only used in multiplication. The summation still uses 64-bit. I hope this can help someone in the future.

Upvotes: 3

Related Questions