从C 调用Fortran子例程会产生非法参数值
Calling FORTRAN subroutine from c++ yields illegal parameter value
我目前正在尝试解决特定的特征值问题(所谓的陀螺仪特征值问题)具有较大的稀疏矩阵(来自fem dsicretrization)。编程语言是C 。
EVP的标准参考为ARPACK。las,它仅实现"古典" Arnoldi过程,这不适合此类问题(C.F.结构保存方法)。
最近,我发现了此算法961参考,它也提供了一些代码 - 在Fortran中!因此,我尝试将DGHUTR例程包括在C 中,但无济于事。以下是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
a,每当我运行由此产生的可执行文件时,它就会给我:
** On entry to DGHUTR parameter number 8 had an illegal value
我的Fortran知识非常有限,上述代码主要使用Yolinux教程混合Fortran和C和克雷文档作为参考。据我了解,例程报告了ldde
变量的错误。我不知道为什么。
任何人都可以为我提供一些启示吗?
n.b。根据eigen文档:存储顺序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
,只要我用 -O2
集编译,就可以看到字符和LIWORK
一样正确接收。我猜有g++
会破坏参数的事情。试图从gcc-5
恢复到gcc-4.8
并不能解决问题。没有优化,LDA
值似乎在程序的每次运行上都在变化,而使用-O2
编译时,它仍然固定。
我想我发现了困扰我的问题的根源。Fortran例程对优化标志收到的值的依赖性是一种一个暗示,如何通过C 解释存储的变量有问题Fortran。在寻找17179869188
的特定值并找到此之后所以发帖我尝试使用库的编译器标志。
当我获取切片时,我将源和图书馆与Gfortran进行了预编译,用于Linux(slicot_linux_gfortran.tar.gz
)。后一个带有 OPTS = -O2 -fpic -fdefault-integer-8
的制作SHHEVP例程包含以下评论。
IMPORTANT: Use the options -fPIC -fdefault-integer-8 for 64bit
architectures.
所以我按照建议做了 - 这就是问题!
删除-fdefault-integer-8
并重新编译切片和Dghutr解决了我的问题。现在上面给出的代码编译和Fortran子例程接收到正确的值。计算的结果是与DGHUTR源提供的参考结果一致。
顺便说一句,大多数切片测试现在正在起作用。用旧的旗帜汇编示例停在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)
,但现在我不关心。
- 如何反转整数参数包
- 使用C++库在Android项目中修改gradle中的cmake参数,用于插入指令的测试
- 如何使用默认参数等选择模板专业化
- 模板参数替换失败,并且未完成隐式转换
- 具有默认模板参数的多态类的模板推导失败
- lambda参数转换为constexpr技巧,然后获取带链接的数组
- 将数组作为参数传递给函数安全吗?作为第三方职能部门,可以探索他们想要的之外的其他元素
- 函数调用中参数的顺序重要吗
- 部分定义/别名模板模板参数
- 模板-模板参数推导:三个不同的编译器三种不同的行为
- 使用不带参数的函数访问结构元素
- 基于另一个成员参数将函数调用从类传递给它的一个成员
- 无法编译 rtmidi 测试 cmidiin.cpp 文件, 非法指令
- 如何在OMNET++中指定与命令行参数组合的输出文件名
- 如何使用Luacneneneba API正确读取字符串和表参数
- 无法让我的模板工作.非法使用显式模板参数
- 为什么非模板化函数具有相同的名称和参数但不同的返回类型是非法的?(但对模板函数合法吗?)
- C++ 模板专用化 - 非类型模板参数 '__formal 的非法类型
- VS2015 - "F"的匹配错误:非类型模板参数"F"的非法类型
- 计算结果为"void..."的非类型参数包不是非法的吗?