интерфейс — вызов подпрограммы FORTRAN из c ++ дает недопустимое значение параметра

В настоящее время я пытаюсь решить конкретную проблему собственных значений (так называемая гироскопическая проблема собственных значений).
с большими разреженными матрицами (из сикретизации FEM). Язык программирования — C ++.

Стандартным эталоном для EVP является ARPACK. Увы, он реализует только «классические» процессы Арнольди,
который не подходит для таких проблем (ср. Методы сохранения структуры).

Недавно я нашел это Алгоритм 961 ссылка, которая также предоставляет некоторый код — на Фортране!
Поэтому я попытался включить подпрограмму DGHUTR в C ++, но безрезультатно.
Ниже приведен MWE, который является адаптацией теста для DGHUTR (TDGHUTR.f) в C ++:

#include <Eigen/Dense>
#include <Eigen/Sparse>
//definition stolen from ARPACK++
#define F77NAME(x) x ## _

//Interface to the SHEIG library function DGHUTR
#ifdef __cplusplus
extern "C"{
#endif
void F77NAME(dghutr)( char* JOB, char* COMPQ1, char* COMPQ2, int* N, double* A, int* LDA,
double* DE, int* LDDE, double* C1, int* LDC1, double* VW, int* LDVW,
double* Q1, int* LDQ1, double* Q2, int* LDQ2, double* B, int* LDB,
double* F, int* LDF, double* C2, int* LDC2, double* ALPHAR, double* ALPHAI,
double* BETA, int* IWORK, int* LIWORK, double* DWORK,int* LDWORK, int* INFO );
#ifdef __cplusplus
}
#endifint main(void){
// define system sizes
int N(8),  M(N/2);
std::cout << "Sizes: " << N << '\t' << M << std::endl;char job('E'),  compq1('I'),  compq2('I');
int lda(M),  ldde(M),  ldq1(N),  ldq2(N),  ldb(M),  ldc1(M),  ldc2(M),  ldf(M),  ldvw(M);

int ldwork = 2*N*N+std::max(4*N+4, 32);
int liwork = N+12;// workspace arrays
int* iwork = new int[liwork];
double* dwork = new double[ldwork];
int info(0);
// auxiliary matrices and  vectors
Eigen::MatrixXd F(ldf, M),  C2(ldc2, M),  Q1(ldq1, N),  Q2(ldq2, N),  B(ldb, M);
Eigen::VectorXd alphaR(M),  alphaI(M),  beta(M);

//matrices with data
Eigen::MatrixXd A(lda,M), DE(ldde,M+1), C1(ldc1,M), VW(ldvw,M+1);

A << 3.1472,   1.3236,   4.5751,   4.5717,
4.0579,  -4.0246,   4.6489,  -0.1462,
-3.7301,  -2.2150,  -3.4239,   3.0028,
4.1338,   0.4688,   4.7059,  -3.5811;

DE << 0.0000,   0.0000,  -1.5510,  -4.5974,  -2.5127,
3.5071,   0.0000,   0.0000,   1.5961,   2.4490,
-3.1428,   2.5648,   0.0000,   0.0000,  -0.0596,
3.0340,   2.4892,  -1.1604,   0.0000,   0.0000;

C1 <<  0.6882,  -3.3782,  -3.3435,   1.8921,
-0.3061,   2.9428,   1.0198,   2.4815,
-4.8810,  -1.8878,  -2.3703,  -0.4946,
-1.6288,   0.2853,   1.5408,  -4.1618;

VW <<  -2.4013,  -2.7102,   0.3834,  -3.9335,   3.1730,
-3.1815,  -2.3620,   4.9613,   4.6190,   3.6869,
3.6929,   0.7970,  0.4986,  -4.9537,  -4.1556,
3.5303,   1.2206,  -1.4905,   0.1325,  -1.0022;

/* outputs of each parameter save for dwork,iwork to check correctness. */

F77NAME(dghutr)( &job, &compq1, &compq2, &N, A.data(), &lda, DE.data(), &ldde,  C1.data(), &ldc1, VW.data(), &ldvw,
Q1.data(), &ldq1,  Q2.data(), &ldq2,  B.data(), &ldb,
F.data(), &ldf,  C2.data(), &ldc2, alphaR.data(),  alphaI.data(),
beta.data(), iwork, &liwork, dwork, &ldwork, &info );
std::cout << "result: " << info << std::endl;
delete[] iwork;
delete[] dwork;
}

Компиляция выполняется с помощью (она использует много других вещей):

g++ -o eigensolver EigenSHEIGSolver.cpp -I/home/shared/eigen-eigen-1306d75b4a21  /home/shared/SHIRA/SHEVP/src/shheig64.a /home/shared/SHIRA/SLICOT_Lib/slicot64.a /home/shared/SHIRA/SLICOT_Lib/lpkaux64.a /home/shared/ATLAS/builddir/lib/libptlapack.a /home/shared/ATLAS/builddir/lib/libptcblas.a /home/shared/ATLAS/builddir/lib/libptf77blas.a /home/shared/ATLAS/builddir/lib/libatlas.a /home/shared/ATLAS/builddir/lib/libptcblas.a -lgfortran -lpthread

Увы, всякий раз, когда я запускаю полученный исполняемый файл, он дает мне:

 ** On entry to DGHUTR parameter number  8 had an illegal value

Мои знания ФОРТРАНА чрезвычайно ограничены, и приведенный выше код был написан в основном с использованием
YoLinux Учебник по смешиванию Фортрана и Си
а также
CRAY Docs
в качестве ссылок.
Как я понимаю, процедура сообщает об ошибке с ldde переменная. Я понятия не имею, почему, хотя.

Кто-нибудь может пролить свет на это для меня, пожалуйста?

Нотабене Согласно Eigen Docs: порядок хранения
Eigen хранит матрицы по умолчанию в порядке кол-мажора, поэтому он должен быть совместим с FORTRAN.
И подпрограмма Фортран DGHUTR является

SUBROUTINE DGHUTR( JOB, COMPQ1, COMPQ2, N, A, LDA, DE, LDDE, C1,
$                   LDC1, VW, LDVW, Q1, LDQ1, Q2, LDQ2, B, LDB, F,
$                   LDF, C2, LDC2, ALPHAR, ALPHAI, BETA, IWORK,
$                   LIWORK, DWORK, LDWORK, INFO )

Обновить: Вот вывод измененной подпрограммы DGHUTR (в основном добавлена ​​печать):

 JOB T
COMPQ1 I
COMPQ2 I
LDA          17179869188
LDDE          34359738372
LDC1          17179869188
LDVW         704374636548
LDQ1          34359738376
LDB          17179869188
LDF          17179869188
LDC2          17179869188
LIWORK                   20
LDWORK          85899346084
N          17179869192
LDDE          34359738372
INFO  6227620798727716864

Как видно, символы получены правильно, как есть LIWORKпри условии, что я скомпилирую с -O2 задавать. Я предполагаю, что есть что-то g++ делает что ломает параметры. Пытаясь отойти от gcc-5 в gcc-4.8 не решил проблему. Без оптимизации LDA Значение, похоже, меняется при каждом запуске программы, тогда как оно остается фиксированным при компиляции с -O2,