Вот рабочая версия примера Спенсера Нельсона выше. Одна загадка об этом заключается в том, что входная матрица находится в строчном порядке, хотя она, как представляется, вызывает базовую процедуру fortran dgetri. Я убежден, что все базовые подпрограммы fortran требуют порядка столбцов, но я не эксперт по LAPACK, на самом деле, я использую этот пример, чтобы помочь мне изучить его. Но, эта одна тайна в стороне:
Входная матрица в примере является единственной. LAPACK пытается сказать вам, что, вернув 3 в errorHandler. Я изменил 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
Я сделал nm на библиотеку LAPACKE и обнаружил следующее:
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, но я, вероятно, не собираюсь обходить его, потому что у меня есть путь вперед.
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
EDIT
тайна больше не появляется, чтобы быть тайной. Я думаю, что вычисления выполняются в порядке столбцов, как и должно быть, но я ввожу и печатаю матрицы так, как если бы они были в порядке строк. У меня есть две ошибки, которые отменяют друг друга, поэтому вещи выглядят рядами, даже если они столбцы.
Относительно 1x9 или 3x3. Там нет никакой разницы в плане макета памяти. На самом деле процедуры BLAS/LAPACK не принимают 2d-массивы. Они берут 1d массивы и делают предположения о том, как вы это положили. Остерегайтесь, хотя BLAS и LAPACK будут принимать заказы FORTRAN (быстрее меняются строки), а не заказывать C. – MRocklin
Возможно, вам понадобится 'LAPACKE_dgetrf', а не программа fortran' dgetrf_'. То же 'dgetri'. Иначе вы должны позвонить всем с адресами. –