Add `nev` argument to `eigs`, `eighs`, and `svds`

At the moment, we have a very basic implementation for eigs and others. We compute the whole Krylov basis (or only part of it if we are lucky with the invariant subspace) and then compute the eigenvalues/eigenvectors. No restarting procedure whatsoever nor any control over how many eigenvalues the user would like to converge.

We thus need two things:

an on-the-fly computation of the eigenvalue residuals so that we can stop the Arnoldi iteration once we have reached n_converged == nev.
a restart procedure such as Krylov-Schur to make sure we keep real arithmetic (in contrast to implicitly restarted arnoldi which can use complex shift for the restart).

I have seen that you have already sketched a krylov_schur subroutine. I am not entirely sure we need a dedicated subroutine in KrylovDecomp.f90 for that or simply a subroutine like krylov_schur_restart which would compute the Schur decomposition of the Hessenberg matrix along with selecting the desired eigenvalues and rotate/truncate the Krylov basis accordingly. If so, the eigs subroutine would read

subroutine eigs(A, X, eigvecs, eigvals, residuals, nev, info, opts)
    !> Standard arguments
    ....
    !> Number of desired eigenpairs.
    integer, optional, intent(in) :: nev
    integer :: nev_

    nev_ = optval(nev, 6) ! 6 is arbitrary but I think it is the default in scipy.sparse.linalg.eigs
    kstart = 1

1   do k = kstart, kdim
        !> Arnoldi step.
        call arnoldi_factorization(A, X, H, info, kstart=k, kend=k+1) ; beta = H(k+1, k)
        !> Compute eigenvalues of H.
        call lapack_routine(H(1:k, 1:k), eigvals(1:k), eigvecs(1:k, 1:k)) ! Eigenvalues need to be sorted as well.
        !> Compute residuals.
        residuals(1:k) = compute_residual(beta, abs(eigvecs(k, :))
        !> Check for the number of converged eigenvalues.
        has_converged = check_convergence(residuals, nev_)
        if (has_converged) exit
    enddo

    !> Check for convergence.
    if (has_converged) then
        continue
    else
        call krylov_schur_restart(X, H, kstart)
GOTO 1
    endif

    ....

    return
end subroutine eigs

I don't quite like the GOTO statement but you get the point. We would also need to keep track of the total number of matrix-vector products computed to stop the computation if we exceed a user-defined limit. This limit could be included in the opts (see #25).

Defining a krylov_schur_restart procedure rather than a dedicated krylov_schur_factorization like we did in nekStab would actually be more re-usable, e.g. we can directly re-use krylov_schur_restart for the Hermitian case (eighs) without having to write krylov_schur_spd_factorization for instance.

What do you think?

nekStab / LightKrylov

Add `nev` argument to `eigs`, `eighs`, and `svds` #27