lapack — LAPACKE_zheevx () не удалось сойтись — как увеличить ABSTOL с 2 * DLAMCH (‘S’) в C ++?
Это вопрос установки правильного допуска («abstol») для сходимости вычисления собственных значений с помощью функции LAPACKE_zheevx () в C ++.
Когда LAPACKE_zheev () не может сойтись при вычислении собственных значений / собственных векторов со значением по умолчанию «abstol» (то есть abstol = -1), в руководстве LAPACK сказано установить abstol = 2 * DLAMCH (‘S’). Тем не менее, DLAMCH — это функция Фортрана, и я использую C ++, который не распознает ее как допустимую функцию C ++. Может ли кто-нибудь помочь мне, как правильно установить «abstol = 2 * DLAMCH (‘S’)» при использовании LAPACK с C ++ (т.е. при использовании LAPACKE)?
Большое спасибо заранее!!
Фон:
LAPACKE — это интерфейс C ++ для LAPACK (библиотека Фортрана для числовой алгебры). LAPACKE_zheevx () — это интерфейс C ++ LAPACKE для функции LAPACK ZHEEVX ().
Ключевые слова:
LAPACK, LAPACKE, C ++, ABSTOL, DLAMCH, конвергенция, собственные значения, собственные значения
Решение
Я ничего не знаю об этих библиотеках, но немного Google показывает, что есть соответствующие LAPACKE_dlamch() функция:
double LAPACKE_dlamch(char cmach)
Таким образом, вы должны быть в состоянии просто пройти LAPACKE_dlamch('S') как abstol (12-й) параметр LAPACKE_zheevx():
lapack_int LAPACKE_zheevx (
int matrix_layout,
char jobz,
char range,
char uplo,
lapack_int n,
lapack_complex_double *a,
lapack_int lda,
double vl,
double vu,
lapack_int il,
lapack_int iu,
double abstol,
// ^^^^^^^^^^^^^
lapack_int *m,
double *w,
lapack_complex_double *z,
lapack_int ldz,
lapack_int *ifail
)
Другие решения
В дополнение к ответу Джона Пурди (который практически подходит):
Вот это определение dlamch.f в Фортране. Проходя через эту функцию, и получая внутренние функции Лапака от Вот, видно, что для ввода 's' функция, переведенная на C ++, становится:
auto dlamch_s()
{
auto one = 1.0;
auto rnd = one;
auto eps = (one == rnd ? 0.5 : 1.0) * std::numeric_limits<double>::epsilon();
auto sfmin = std::numeric_limits<double>::min();
auto small = one / std::numeric_limits<double>::max();
if(small>=sfmin)
{
sfmin = small*(one + eps);
}
return sfmin;
}
Не спрашивайте меня, почему one == rnd шаг нужен.