从 c++ 调用 FORTRAN 子例程会产生非法参数值

Posted 2023-02-16

【中文标题】从 c++ 调用 FORTRAN 子例程会产生非法参数值

【英文标题】：Calling FORTRAN subroutine from c++ yields illegal parameter value

【发布时间】：2017-02-23 10:03:56

【问题描述】：

我目前正在尝试解决一个特定的特征值问题（所谓的陀螺特征值问题） 具有大型稀疏矩阵（来自 FEM 离散化）。编程语言是 C++。

EVP 的标准参考是 ARPACK。唉，它只实现了“经典”的阿诺尔迪流程， 这不适合此类问题（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 变量的错误。不过，我不知道为什么。

谁能帮我解释一下？

注意根据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 )

更新：这是修改后的 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 编译时它保持不变。

【问题讨论】：

参数号8为LDDE。值得将 print 语句添加到 DGHUTR 中，以查看它实际看到的值以及它与您尝试传递的值有何不同。还要检查您从 C++ 传递的值在 DGHUTR 的文档中是否正确。

完成，输出为（n，lda，ldde ...）：1371798691817129986918171799861881717998691817179986988

好的。有趣的是：我重新编译了设置-O2的代码，突然LIWORK的正确值20，其余的仍然是垃圾，例程立即报告参数27（LIWORK）具有非法值

很奇怪。当我构建它时（即使没有-O2），我得到了预期的值。如果有任何垃圾，那么会出现一些随机错误，因为DGHUTR 所做的第一件事就是尝试验证参数。我也看不到你如何/在哪里有任何稀疏矩阵。

还没有稀疏矩阵。我只是想为此使用这个例程，但据我所知，这个例程需要密集矩阵，所以这将是一个不可行的 - 可能。我已经使用系统 BLAS 和 Lapack (-llapack -lcblas -lblas) 重新编译了库，但仍然产生相同的错误。

【参考方案1】：

我想我已经找到了困扰我的问题的根源。 fortran 例程接收到的值对优化标志的依赖性有点像 提示 C++ 解释存储变量的方式可能存在问题，并且 福兰。 在寻找17179869188的具体值并找到这个之后 SO post 我尝试使用库的编译器标志。

当我获取 SLICOT 时，我获取了源代码和一个用 gfortran for linux (slicot_linux_gfortran.tar.gz) 预编译的库。 后一个带有带有OPTS = -O2 -fpic -fdefault-integer-8的make.inc 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)

但是，目前，这与我无关。