我目前正在尝试使用大型稀疏矩阵(来自FEM离散化)来解决特定的特征值问题(所谓的陀螺仪特征值问题)。编程语言是C ++。
EVP的标准参考是ARPACK。las,它仅实现“经典” Arnoldi流程,不适用于此类问题(参见“结构保留方法”)。
最近,我发现了这个Algorithm 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
,每当我运行生成的可执行文件时,它就会给我:
** On entry to DGHUTR parameter number 8 had an illegal value
我的FORTRAN知识非常有限,上面的代码主要是使用YoLinux教程编写的,其中混合了FORTRAN和C和CRAY Docs作为参考。据我了解,例程会报告ldde
变量错误。我不知道为什么。
请问有人可以帮我一下吗?
注意:根据Eigen Docs:存储顺序Eigen默认以大写顺序存储矩阵,因此它应可与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
set编译,就可以正确接收到字符。我猜想有什么事情g++
会破坏参数。尝试从还原gcc-5
为gcc-4.8
不能解决问题。在没有优化的情况下,该LDA
值似乎在程序的每次运行时都在变化,而在使用编译时,该值保持不变-O2
。
我想我已经找到困扰我的问题根源。fortran例程接收到的值对优化标志的依赖性有点暗示,C ++和FORTRAN解释存储的变量的方式可能有问题。在查找的特定值17179869188
并找到此SO帖子之后,我尝试使用这些库的编译器标志。
当我获取SLICOT时,我获取了源代码和一个使用gfortran针对Linux(slicot_linux_gfortran.tar.gz
)预编译的库。后一个带有make.incOPTS = -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)
但是,就目前而言,我不关心。
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句