I'm currently trying to solve a particular eigenvalue problem (so-called gyroscopic eigenvalue problem) with large, sparse matrices (from FEM dsicretization). The programming language is C++.
The standard reference for EVP is ARPACK. Alas, it implements only "classical" Arnoldi processes, which is not suited for such problems (c.f. Structure Preserving Methods).
Recently I've found this Algorithm 961 reference, which also provides some code - in FORTRAN! So i've tried to include the DGHUTR routine in C++, to no avail. Below is the MWE, which is an adaptation of the test for DGHUTR (TDGHUTR.f) in 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
}
#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;
}
Compilation is done with (it uses a lot of other stuff):
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
Alas, whenever I run the resulting executable it gives me:
** On entry to DGHUTR parameter number 8 had an illegal value
My FORTRAN knowledge is extremely limited and the above code was written mainly using
YoLinux Tutorial mixing FORTRAN and C
and
CRAY Docs
as references.
As I understand it the routine reports an error with the ldde variable. I've got no clue why, though.
Can anyone shed some light on this for me, please?
N.B. As per Eigen Docs: storage order Eigen stores matrices per default in col-major order, thus it should be interfaceable with FORTRAN. And the FORTRAN subroutine DGHUTR is
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: Here's the output of the modified DGHUTR subroutine (basically added printing):
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
As one can see the characters are received correctly, as is LIWORK, provided I compile with -O2 set. I'm guessing there's something g++ does which breaks the parameters. Trying to revert from gcc-5 to gcc-4.8 did not solve the issue. With no optimization the LDA value seems to be changing on each run of the program, whereas it remains fixed when compiled with -O2.
I think I have found the source of the problem which has been plaguing me.
The dependence of the values received by the fortran routine on the optimization flags was kind of
a hint that there might be something wrong with how the stored variables are interpreted by C++ and
FORTRAN.
After looking for the specific value of 17179869188 and finding this
SO post
I tried playing with the compiler flags for the libraries.
When I fetched SLICOT I took the source and a library precompiled with gfortran for linux (slicot_linux_gfortran.tar.gz).
This latter one came with a make.inc with OPTS = -O2 -fpic -fdefault-integer-8
The SHHEVP routines contained the following comment in the make.inc
IMPORTANT: Use the options -fPIC -fdefault-integer-8 for 64bit
architectures.
So I did as advised - and that was the problem!
Removing -fdefault-integer-8 and recompiling both SLICOT and DGHUTR solved my problem. Now the code given above
compiles and the FORTRAN subroutine receives the correct values. The results of the computation are
consistent with the reference results provided with the DGHUTR source.
Incidentally, the majority of the SLICOT tests are now working. With the old flags the compilation of examples stopped at TAB01ND, which would always hang. Now I get down to TMB03LD, whose compilation fails with
IF( LSAME( COMPQ, 'C' ) .AND. NEIG.GT.0 ) THEN
1
Error: Operands of logical operator '.and.' at (1) are INTEGER(4)/LOGICAL(4)
But that is, for now, of no concern to me.