Incorrect result with `cblas_dgemv` vs reference netlib and other libraries
メンテナーはふだん 1 日以内に返信
まだ誰も着手していません。
評価
- 難易度
- 4/5
- 見積もり時間
- 3〜5日
- 初心者へのやさしさ
- 28/100
- issue の種類
- バグ
- 明瞭さ
- おおむね明確
- 活発さ
- 停滞
- 技術スタック
- cpp
- 領域
- performance
調査の方向性
まず、issue に記載されているコマンドを使って、添付された C++ の再現コードを OpenBLAS および netlib CBLAS でコンパイルします。次に、cblas_dgemv と cblas_ddot の結果を単純な乗算と比較して調査します。入力には reproduction.txt を使用し、報告されている環境間で動作を比較します。最後の要素の不一致の原因を特定し、示されている数値精度を失わずに修正後の結果を実証できれば完了です。
索引モデルが issue の本文から書いたものです。
説明
We recently switched to testing openBLAS on a project and are noticing some test case failures due to a matrix multiplication operation returning an incorrect result.
This issue has been observed on a variety of platforms (Ubuntu 22.04, RHEL7, RHEL9, MSYS2 mingw), a variety of compilers (clang-15, mingw-13, gcc-12, gcc-11, gcc-9), as well as a variety of openblas versions(0.3.3, 0.3.20, 0.3.21, 0.3.24), and a variety of CPUs:
- Intel(R) Xeon(R) CPU E5-4650 0 @ 2.70GHz
- Intel(R) Xeon(R) Gold 6226R CPU @ 2.90GHz
- Intel(R) Core(TM) i7-9850H CPU @ 2.60GHz 2.59 GHz
Reproduction
I have attached a minimally reproducible example (in C++) showing the problem
Reproduction Code
#include <cblas.h>
#include <cassert>
#include <iostream>
#include <fstream>
constexpr static size_t size = 16;
void get_values(double* A, double* b, const char* binfile) {
std::fstream binaryReader;
binaryReader.open(binfile, std::ios::in | std::ios::binary);
assert(binaryReader.is_open());
for (size_t i = 0; i < (size * size); i++) {
binaryReader.read(reinterpret_cast<char *>(&A[i]), sizeof(double));
}
for (size_t i = 0; i < size; i++) {
binaryReader.read(reinterpret_cast<char *>(&b[i]), sizeof(double));
}
}
void matrixMultiply(const double* A, const double* b, double* c) {
const size_t N = size;
for (size_t i = 0; i < N; ++i) {
for (size_t j = 0; j < N; ++j){
c[i] += A[ (i * N) + j] * b[j];
}
}
}
void CBlasMatrixMultiply(const double* a, const double* b, double* c){
int32_t M = size;
int32_t N = size;
cblas_dgemv(CblasRowMajor,
CblasNoTrans,
M, N,
1.0, a, N,
b, 1,
1.0, c, 1);
}
int main(int argc, char* argv[]) {
auto* A = new double[256]();
auto* b = new double[16]();
const char* binfile = argc > 1 ? argv[1] : BIN_FILE;
get_values(A, b, binfile);
auto* c_blas = new double[16]();
CBlasMatrixMultiply(A, b, c_blas);
auto *c_mat = new double[16]();
matrixMultiply(A, b, c_mat);
printf(" BLAS MAT\n");
for (size_t i = 0; i < size; i++) {
printf("%0.18e\t%0.18e\n", c_blas[i], c_mat[i]);
}
// dot product?
auto c_ddot = cblas_ddot(size, A + (size * 15), 1, b, 1);
printf("ddot: %0.18e\n", c_ddot);
printf("dgemv: %0.18e\n", c_blas[15]);
printf("ddot == dgemv? %s\n", c_ddot == c_blas[15] ? "YES" : "NO");
printf("dgemv is 0.0? %s\n", c_blas[15] == 0.0 ? "YES" : "NO");
delete[] A;
delete[] b;
delete[] c_blas;
delete[] c_mat;
return 0;
}
Compile this code with:
g++ -o blas_test blas_test.cpp -DBIN_FILE=\"path/to/bin\" $(pkg-config --libs --cflags openblas)
And observe the following output:
OpenBLAS result
BLAS MAT
1.175201193643801378e+00 1.175201193643801822e+00
1.103638323514327002e+00 1.103638323514327224e+00
3.578143506473725477e-01 3.578143506473724922e-01
7.045563366848892062e-02 7.045563366848883735e-02
9.965128148869371871e-03 9.965128148869309421e-03
1.099586127207577424e-03 1.099586127207556390e-03
9.945433911373591229e-05 9.945433911360671605e-05
7.620541308983597162e-06 7.620541308896932986e-06
5.064719745540013918e-07 5.064719744437437483e-07
2.971814122565419325e-08 2.971814140421127963e-08
1.560886642160141946e-09 1.560886564391797127e-09
7.419920233786569952e-11 7.419902813631537848e-11
3.221201083647429186e-12 3.221406076314455686e-12
1.292299600663682213e-13 1.289813102624490508e-13
4.440892098500626162e-15 4.440892098500626162e-15
0.000000000000000000e+00 -4.440892098500626162e-16
ddot: 0.000000000000000000e+00
dgemv: 0.000000000000000000e+00
ddot == dgemv? YES
dgemv is 0.0? YES
It is worth noting that only the last value is different outside of acceptable numerical precision, and that every other value passes within 1e-16. Furthermore, a value of exactly 0.0 is, in itself, suspicious, as there's no real circumstance the value could be that.
Change the compile command to:
# cblas here is netlib
g++ -o blas_test blas_test.cpp -DBIN_FILE=\"path/to/bin\" $(pkg-config --libs --cflags cblas)
and observe this result:
netlib result
BLAS MAT
1.175201193643801822e+00 1.175201193643801822e+00
1.103638323514327224e+00 1.103638323514327224e+00
3.578143506473724922e-01 3.578143506473724922e-01
7.045563366848883735e-02 7.045563366848883735e-02
9.965128148869309421e-03 9.965128148869309421e-03
1.099586127207556390e-03 1.099586127207556390e-03
9.945433911360671605e-05 9.945433911360671605e-05
7.620541308896932986e-06 7.620541308896932986e-06
5.064719744437437483e-07 5.064719744437437483e-07
2.971814140421127963e-08 2.971814140421127963e-08
1.560886564391797127e-09 1.560886564391797127e-09
7.419902813631537848e-11 7.419902813631537848e-11
3.221406076314455686e-12 3.221406076314455686e-12
1.289813102624490508e-13 1.289813102624490508e-13
4.440892098500626162e-15 4.440892098500626162e-15
-4.440892098500626162e-16 -4.440892098500626162e-16
ddot: -4.440892098500626162e-16
dgemv: -4.440892098500626162e-16
ddot == dgemv? YES
dgemv is 0.0? NO
Here is the binary file that contains a 16x16 matrix and a 16x1 vector:
(NOTE: This is a binary data file, extension changed to make github happy)
reproduction.txt
Other Notes
We have done extensive testing in other BLAS-like environments to get a result close to the expected -4e-16 result, which passes our test. Both MATLAB (2023a) and numpy (1.26 w/ MKL) return a result very close to what we expect, and pass our test. And, obviously, our naive matrix multiplication in the reproduction code gives
The matrix in question is not overly ill-conditioned, it has a condition number of ~10.
- 主要言語
- C
- スター
- 7.6k
- フォーク
- 1.7k
- 平均マージ
- 1日 6時間
- マージ済み PR(30日)
- 46
環境構築
このプロジェクトには開発コンテナ、Dockerfile、コントリビューションガイドがありません。まず README を読み、一般的な手順ははじめてのコントリビューションガイドを参照してください。
はじめの一歩
- issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
- 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
- リポジトリをフォークし、ブランチを切って変更します。
- issue 番号を参照したプルリクエストを送ります。
OpenMathLib/OpenBLAS のほかの issue
-
難易度 1/5 1時間未満 初心者へのやさしさ 88/100
OpenMathLib/OpenBLAS#6062 · コメント 2 件 ·
メンテナーはふだん 1 日以内に返信
-
難易度 4/5 3〜5日 初心者へのやさしさ 52/100
OpenMathLib/OpenBLAS#6059 ·
メンテナーはふだん 1 日以内に返信
-
難易度 4/5 3〜5日 初心者へのやさしさ 48/100
OpenMathLib/OpenBLAS#6029 · コメント 21 件 ·
メンテナーはふだん 1 日以内に返信
-
難易度 3/5 1〜2日 初心者へのやさしさ 68/100
OpenMathLib/OpenBLAS#6028 · コメント 1 件 ·
メンテナーはふだん 1 日以内に返信
-
難易度 4/5 3〜5日 初心者へのやさしさ 35/100
OpenMathLib/OpenBLAS#6005 · コメント 21 件 · リアクション 2 件 ·
メンテナーはふだん 1 日以内に返信
OpenMathLib/OpenBLAS の issue をすべて見る
似ている issue
-
難易度 2/5 1〜3時間 初心者へのやさしさ 88/100
ARM-software/sysarch-acs#556 · コメント 1 件 ·
-
難易度 2/5 1〜3時間 初心者へのやさしさ 68/100
メンテナーはふだん 1 日以内に返信
-
難易度 2/5 1〜3時間 初心者へのやさしさ 78/100
メンテナーはふだん 1 日以内に返信
-
bug needs triage
難易度 2/5 1〜3時間 初心者へのやさしさ 78/100
netdata/netdata#24062 · コメント 1 件 ·
メンテナーはふだん 1 日以内に返信
-
難易度 2/5 1〜3時間 初心者へのやさしさ 76/100