这是Spencer Nelson上面例子的可用版本。其中一个谜团是输入矩阵按行主序排列,尽管它似乎调用了底层的Fortran例程dgetri。我认为所有底层的Fortran例程都需要按列主序排列,但我对LAPACK并不是很了解,实际上,我正在使用这个例子来帮助我学习它。但是,除了这个谜团之外:
该示例中的输入矩阵是奇异的。LAPACK试图通过在errorHandler中返回3来告诉您。我将该矩阵中的9更改为19,得到了errorHandler为0表示成功的结果,并将其与Mathematica的结果进行了比较。比较也成功了,并确认了示例中的矩阵应按行主序排列。
以下是工作代码:
#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];
int errorHandler;
double lapackWorkspace[9];
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; }
我在 Mac 上按以下方式构建并运行它:
gcc main.c -llapacke -llapack
./a.out
我对LAPACKE库进行了nm
操作,发现以下结果:
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
我看到有一个LAPACKE [sic]包装器,可以方便地为Fortran提供地址,但我可能不会尝试它,因为我有一种前进的方式。
编辑
这是一个可行的版本,绕过LAPACKE [sic],直接使用LAPACK Fortran例程。 我不明白为什么行主输入会产生正确的结果,但我在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];
int errorHandler;
double lapackWorkspace[9];
extern void dgetrf_ (int * m, int * n, double * A, int * LDA, int * IPIV,
int * INFO);
extern void dgetri_ (int * n, double * A, int * LDA, int * IPIV,
double * WORK, int * LWORK, int * INFO);
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; }
像这样构建和运行:
$ 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
编辑
这个谜团似乎不再是谜团了。我认为计算是按列主序进行的,因为必须这样做,但我输入和输出矩阵时却像它们是按行主序排列的。我有两个错误互相抵消,所以事情看起来像是按行排列,即使它们实际上是按列排列的。
delete[] IPIV
和delete [] work
。 - Reb.Cabinnew
调用更改为malloc
,将delete[]
调用更改为frees
(并且摆脱extern "C")。 - alfalfasprout