Đây là phiên bản hoạt động của ví dụ của Spencer Nelson ở trên. Một bí ẩn về nó là ma trận đầu vào nằm trong thứ tự hàng lớn, mặc dù nó xuất hiện để gọi thường trình fortran cơ bản là dgetri. Tôi tin rằng tất cả các thói quen fortran cơ bản đòi hỏi thứ tự cột lớn, nhưng tôi không có chuyên gia về LAPACK, trên thực tế, tôi đang sử dụng ví dụ này để giúp tôi tìm hiểu nó. Nhưng, một điều bí ẩn đó sang một bên:
Ma trận đầu vào trong ví dụ này là số ít. LAPACK cố gắng cho bạn biết rằng bằng cách trả lại 3 trong số errorHandler. Tôi đã thay đổi 9 trong ma trận đó thành một số 19, nhận được một số errorHandler của 0 báo hiệu thành công và so sánh kết quả với điều đó từ Mathematica. So sánh cũng thành công và xác nhận rằng ma trận trong ví dụ phải theo thứ tự hàng lớn, như đã trình bày.
Dưới đây là đoạn code làm việc:
#include <stdio.h>
#include <stddef.h>
#include <lapacke.h>
int main() {
int N = 3;
int NN = 9;
double M[3][3] = { {1 , 2 , 3},
{4 , 5 , 6},
{7 , 8 , 9} };
int pivotArray[3]; //since our matrix has three rows
int errorHandler;
double lapackWorkspace[9];
// dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix
// called A, sending the pivot indices to IPIV, and spitting error information
// to INFO. also don't forget (like I did) that when you pass a two-dimensional
// array to a function you need to specify the number of "rows"
dgetrf_(&N, &N, M[0], &N, pivotArray, &errorHandler);
printf ("dgetrf eh, %d, should be zero\n", errorHandler);
dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler);
printf ("dgetri eh, %d, should be zero\n", errorHandler);
for (size_t row = 0; row < N; ++row)
{ for (size_t col = 0; col < N; ++col)
{ printf ("%g", M[row][col]);
if (N-1 != col)
{ printf (", "); } }
if (N-1 != row)
{ printf ("\n"); } }
return 0; }
tôi đã xây dựng và chạy nó như sau trên máy Mac:
gcc main.c -llapacke -llapack
./a.out
Tôi đã làm một nm trên thư viện LAPACKE và thấy như sau:
liblapacke.a(lapacke_dgetri.o):
U _LAPACKE_dge_nancheck
0000000000000000 T _LAPACKE_dgetri
U _LAPACKE_dgetri_work
U _LAPACKE_xerbla
U _free
U _malloc
liblapacke.a(lapacke_dgetri_work.o):
U _LAPACKE_dge_trans
0000000000000000 T _LAPACKE_dgetri_work
U _LAPACKE_xerbla
U _dgetri_
U _free
U _malloc
và có vẻ như có một trình bao bọc LAPACKE [sic] có lẽ sẽ giảm bớt chúng tôi phải có địa chỉ ở khắp mọi nơi để thuận tiện cho fortran, nhưng tôi có lẽ sẽ không đi xung quanh để thử nó bởi vì tôi có một con đường phía trước.
EDIT
Đây là một phiên bản làm việc mà bỏ qua LAPACKE [sic], thói quen sử dụng fortran LAPACK trực tiếp. Tôi không hiểu tại sao một đầu vào hàng lớn tạo ra kết quả chính xác, nhưng tôi đã xác nhận nó một lần nữa trong Mathematica.
#include <stdio.h>
#include <stddef.h>
int main() {
int N = 3;
int NN = 9;
double M[3][3] = { {1 , 2 , 3},
{4 , 5 , 6},
{7 , 8 , 19} };
int pivotArray[3]; //since our matrix has three rows
int errorHandler;
double lapackWorkspace[9];
/* from http://www.netlib.no/netlib/lapack/double/dgetrf.f
SUBROUTINE DGETRF(M, N, A, LDA, IPIV, INFO)
*
* -- LAPACK routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
* ..
* .. Array Arguments ..
INTEGER IPIV(*)
DOUBLE PRECISION A(LDA, *)
*/
extern void dgetrf_ (int * m, int * n, double * A, int * LDA, int * IPIV,
int * INFO);
/* from http://www.netlib.no/netlib/lapack/double/dgetri.f
SUBROUTINE DGETRI(N, A, LDA, IPIV, WORK, LWORK, INFO)
*
* -- LAPACK routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, N
* ..
* .. Array Arguments ..
INTEGER IPIV(*)
DOUBLE PRECISION A(LDA, *), WORK(*)
*/
extern void dgetri_ (int * n, double * A, int * LDA, int * IPIV,
double * WORK, int * LWORK, int * INFO);
// dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix
// called A, sending the pivot indices to IPIV, and spitting error information
// to INFO. also don't forget (like I did) that when you pass a two-dimensional
// array to a function you need to specify the number of "rows"
dgetrf_(&N, &N, M[0], &N, pivotArray, &errorHandler);
printf ("dgetrf eh, %d, should be zero\n", errorHandler);
dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler);
printf ("dgetri eh, %d, should be zero\n", errorHandler);
for (size_t row = 0; row < N; ++row)
{ for (size_t col = 0; col < N; ++col)
{ printf ("%g", M[row][col]);
if (N-1 != col)
{ printf (", "); } }
if (N-1 != row)
{ printf ("\n"); } }
return 0; }
xây dựng và chạy như thế này:
$ gcc foo.c -llapack
$ ./a.out
dgetrf eh, 0, should be zero
dgetri eh, 0, should be zero
-1.56667, 0.466667, 0.1
1.13333, 0.0666667, -0.2
0.1, -0.2, 0.1
EDIT
Các bí ẩn không còn dường như là một điều bí ẩn. Tôi nghĩ rằng các tính toán đang được thực hiện theo thứ tự cột lớn, vì chúng phải, nhưng tôi vừa nhập và in các ma trận như thể chúng nằm trong thứ tự lớn. Tôi có hai lỗi mà hủy bỏ lẫn nhau để mọi thứ trông hàng-ish mặc dù họ đang cột-ish.
Về 1x9 hoặc 3x3. Có thực sự không có bất kỳ sự khác biệt về bố trí bộ nhớ. Trong thực tế, các thói quen BLAS/LAPACK không nhận mảng 2d. Họ lấy 1d mảng và đưa ra giả định về cách bạn đã đặt nó ra. Hãy coi chừng mặc dù BLAS và LAPACK sẽ giả định thứ tự FORTRAN (hàng thay đổi nhanh nhất) thay vì đặt hàng C. – MRocklin
Bạn có thể cần 'LAPACKE_dgetrf' thay vì thường trình fortran' dgetrf_'. Ditto 'dgetri'. Khác bạn phải gọi tất cả mọi thứ với địa chỉ. –