下面是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]; //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; }
我建立并运行它在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 [原文如此]包装,可能会缓解我们不得不为了fortran的方便而到处寻址,但我可能不会试着去尝试它,因为我有前进的方向。
EDIT
这里是绕过LAPACKE [原文如此],直接使用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]; //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; }
构建并运行这样的:
$ 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
编辑
谜似乎不再是一个谜。我认为计算是按照列主要的顺序完成的,因为它们必须这样做,但我同时输入和打印这些矩阵,就好像它们是按照主要顺序排列的一样。我有两个互相抵消的错误,所以事情看起来很划算,尽管它们是列式的。
关于1X9或3×3。内存布局方面确实没有任何区别。事实上,BLAS/LAPACK例程不需要2d阵列。他们采用一维数组,并假设你如何布置它。请注意BLAS和LAPACK将假设FORTRAN排序(行更改最快)而不是C排序。 – MRocklin 2012-10-18 18:28:09
您可能需要'LAPACKE_dgetrf'而不是fortran例程'dgetrf_'。同上'dgetri'。否则你必须使用地址来调用所有内容。 – 2016-05-20 02:32:48