Skip to content

Commit

Permalink
Merge pull request #34 from nekStab/two_sided_eigs
Browse files Browse the repository at this point in the history 
Two sided eigs
  • Loading branch information
@loiseaujc
loiseaujc authored Feb 6, 2024
2 parents e4be6e6 + 87651b9 commit 199fadb
Show file tree
Hide file tree
Showing 3 changed files with 223 additions and 45 deletions.
42 changes: 28 additions & 14 deletions src/BaseKrylov.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,22 @@ module BaseKrylov
public :: arnoldi_factorization, &
lanczos_tridiagonalization, &
lanczos_bidiagonalization, &
nonsymmetric_lanczos_tridiagonalization
nonsymmetric_lanczos_tridiagonalization, &
initialize_krylov_subspace

contains

subroutine initialize_krylov_subspace(X)
!> Krylov subspace to be zeroed-out.
class(abstract_vector), intent(inout) :: X(:)
!> Internal variables.
integer :: i
do i = 1, size(X)
call X(i)%zero()
enddo
return
end subroutine initialize_krylov_subspace

!---------------------------------------------------------------------
!----- -----
!----- ARNOLDI FACTORIZATION FOR GENERAL SQUARE MATRICES -----
Expand Down Expand Up @@ -397,11 +409,13 @@ subroutine nonsymmetric_lanczos_tridiagonalization(A, V, W, T, info, kstart, ken
k_start = optval(kstart, 1)
k_end = optval(kend, kdim)
verbose = optval(verbosity, .false.)
tolerance = optval(tol, atol)
tolerance = optval(tol, rtol)

! --> Bi-orthogonalize the left and right starting vectors.
tmp = V(1)%dot(W(1)) ; beta = sqrt(abs(tmp)) ; gamma = sign(beta, tmp)
call V(1)%scal(1.0_wp / beta) ; call W(1)%scal(1.0_wp / gamma)
if (k_start .eq. 1) then
tmp = V(1)%dot(W(1)) ; beta = sqrt(abs(tmp)) ; gamma = sign(beta, tmp)
call V(1)%scal(1.0_wp / beta) ; call W(1)%scal(1.0_wp / gamma)
endif

! --> Nonsymmetric Lanczos iterations.
lanczos : do k = k_start, k_end
Expand All @@ -410,14 +424,20 @@ subroutine nonsymmetric_lanczos_tridiagonalization(A, V, W, T, info, kstart, ken
call A%rmatvec(W(k), W(k+1))

! --> Update diagonal entry of the nonsymmetric tridiagonal matrix.
alpha = W(k)%dot(V(k+1)) ; T(k, k) = alpha
alpha = V(k+1)%dot(W(k)) ; T(k, k) = alpha

! --> Lanczos three term recurrence.
call V(k+1)%axpby(1.0_wp, V(k), -alpha)
if (k > 1) call V(k+1)%axpby(1.0_wp, V(k-1), -gamma)
if (k > 1) call V(k+1)%axpby(1.0_wp, V(k-1), -T(k-1, k))

call W(k+1)%axpby(1.0_wp, W(k), -alpha)
if (k > 1) call W(k+1)%axpby(1.0_wp, W(k-1), -beta)
if (k > 1) call W(k+1)%axpby(1.0_wp, W(k-1), -T(k, k-1))

! --> Full re-biorthogonalization.
do i = 1, k
alpha = V(k+1)%dot(W(i)) ; call V(k+1)%axpby(1.0_wp, V(i), -alpha)
alpha = W(k+1)%dot(V(i)) ; call W(k+1)%axpby(1.0_wp, W(i), -alpha)
enddo

! --> Update the off-diagonal entries of the nonsymmetric tridiagonal matrix.
tmp = V(k+1)%dot(W(k+1)) ; beta = sqrt(abs(tmp)) ; gamma = sign(beta, tmp)
Expand All @@ -428,16 +448,10 @@ subroutine nonsymmetric_lanczos_tridiagonalization(A, V, W, T, info, kstart, ken
write(*, *) "INFO : Invariant subspaces have been computed (beta, gamma) = (", beta, ",", gamma, ")."
endif
else
! --> Full re-biorthogonalization.
do i = 1, k
alpha = V(k+1)%dot(W(i)) ; call V(k+1)%axpby(1.0_wp, V(i), -alpha)
alpha = W(k+1)%dot(V(i)) ; call W(k+1)%axpby(1.0_wp, W(i), -alpha)
enddo

! --> Normalization step.
call V(k+1)%scal(1.0_wp / beta) ; call W(k+1)%scal(1.0_wp / gamma)
endif

enddo lanczos

if (verbose) then
Expand Down
223 changes: 193 additions & 30 deletions src/IterativeSolvers.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ module IterativeSolvers

private
public :: eigs, eighs, gmres, save_eigenspectrum, svds, cg
public :: two_sided_eigs
public :: abstract_linear_solver

!------------------------------------------------
Expand Down Expand Up @@ -292,6 +293,130 @@ subroutine eighs(A, X, eigvecs, eigvals, residuals, info, nev, tolerance, verbos
return
end subroutine eighs

subroutine two_sided_eigs(A, V, W, rvecs, lvecs, eigvals, residuals, info, nev, tolerance, verbosity)
!> Linear Operator.
class(abstract_linop), intent(in) :: A
!> Left and right Krylov basis.
class(abstract_vector), dimension(:), intent(inout) :: V
class(abstract_vector), dimension(:), intent(inout) :: W
!> Coordinates of eigenvectors in Krylov bases, eigenvalues and associated residuals.
complex(kind=wp), dimension(size(V)-1, size(V)-1), intent(out) :: rvecs
complex(kind=wp), dimension(size(W)-1, size(W)-1), intent(out) :: lvecs
complex(kind=wp), dimension(size(V)-1) , intent(out) :: eigvals
real(kind=wp) , dimension(size(V)-1) , intent(out) :: residuals
!> Information flag.
integer, intent(out) :: info
!> Number of desired eigenvalues.
integer, optional, intent(in) :: nev
integer :: nev_, conv
!> Tolerance control.
real(kind=wp), optional, intent(in) :: tolerance
real(kind=wp) :: tol
!> Verbosity control.
logical, optional, intent(in) :: verbosity
logical :: verbose

!> Tridiagonal matrix.
real(kind=wp), dimension(size(V), size(V)) :: T
real(kind=wp) :: beta, gamma
!> Krylov subspace dimension.
integer :: kdim
!> Miscellaneous.
integer :: i, j, k
integer(int_size), dimension(size(V)-1) :: indices
real(kind=wp) , dimension(size(V)-1) :: abs_eigvals
real(kind=wp) :: alpha, tmp

!> Gram matrix.
complex(kind=wp), dimension(size(V)-1, size(V)-1) :: G
complex(kind=wp), dimension(size(V)-1, size(V)-1) :: H
real(kind=wp) :: direct_residual, adjoint_residual

!> Dimension of the Krylov subspaces.
kdim = size(V)-1

!> Deals with optional args.
verbose = optval(verbosity, .false.)
nev_ = optval(nev, kdim)
tol = optval(tolerance, rtol)

!> Initialize variables.
T = 0.0_wp ; residuals = 0.0_wp ; eigvals = cmplx(0.0_wp, 0.0_wp, kind=wp)
lvecs = cmplx(0.0_wp, 0.0_wp, kind=wp) ; rvecs = cmplx(0.0_wp, 0.0_wp, kind=wp)

!> Make sure the first Krylov vectors are bi-orthogonal.
tmp = V(1)%dot(W(1)) ; beta = sqrt(abs(tmp)) ; gamma = sign(beta, tmp)
call V(1)%scal(1.0_wp / beta) ; call W(1)%scal(1.0_wp / gamma)

do i = 2, size(V)
call V(i)%zero() ; call W(i)%zero()
enddo

lanczos : do k = 1, kdim
!> Lanczos step.
call nonsymmetric_lanczos_tridiagonalization(A, V, W, T, info, kstart=k, kend=k, verbosity=verbose)
if (info < 0) then
if (verbose) then
write(*, *) "INFO : Lanczos iteration failed. Exiting the eig subroutine."
write(*, *) " Lanczos exit code :", info
endif
info = -1
return
endif

!> Spectral decomposition of the k x k tridiagonal matrix.
eigvals = cmplx(0.0_wp, 0.0_wp, kind=wp)
rvecs = cmplx(0.0_wp, 0.0_wp, kind=wp) ; lvecs = cmplx(0.0_wp, 0.0_wp, kind=wp)
call evd(T(1:k, 1:k), rvecs(1:k, 1:k), eigvals(1:k), k)
lvecs(1:k, 1:k) = rvecs(1:k, 1:k) ; call cinv(lvecs(1:k, 1:k)) ; lvecs(1:k, 1:k) = transpose(conjg(lvecs(1:k, 1:k)))

!> Sort eigenvalues.
abs_eigvals = abs(eigvals) ; call sort_index(abs_eigvals, indices, reverse=.true.)
eigvals = eigvals(indices) ; rvecs = rvecs(:, indices) ; lvecs = lvecs(:, indices)

!> Compute residuals.
beta = abs(T(k+1, k)) ; gamma = abs(T(k, k+1)) ; alpha = V(k+1)%norm()
residuals(1:k) = compute_residual(beta*alpha, abs(rvecs(k, 1:k)))

G = 0.0_wp ; H = 0.0_wp
do i = 1, k
G(i, i) = V(i)%dot(V(i)) ; H(i, i) = W(i)%dot(W(i))
do j = i+1, k
G(i, j) = V(i)%dot(V(j))
G(j, i) = G(i, j)

H(i, j) = W(i)%dot(W(j))
H(j, i) = H(i, j)
enddo
enddo
residuals = 100.0_wp
do i = 1, k
alpha = V(k+1)%norm()
direct_residual = compute_residual(beta*alpha, abs(rvecs(k, i))) &
/ sqrt( real(dot_product(rvecs(:, k), matmul(G, rvecs(:, k)))) )
alpha = W(k+1)%norm()
adjoint_residual = compute_residual(gamma*alpha, abs(lvecs(k, i))) &
/ sqrt( real(dot_product(lvecs(:, k), matmul(H, lvecs(:, k)))) )
residuals(i) = max(direct_residual, adjoint_residual)
enddo

!> Check convergence.
conv = count(residuals(1:k) < tol)
write(*, *) "--- Iteration ", k, "---"
write(*, *) conv, "eigentriplets have converged."
write(*, *)
if (conv >= nev_) then
info = k
exit lanczos
endif

end do lanczos

call sort_index(residuals, indices, reverse=.false.)
eigvals = eigvals(indices) ; rvecs = rvecs(:, indices) ; lvecs = lvecs(:, indices)
return
end subroutine two_sided_eigs

subroutine evd(A, vecs, vals, n)
!> Lapack job.
character*1 :: jobvl = "N", jobvr = "V"
Expand Down Expand Up @@ -366,6 +491,38 @@ end subroutine dsyev
return
end subroutine hevd

subroutine rinv(A)
!> Matrix to invert (in-place)
real(kind=wp), intent(inout) :: A(:, :)
!> Lapack-related.
integer :: n, info
real(kind=wp) :: work(size(A, 1))
integer :: ipiv(size(A, 1))

!> Compute A = LU (in-place).
n = size(A, 1) ; call dgetrf(n, n, A, n, ipiv, info)
!> Compute inv(A).
call dgetri(n, A, n, ipiv, work, n, info)

return
end subroutine rinv

subroutine cinv(A)
!> Matrix to invert (in-place)
complex(kind=wp), intent(inout) :: A(:, :)
!> Lapack-related.
integer :: n, info
complex(kind=wp) :: work(size(A, 1))
integer :: ipiv(size(A, 1))

!> Compute A = LU (in-place).
n = size(A, 1) ; call zgetrf(n, n, A, n, ipiv, info)
!> Compute inv(A).
call zgetri(n, A, n, ipiv, work, n, info)

return
end subroutine cinv

!----------------------------------------------
!----- -----
!----- SINGULAR VALUE COMPUTATION -----
Expand Down Expand Up @@ -450,6 +607,8 @@ subroutine svds(A, U, V, uvecs, vvecs, sigma, residuals, info, nev, tolerance, v

enddo lanczos

info = k

return
end subroutine svds

Expand Down Expand Up @@ -613,6 +772,8 @@ subroutine gmres(A, b, x, info, preconditioner, options, transpose)
allocate(y(1:k_dim)) ; y = 0.0_wp
allocate(e(1:k_dim+1)) ; e = 0.0_wp

info = 0

! --> Initial Krylov vector.
if (trans) then
call A%rmatvec(x, V(1))
Expand All @@ -621,33 +782,25 @@ subroutine gmres(A, b, x, info, preconditioner, options, transpose)
endif
call V(1)%sub(b) ; call V(1)%scal(-1.0_wp)
beta = V(1)%norm() ; call V(1)%scal(1.0_wp / beta)

gmres_iterations : do i = 1, maxiter
! --> Zero-out variables.
H = 0.0_wp ; y = 0.0_wp ; e = 0.0_wp ; e(1) = beta
do j = 2, size(V)
call V(j)%zero()
enddo


! --> Arnoldi factorization.
arnoldi : do k = 1, k_dim
! --> Arnoldi factorization.
if (has_precond) then
! --> Apply preconditioner.
wrk = V(k) ; call precond%apply(wrk)
! --> Matrix-vector product.
if (trans) then
call A%rmatvec(wrk, V(k+1))
else
call A%matvec(wrk, V(k+1))
endif
wrk = V(k) ; if (has_precond) call precond%apply(wrk)

! --> Matrix-vector product.
if (trans) then
call A%rmatvec(wrk, V(k+1))
else
if (trans) then
call A%rmatvec(V(k), V(k+1))
else
call A%matvec(V(k), V(k+1))
endif
call A%matvec(wrk, V(k+1))
endif

! --> Gram-Schmid orthogonalization (twice is enough).
do j = 1, k
alpha = V(k+1)%dot(V(j)) ; call V(k+1)%axpby(1.0_wp, V(j), -alpha)
H(j, k) = alpha
Expand All @@ -656,31 +809,36 @@ subroutine gmres(A, b, x, info, preconditioner, options, transpose)
alpha = V(k+1)%dot(V(j)) ; call V(k+1)%axpby(1.0_wp, V(j), -alpha)
H(j, k) = H(j, k) + alpha
enddo

! --> Update Hessenberg matrix and normalize new Krylov vector.
H(k+1, k) = V(k+1)%norm()
if (H(k+1, k) > tol) then
call V(k+1)%scal(1.0_wp / H(k+1, k))
endif

! --> Least-squares problem.
call lstsq(H(1:k+1, 1:k), e(1:k+1), y(1:k))

! --> Compute residual.
beta = norm2(e(1:k+1) - matmul(H(1:k+1, 1:k), y(1:k)))
if (verbose) then
write(*, *) "INFO : GMRES residual after ", (i-1)*k_dim + k, "iteration : ", beta
! write(*, *) "INFO : GMRES residual after ", (i-1)*k_dim + k, "iteration : ", beta
write(*, *) "INFO : GMRES residual after ", info+1, "iteration : ", beta
endif

! --> Current number of iterations performed.
info = info + 1

! --> Check convergence.
if (beta**2 .lt. tol) then
if (beta .lt. tol) then
exit arnoldi
endif
enddo arnoldi

! --> Update solution.
k = min(k, k_dim)
if (allocated(dx) .eqv. .false.) allocate(dx, source=x)
call dx%zero() ; call get_vec(dx, V(1:k), y(1:k))
if (has_precond) then
call precond%apply(dx)
endif
if (allocated(dx) .eqv. .false.) allocate(dx, source=x) ; call dx%zero()
call get_vec(dx, V(1:k), y(1:k)) ; if (has_precond) call precond%apply(dx)
call x%add(dx)

! --> Recompute residual for sanity check.
Expand All @@ -695,14 +853,14 @@ subroutine gmres(A, b, x, info, preconditioner, options, transpose)
! --> Initialize new starting Krylov vector if needed.
beta = V(1)%norm() ; call V(1)%scal(1.0_wp / beta)

if (beta**2 .lt. tol) then
exit gmres_iterations
endif

! --> Exit GMRES if desired accuracy is reached.
if (beta .lt. tol) exit gmres_iterations

enddo gmres_iterations

! --> Deallocate variables.
deallocate(V, H, y, e)
! --> Convergence information.
write(*, *) "INFO : GMRES converged with residual ", beta
write(*, *) " Computation required ", info, "iterations"

return
end subroutine gmres
Expand Down Expand Up @@ -741,6 +899,11 @@ end subroutine dgels
!> Solve the least-squares problem.
call dgels(trans, m, n, nrhs, A_tilde, lda, b_tilde, ldb, work, lwork, info)

if (info > 0) then
write(*, *) "INFO: Error in LAPACK least-squares solver. Stopping the computation."
call exit()
endif

!> Return solution.
x = b_tilde(1:n)

Expand Down
Loading

0 comments on commit 199fadb

Please sign in to comment.