В настоящее время я пытаюсь решить конкретную проблему собственных значений (так называемую гироскопическую проблему собственных значений) с большими разреженными матрицами (из дискретизации FEM). Язык программирования - C ++.
Стандартным эталоном для EVP является ARPACK. Увы, он реализует только «классические» процессы Arnoldi, что не подходит для таких задач (см. Методы сохранения структуры).
Недавно я нашел ссылку на алгоритм 961, который также предоставляет некоторый код - на ФОРТРАНЕ! Итак, я попытался включить подпрограмму DGHUTR в C ++, но безуспешно. Ниже представлен MWE, который представляет собой адаптацию теста для DGHUTR (TDGHUTR.f) на 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;
}
Компиляция выполняется с помощью (он использует много других вещей):
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
. Но я понятия не имею, почему.
Кто-нибудь может пролить свет на это для меня, пожалуйста?
N.B. Согласно 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
. Я предполагаю, что есть что-то g++
, что нарушает параметры. Попытка вернуться с gcc-5
на gcc-4.8
не решила проблему. Без оптимизации значение LDA
, кажется, меняется при каждом запуске программы, тогда как оно остается неизменным при компиляции с -O2
.
LDDE
. Было бы целесообразно просто добавить операторprint
вDGHUTR
, чтобы увидеть, какое значение он на самом деле видит и чем отличается от значения, которое вы пытаетесь передать. Также проверьте, что значение, которое вы передаете из C ++, действительно правильно в документацииDGHUTR
. - person Vladimir F   schedule 23.02.2017-O2
, и внезапноLIWORK
имеет правильное значение20
, остальное по-прежнему мусор, и подпрограмма немедленно сообщает о параметре 27 (LIWORK), имеющем недопустимое значение. - person Nox   schedule 23.02.2017-O2
), я получаю ожидаемые значения. Если есть какой-либо мусор, ожидайте случайной ошибки, поскольку первое, что делаетDGHUTR
, - это пытается проверить параметры. Я также не понимаю, как / где у вас есть разреженные матрицы. - person Avi Ginsburg   schedule 23.02.2017-llapack -lcblas -lblas
), и она по-прежнему дает ту же ошибку. - person Nox   schedule 23.02.2017DGHUTR
ожидает в качестве параметра (-типы) в FORTRAN? - person Alex   schedule 24.02.2017