diff --git a/2025.09.28_cpp_openblas_test/openblas_test.cpp b/2025.09.28_cpp_openblas_test/openblas_test.cpp new file mode 100644 index 0000000..2fc789e --- /dev/null +++ b/2025.09.28_cpp_openblas_test/openblas_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include + +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 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 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 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; +} \ No newline at end of file diff --git a/2025.09.28_cpp_openblas_test/openblas_test_wtih_eigen.cpp b/2025.09.28_cpp_openblas_test/openblas_test_wtih_eigen.cpp new file mode 100644 index 0000000..a58ef38 --- /dev/null +++ b/2025.09.28_cpp_openblas_test/openblas_test_wtih_eigen.cpp @@ -0,0 +1,19 @@ +#define EIGEN_USE_BLAS +#include +#include "Eigen/Dense" // #include + +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; +} \ No newline at end of file diff --git a/2025.09.28_fortran_openblas_test/openblas_test.f90 b/2025.09.28_fortran_openblas_test/openblas_test.f90 new file mode 100644 index 0000000..0904b4c --- /dev/null +++ b/2025.09.28_fortran_openblas_test/openblas_test.f90 @@ -0,0 +1,25 @@ +program invert_simple + implicit none + integer, parameter :: n = 2 + double precision :: a(n,n) + integer :: ipiv(n) + double precision :: work(n) ! 最小工作空间(LWORK = N) + integer :: info + + ! 初始化一个简单的可逆矩阵: [[2, 1], [1, 2]] + a(1,1) = 2.0d0; a(1,2) = 1.0d0 + a(2,1) = 1.0d0; a(2,2) = 2.0d0 + + ! LU 分解 + call DGETRF(n, n, a, n, ipiv, info) + if (info /= 0) stop 'DGETRF failed' + + ! 求逆(使用最小工作空间 LWORK = N) + call DGETRI(n, a, n, ipiv, work, n, info) + if (info /= 0) stop 'DGETRI failed' + + ! 输出结果(应为 [[0.6667, -0.3333], [-0.3333, 0.6667]]) + print *, 'Inverse:' + print *, a(1,1), a(1,2) + print *, a(2,1), a(2,2) +end program \ No newline at end of file