CODEWITHSUNDEEP
c++lapacklapacke
Boki
Boki

Reputation: 193

LAPACKE_zheevx() failed to converge -- how to increase ABSTOL with 2*DLAMCH('S') in C++?

This is a question about setting proper tolerance ("abstol") for convergence of eigenvalue calculation with LAPACKE_zheevx() function in C++.

When LAPACKE_zheev() fails to converge when calculating eigenvalues/eigenvectors with the default value of "abstol" (i.e. abstol=-1), the LAPACK manual says to set abstol=2*DLAMCH('S'). However, DLAMCH is Fortran function and I use C++ which does not recognize it as a valid C++ function. Could anyone please help me how to properly set "abstol=2*DLAMCH('S')" when using LAPACK with C++ (i.e. when using LAPACKE)?

Thanks very much in advance!!

Background: LAPACKE is C++ interface for LAPACK (Fortran library for numerical algebra). LAPACKE_zheevx() is LAPACKE's C++ interface for LAPACK's ZHEEVX() function.

Keywords: LAPACK, LAPACKE, C++, ABSTOL, DLAMCH, CONVERGENCE, EIGENVALUES, EIGENVECTORS

Upvotes: 1

Views: 283

Answers (2)

davidhigh
davidhigh

Reputation: 15508

In addition to the answer of Jon Purdy (which is practically the way to go):

Here is the definition of dlamch.f in Fortran. Going through that function, and getting the Lapack-intrinsic functions from here, one sees that for the input 's' the function translated to C++ becomes:

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;
}

Don't ask me, however, why the one == rnd step is needed.

Upvotes: 2

Jon Purdy
Jon Purdy

Reputation: 55059

I don’t know anything about these libraries, but a bit of Googling reveals there is a corresponding LAPACKE_dlamch() function:

double LAPACKE_dlamch(char cmach)

So you should be able to just pass LAPACKE_dlamch('S') as the abstol (12th) parameter of 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
)

Upvotes: 2

Related Questions