ARPACK называет процедуру SuperLU неправильной?

При использовании хорошо известной библиотеки ARPACK для решения больших, разреженных обобщенных задач на собственные значения (здесь я использую интерфейс C ++), я заметил очень странные решения при попытке вызвать

nev = 10 // number of eigenvalues to calculate
ARluSymGenEig<double> eigenvalueProblem('C', nev, A, B, 5., "SM");
eigenvalueProblem.FindEigenvectors();

здесь математически: Решите задачу Лапласа-Собственного значения в режиме Кэли на сетке 5×5 на единичном квадрате [0,1] ^ 2 с линейными элементами. «А» — соответствующая жесткость, а «В» — матрица масс. Граничные условия Дирхлета не вводятся, так что это проблема Неймана.
Меня интересуют только первые 10 самых маленьких собственных пар, поэтому «SM» (наименьшая величина); «5.» некоторый параметр сдвига.

Я получил этот вывод:

~/LapEVSol$ ./main SquareGrid.txt meshoutput.vtk
Reading eigenvalue problem from file SquareGrid.txt...
** On entry to dgstrs, parameter number  6 had an illegal value
** On entry to dgstrs, parameter number  6 had an illegal value
** On entry to dgstrs, parameter number  6 had an illegal value
...
** On entry to dgstrs, parameter number  6 had an illegal value
** On entry to dgstrs, parameter number  6 had an illegal value
** On entry to dgstrs, parameter number  6 had an illegal value

Done!
Dimension of the system            : 25
Number of 'requested' eigenvalues  : 10
Number of 'converged' eigenvalues  : 10
Number of Arnoldi vectors generated: 21
Number of iterations taken         : 10

Eigenvalues:
lambda[1]: -90.5854
lambda[2]: -16.2358
lambda[3]: -11.0635
lambda[4]: -8.95837
lambda[5]: -5.78369
lambda[6]: 16.5812
lambda[7]: 16.9097
lambda[8]: 36.4766
lambda[9]: 62.1445
lambda[10]: 66.4824

||A*x(0) - lambda(0)*B*x(0)||: 0.252506
||A*x(1) - lambda(1)*B*x(1)||: 0.75909
||A*x(2) - lambda(2)*B*x(2)||: 0.911785
||A*x(3) - lambda(3)*B*x(3)||: 0.901433
||A*x(4) - lambda(4)*B*x(4)||: 1.20555
||A*x(5) - lambda(5)*B*x(5)||: 0.980946
||A*x(6) - lambda(6)*B*x(6)||: 0.985718
||A*x(7) - lambda(7)*B*x(7)||: 0.26326
||A*x(8) - lambda(8)*B*x(8)||: 0.204725
||A*x(9) - lambda(9)*B*x(9)||: 0.113176

Конечно, этот вывод совершенно неверен. Все собственные значения должны быть неотрицательными и составлять около 0, 10, 10, 22, 47, 47 и т. Д. Сообщения об ошибках с помощью dgstrs () приводят к предположению, что ARPACK вызывает подпрограмму SuperLU dgstrs () неправильно для решить некоторую линейную систему. В коде ARPACK пытается вызвать

в arlspen.h в шаблоне функции void ARluSymPencil :: MultInvAsBv (ARTYPE v, ARTYPE * w): *

Create_Dense_Matrix(&RHS, A->nrows(), 1, w, A->nrows(), SLU_DN, SLU_GE);
gstrs(trans, &L, &U, permc, permr, &RHS, &stat, &info);

в суперлуке

inline void gstrs(trans_t trans, SuperMatrix *L, SuperMatrix *U,
int *perm_c, int *perm_r, SuperMatrix *B, SuperLUStat_t* stat, int *info)
{

if (L->Dtype == SLU_D) {       // calling the double precision routine.
dgstrs(trans,L,U,perm_c,perm_r,B,stat,info);
}
...
}inline void Create_Dense_Matrix(SuperMatrix* A, int m, int n, double* x,
int ldx, Stype_t S, Mtype_t M)
{

dCreate_Dense_Matrix(A,m,n,x,ldx,S,SLU_D,M);

} // Create_Dense_Matrix (double).

dgstrs () и dCreate_Dense_Matrix () являются подпрограммами SuperLU, написанными на Фортране, и поэтому их трудно отследить при отладке. Сообщение об ошибке появляется прямо в dgstrs () и говорит, что с шестым параметром что-то не так, то есть SuperMatrix B. B создан выше в функции Create_Dense_Matrix () — и действительно кажется немного странным, потому что с параметрами «.. ., A-> nrows (), 1, …, A-> nrows () «, он в основном создал плотную матрицу» m times 1 «с начальным измерением n (где m — количество строк, а n — количество столбцов входной матрицы A) — я бы ожидал, что это будет другим.

