从 C++ 调用 FORTRAN 子例程产生非法参数值
Calling FORTRAN subroutine from c++ yields illegal parameter value
我目前正在尝试解决一个特定的特征值问题（所谓的陀螺特征值问题）
具有大而稀疏的矩阵（来自 FEM 离散化）。编程语言为C++。
EVP 的标准参考是 ARPACK。唉，它只实现了 "classical" Arnoldi 进程，
不适用于此类问题 (c.f.Structure Preserving Methods).
最近我发现了这个 Algorithm 961 参考，它还提供了一些代码 - 在 FORTRAN 中！
所以我试图在 C++ 中包含 DGHUTR 例程，但无济于事。
下面是 MWE，它是对 C++ 中 DGHUTR (TDGHUTR.f) 测试的改编：
#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
}
#endif
int 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
我的FORTRAN知识极其有限，上面的代码主要是用
YoLinux Tutorial mixing FORTRAN and C
和
CRAY Docs
作为参考。
据我了解，例程报告 ldde 变量错误。不过，我不知道为什么。
任何人都可以帮我解释一下吗？
N.B。根据 Eigen Docs: storage order
Eigen 以 col-major 顺序默认存储矩阵，因此它应该可以与 FORTRAN 交互。
而 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 )
Update：这里是修改后的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.
编译时它保持固定
我想我找到了一直困扰我的问题的根源。
Fortran 例程收到的值对优化标志的依赖性有点像
暗示 C++ 和
语言
查找17179869188的具体值后发现这个
我尝试为库 .
使用编译器标志
当我获取 SLICOT 时，我获取了源代码和一个使用 gfortran 为 linux (slicot_linux_gfortran.tar.gz) 预编译的库。
后一个带有 make.inc 和 OPTS = -O2 -fpic -fdefault-integer-8
SHHEVP 例程在 make.inc
中包含以下注释
IMPORTANT: Use the options -fPIC -fdefault-integer-8 for 64bit
architectures.
所以我按照建议做了 - 这就是问题所在！
删除 -fdefault-integer-8 并重新编译 SLICOT 和 DGHUTR 解决了我的问题。现在上面给出的代码
编译并且 FORTRAN 子例程接收到正确的值。计算结果为
与 DGHUTR 源提供的参考结果一致。
顺便说一句，大多数 SLICOT 测试现在都可以使用。使用旧标志编译示例
停在 TAB01ND，它总是会挂起。现在我进入 TMB03LD，它的编译失败
IF( LSAME( COMPQ, 'C' ) .AND. NEIG.GT.0 ) THEN
1
Error: Operands of logical operator '.and.' at (1) are INTEGER(4)/LOGICAL(4)
但目前我不关心这个。
我目前正在尝试解决一个特定的特征值问题（所谓的陀螺特征值问题） 具有大而稀疏的矩阵（来自 FEM 离散化）。编程语言为C++。
EVP 的标准参考是 ARPACK。唉，它只实现了 "classical" Arnoldi 进程， 不适用于此类问题 (c.f.Structure Preserving Methods).
最近我发现了这个 Algorithm 961 参考，它还提供了一些代码 - 在 FORTRAN 中！ 所以我试图在 C++ 中包含 DGHUTR 例程，但无济于事。 下面是 MWE，它是对 C++ 中 DGHUTR (TDGHUTR.f) 测试的改编：
#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
}
#endif
int 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
我的FORTRAN知识极其有限，上面的代码主要是用
YoLinux Tutorial mixing FORTRAN and C
和
CRAY Docs
作为参考。
据我了解，例程报告 ldde 变量错误。不过，我不知道为什么。
任何人都可以帮我解释一下吗？
N.B。根据 Eigen Docs: storage order Eigen 以 col-major 顺序默认存储矩阵，因此它应该可以与 FORTRAN 交互。 而 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 )
Update：这里是修改后的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.
我想我找到了一直困扰我的问题的根源。
Fortran 例程收到的值对优化标志的依赖性有点像
暗示 C++ 和
语言
查找17179869188的具体值后发现这个
当我获取 SLICOT 时，我获取了源代码和一个使用 gfortran 为 linux (slicot_linux_gfortran.tar.gz) 预编译的库。
后一个带有 make.inc 和 OPTS = -O2 -fpic -fdefault-integer-8
SHHEVP 例程在 make.inc
IMPORTANT: Use the options -fPIC -fdefault-integer-8 for 64bit
architectures.
所以我按照建议做了 - 这就是问题所在！
删除 -fdefault-integer-8 并重新编译 SLICOT 和 DGHUTR 解决了我的问题。现在上面给出的代码
编译并且 FORTRAN 子例程接收到正确的值。计算结果为
与 DGHUTR 源提供的参考结果一致。
顺便说一句，大多数 SLICOT 测试现在都可以使用。使用旧标志编译示例 停在 TAB01ND，它总是会挂起。现在我进入 TMB03LD，它的编译失败
IF( LSAME( COMPQ, 'C' ) .AND. NEIG.GT.0 ) THEN
1
Error: Operands of logical operator '.and.' at (1) are INTEGER(4)/LOGICAL(4)
但目前我不关心这个。