这是之前相关的几篇:
本篇介绍 C++ 代码实现 OpenBLAS 的链接和编译。
一、C++ 直接调用 OpenBLAS
C++ 代码调用 OpenBLAS 矩阵求逆的例子 openblas_test.cpp:
#include <iostream>
#include <vector>
#include <cmath>
extern "C" {
void dgetrf_(int* m, int* n, double* a, int* lda, int* ipiv, int* info);
void dgetri_(int* n, double* a, int* lda, int* ipiv, double* work, int* lwork, int* info);
}
int main() {
const int N = 3;
double A_data[] = {
1, 0, 5, // 第1列
2, 1, 6, // 第2列
3, 4, 0 // 第3列
};
std::vector<double> A(A_data, A_data + 9);
std::cout << "Original matrix (column-major):\n";
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
std::cout << A[j*N + i] << " ";
}
std::cout << "\n";
}
// LU 分解
std::vector<int> ipiv(N);
int info;
int n = N;
int lda = N;
dgetrf_(&n, &n, &A[0], &lda, &ipiv[0], &info);
if (info != 0) {
std::cerr << "LU decomposition failed! (info=" << info << ")\n";
return 1;
}
// 计算逆矩阵
std::vector<double> work(N);
int lwork = N; // 关键:必须是可变变量
dgetri_(&n, &A[0], &lda, &ipiv[0], &work[0], &lwork, &info);
if (info != 0) {
std::cerr << "Matrix inversion failed! (info=" << info << ")\n";
return 1;
}
std::cout << "\nInverse matrix:\n";
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
std::cout << A[j*N + i] << " ";
}
std::cout << "\n";
}
return 0;
}
编译命令和 Fortran 的命令差不多,参考:Fortran程序的OpenBLAS链接和编译。
默认 OpenBLAS 环境,动态链接编译:
g++ -o test openblas_test.cpp -lopenblas
指定 OpenBLAS 路径,静态链接 OpenBLAS,动态链接 pthread 库,编译:
g++ -o test openblas_test.cpp ~/OpenBLAS/lib/libopenblas.a -lpthread
指定 OpenBLAS 路径,完全静态链接编译:
g++ -static -o test openblas_test.cpp ~/OpenBLAS/lib/libopenblas.a -lpthread
可以查看编译后的可执行文件 test 的大小,静态链接编译的文件会更大些:
32K test_1
316K test_2
2.8M test_3
编译后运行:
./test
运行结果:
Original matrix (column-major):
1 2 3
0 1 4
5 6 0
Inverse matrix:
-24 18 5
20 -15 -4
-5 4 1
二、通过 Eigen 调用 OpenBLAS(推荐)
这里通过简单修改代码,通过 Eigen 实现调用 OpenBLAS,而不是使用 Eigen 自带的默认 BLAS。Eigen 环境的安装以及原代码示例参考这篇:C++线性代数库Eigen的下载和代码示例。
代码例子 openblas_test_with_eigen.cpp:
#define EIGEN_USE_BLAS
#include <iostream>
#include "Eigen/Dense" // #include <Eigen/Dense>
int main() {
Eigen::MatrixXd A(4, 4);
A << 2, 1, 1, 0,
1, 3, 1, 1,
1, 1, 4, 1,
0, 1, 1, 5;
Eigen::MatrixXd A_inv = A.inverse();
std::cout << "A:\n" << A << std::endl;
std::cout << "\nA_inv:\n" << A_inv << std::endl;
std::cout << "\nA*A_inv:\n" << A * A_inv << std::endl;
return 0;
}
链接编译(只选一个):
g++ -o test openblas_test_with_eigen.cpp -lopenblas
g++ -o test openblas_test_with_eigen.cpp ~/OpenBLAS/lib/libopenblas.a -lpthread
g++ -static -o test openblas_test_with_eigen.cpp ~/OpenBLAS/lib/libopenblas.a -lpthread
如果编译出现报错,可能是因为编译器默认使用的是 C++98 标准,而软件包要求 C++11 或更高标准。这里考虑指定使用 C++14 标准来编译 C++ 代码,一般会解决问题(只选一个):
g++ -std=c++14 -o test openblas_test_with_eigen.cpp -lopenblas
g++ -std=c++14 -o test openblas_test_with_eigen.cpp ~/OpenBLAS/lib/libopenblas.a -lpthread
g++ -std=c++14 -static -o test openblas_test_with_eigen.cpp ~/OpenBLAS/lib/libopenblas.a -lpthread
编译后运行:
./test
运行结果:
A:
2 1 1 0
1 3 1 1
1 1 4 1
0 1 1 5
A_inv:
0.666667 -0.2 -0.133333 0.0666667
-0.2 0.44 -0.04 -0.08
-0.133333 -0.04 0.306667 -0.0533333
0.0666667 -0.08 -0.0533333 0.226667
A*A_inv:
1 -2.08167e-17 0 6.93889e-18
-1.38778e-16 1 0 -2.77556e-17
-1.38778e-16 0 1 0
0 0 -5.55112e-17 1
【说明:本站主要是个人的一些笔记和代码分享,内容可能会不定期修改。为了使全网显示的始终是最新版本,这里的文章未经同意请勿转载。引用请注明出处:https://www.guanjihuan.com】