Даже в примерах, приведенных ARPACK ++, таких как lsymgreg.cc (который в основном вызывает те же функции, что и в моей задаче), вывод дает неверные отрицательные собственные значения с такими же ошибками dgstrs (). И все тестируемые задачи имеют действительно небольшие размеры (матрица 25×25), поэтому проблем со стабильностью быть не должно, а проблема симметричных собственных значений хорошо обусловлена.

Кто-нибудь когда-нибудь сталкивался с такой же проблемой с ARPACK / ARPACK ++ или знает, в чем причина?

редактировать: Вот еще немного информации о двух основных функциях для этой проблемы:
Первый FindArnoldiBasis () в argsym.h:

template<class ARFLOAT, class ARFOP, class ARFB>
int ARSymGenEig<ARFLOAT, ARFOP, ARFB>::FindArnoldiBasis()
{

ARFLOAT* temp;

if (this->mode != 5) {  // Using base function if not in Cayley mode.
return ARGenEig<ARFLOAT, ARFLOAT, ARFOP, ARFB>::FindArnoldiBasis();
}
else {

temp = new ARFLOAT[this->n+1];

if (!this->BasisOK) this->Restart();

// Changing to auto shift mode.

if (!this->AutoShift) {
ArpackError::Set(ArpackError::CHANGING_AUTOSHIFT, "FindArnoldiBasis");
this->AutoShift=true;
}

// ARPACK main loop.

while (!this->BasisOK) {

// Calling Aupp.

try { this->TakeStep(); }
catch (ArpackError) {
ArpackError(ArpackError::CANNOT_FIND_BASIS, "FindArnoldiBasis");
delete[] temp;
return 0;
}

switch (this->ido) {
case -1:

// Performing y <- B*x for the first time.

this->ipntr[3] = this->ipntr[2]+this->n; // not a clever idea, but...
(this->objB->*(this->MultBx))(&this->workd[this->ipntr[1]],&this->workd[this->ipntr[3]]);

case  1:

// Performing y <- OP*(A+sigma*B)*x, B*x is already available.

(this->objB->*MultAx)(&this->workd[this->ipntr[1]], temp);
axpy(this->n, this->sigmaR, &this->workd[this->ipntr[3]], 1, temp, 1);
(this->objOP->*(this->MultOPx))(temp, &this->workd[this->ipntr[2]]);
break;

case  2:

// Performing y <- B*x.

(this->objB->*(this->MultBx))(&this->workd[this->ipntr[1]],&this->workd[this->ipntr[2]]);

}
}

delete[] temp;

return this->nconv;
}

// FindArnoldiBasis.

Этот метод вызывает (в случае 1) MultOPx (), какие звонки MultInvAsBv () и поэтому gstrs ():

template<class ARTYPE>
void ARluSymPencil<ARTYPE>::MultInvAsBv(ARTYPE* v, ARTYPE* w)
{

// Quitting the function if AsB was not factored.

if (!IsFactored()) {
throw ArpackError(ArpackError::NOT_FACTORED_MATRIX,
"ARluSymPencil::MultInvAsBv");
}

// Solving AsB.w = v.

int         info;
SuperMatrix RHS;

copy(A->nrows(), v, 1, w, 1);
Create_Dense_Matrix(&RHS, A->nrows(), 1, w, A->nrows(), SLU_DN, SLU_GE);
//  gstrs("N", &L, &U, permr, permc, &RHS, &info);
trans_t trans = NOTRANS;
StatInit(&stat);

gstrs(trans, &L, &U, permc, permr, &RHS, &stat, &info);

Destroy_SuperMatrix_Store(&RHS); // delete RHS.Store;

} // MultInvAsBv.

Переменные L, U, permc и permr находятся в классе ARluSymPencil, в который включена вышеуказанная функция-член:

template<class ARTYPE>
class ARluSymPencil
{

protected:

bool                   factored;
int*                   permc;
int*                   permr;
char                   part;
char                   uplo;
ARluSymMatrix<ARTYPE>* A;
ARluSymMatrix<ARTYPE>* B;
SuperMatrix            L;
SuperMatrix            U;
SuperLUStat_t stat;

(SuperMatrix и некоторые другие типы определены в SuperLU-библиотеке)