Fixed bug in eigs.

nekStab · Jul 2, 2024 · 3e70b56 · 3e70b56
1 parent b5e0dce
commit 3e70b56
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 5 deletions.
diff --git a/src/IterativeSolvers.f90 b/src/IterativeSolvers.f90
@@ -363,6 +363,8 @@ subroutine eigs_rsp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         integer :: indices(kdim_)
         real(sp) :: abs_eigvals(kdim_)
 
+        ! Re-compute eigenvalues and eigenvectors.
+        k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k))
         ! Sort eigenvalues.
         abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.)
         eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices)
@@ -373,7 +375,6 @@ subroutine eigs_rsp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         end block
 
         ! Construct eigenvectors.
-        k = min(k, kdim_)
         do i = 1, nev
             call X(i)%zero()
             do j = 1, k
@@ -493,6 +494,8 @@ subroutine eigs_rdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         integer :: indices(kdim_)
         real(dp) :: abs_eigvals(kdim_)
 
+        ! Re-compute eigenvalues and eigenvectors.
+        k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k))
         ! Sort eigenvalues.
         abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.)
         eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices)
@@ -503,7 +506,6 @@ subroutine eigs_rdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         end block
 
         ! Construct eigenvectors.
-        k = min(k, kdim_)
         do i = 1, nev
             call X(i)%zero()
             do j = 1, k
@@ -613,6 +615,8 @@ subroutine eigs_csp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         integer :: indices(kdim_)
         real(sp) :: abs_eigvals(kdim_)
 
+        ! Re-compute eigenvalues and eigenvectors.
+        k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k))
         ! Sort eigenvalues.
         abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.)
         eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices)
@@ -623,7 +627,6 @@ subroutine eigs_csp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         end block
 
         ! Construct eigenvectors.
-        k = min(k, kdim_)
         do i = 1, nev
             call X(i)%zero()
             do j = 1, k
@@ -733,6 +736,8 @@ subroutine eigs_cdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         integer :: indices(kdim_)
         real(dp) :: abs_eigvals(kdim_)
 
+        ! Re-compute eigenvalues and eigenvectors.
+        k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k))
         ! Sort eigenvalues.
         abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.)
         eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices)
@@ -743,7 +748,6 @@ subroutine eigs_cdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver
         end block
 
         ! Construct eigenvectors.
-        k = min(k, kdim_)
         do i = 1, nev
             call X(i)%zero()
             do j = 1, k

diff --git a/src/IterativeSolvers.fypp b/src/IterativeSolvers.fypp
@@ -264,6 +264,8 @@ contains
         integer :: indices(kdim_)
         real(${kind}$) :: abs_eigvals(kdim_)
 
+        ! Re-compute eigenvalues and eigenvectors.
+        k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k))
         ! Sort eigenvalues.
         abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.)
         eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices)
@@ -274,7 +276,6 @@ contains
         end block
 
         ! Construct eigenvectors.
-        k = min(k, kdim_)
         do i = 1, nev
             call X(i)%zero()
             do j = 1, k