Skip to content

Commit

Permalink
Fix generalized eig rwork size (#929)
Browse files Browse the repository at this point in the history
  • Loading branch information
@perazz
perazz authored Jan 30, 2025
2 parents 399f9a1 + 8a2c6e8 commit 036b091
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 22 deletions.
20 changes: 10 additions & 10 deletions src/stdlib_linalg_eigenvalues.fypp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,18 +139,18 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
module function stdlib_linalg_eigvals_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,err) result(lambda)
!! Return an array of eigenvalues of matrix A.
!> Input matrix A[m,n]
${rt}$, intent(in), dimension(:,:), target :: a
${rt}$, intent(in), target :: a(:,:)
#:if ei=='ggev'
!> Generalized problem matrix B[n,n]
${rt}$, intent(inout), dimension(:,:), target :: b
${rt}$, intent(inout), target :: b(:,:)
#:endif
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out) :: err
!> Array of eigenvalues
complex(${rk}$), allocatable :: lambda(:)

!> Create
${rt}$, pointer, dimension(:,:) :: amat#{if ei=='ggev'}#, bmat #{endif}#
${rt}$, pointer :: amat(:,:)#{if ei=='ggev'}#, bmat(:,:) #{endif}#
integer(ilp) :: m,n,k

!> Create an internal pointer so the intent of A won't affect the next call
Expand All @@ -172,16 +172,16 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
module function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#) result(lambda)
!! Return an array of eigenvalues of matrix A.
!> Input matrix A[m,n]
${rt}$, intent(in), dimension(:,:), target :: a
${rt}$, intent(in), target :: a(:,:)
#:if ei=='ggev'
!> Generalized problem matrix B[n,n]
${rt}$, intent(inout), dimension(:,:), target :: b
${rt}$, intent(inout), target :: b(:,:)
#:endif
!> Array of eigenvalues
complex(${rk}$), allocatable :: lambda(:)
!> Create
${rt}$, pointer, dimension(:,:) :: amat#{if ei=='ggev'}#, bmat #{endif}#
${rt}$, pointer :: amat(:,:)#{if ei=='ggev'}#, bmat(:,:) #{endif}#
integer(ilp) :: m,n,k
!> Create an internal pointer so the intent of A won't affect the next call
Expand All @@ -205,10 +205,10 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
!! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues,
!! and optionally right or left eigenvectors.
!> Input matrix A[m,n]
${rt}$, intent(inout), dimension(:,:), target :: a
${rt}$, intent(inout), target :: a(:,:)
#:if ei=='ggev'
!> Generalized problem matrix B[n,n]
${rt}$, intent(inout), dimension(:,:), target :: b
${rt}$, intent(inout), target :: b(:,:)
#:endif
!> Array of eigenvalues
complex(${rk}$), intent(out) :: lambda(:)
Expand All @@ -232,7 +232,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
character :: task_u,task_v
${rt}$, target :: work_dummy(1),u_dummy(1,1),v_dummy(1,1)
${rt}$, allocatable :: work(:)
${rt}$, dimension(:,:), pointer :: amat,umat,vmat#{if ei=='ggev'}#,bmat#{endif}#
${rt}$, pointer :: amat(:,:),umat(:,:),vmat(:,:)#{if ei=='ggev'}#,bmat(:,:)#{endif}#
#:if rt.startswith('complex')
real(${rk}$), allocatable :: rwork(:)
#:else
Expand Down Expand Up @@ -353,7 +353,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues

! Compute workspace size
#:if rt.startswith('complex')
allocate(rwork(2*n))
allocate(rwork( #{if ei=='ggev'}# 8*n #{else}# 2*n #{endif}# ))
#:else
allocate(lreal(n),limag(n))
#:endif
Expand Down
46 changes: 34 additions & 12 deletions test/linalg/test_linalg_eigenvalues.fypp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
module test_linalg_eigenvalues
use stdlib_linalg_constants
use stdlib_linalg_state
use stdlib_linalg, only: eig, eigh, eigvals, eigvalsh, diag
use stdlib_linalg, only: eig, eigh, eigvals, eigvalsh, diag, eye
use testdrive, only: error_type, check, new_unittest, unittest_type

implicit none (type,external)
Expand All @@ -21,27 +21,23 @@ module test_linalg_eigenvalues
allocate(tests(0))

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if rk!="xdp"
tests = [tests,new_unittest("test_eig_real_${ri}$",test_eig_real_${ri}$), &
new_unittest("test_eigvals_identity_${ri}$",test_eigvals_identity_${ri}$), &
new_unittest("test_eigvals_diagonal_B_${ri}$",test_eigvals_diagonal_B_${ri}$), &
new_unittest("test_eigvals_nondiagonal_B_${ri}$",test_eigvals_nondiagonal_B_${ri}$), &
new_unittest("test_eigh_real_${ri}$",test_eigh_real_${ri}$)]
#:endif
#: endfor

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if ck!="xdp"
tests = [tests,new_unittest("test_eig_complex_${ci}$",test_eig_complex_${ci}$), &
new_unittest("test_eig_generalized_complex_${ci}$",test_eigvals_generalized_complex_${ci}$)]
#:endif
new_unittest("test_eig_generalized_complex_${ci}$",test_eigvals_generalized_complex_${ci}$), &
new_unittest("test_eig_issue_927_${ci}$",test_issue_927_${ci}$)]
#: endfor

end subroutine test_eig_eigh

!> Simple real matrix eigenvalues
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if rk!="xdp"
subroutine test_eig_real_${ri}$(error)
type(error_type), allocatable, intent(out) :: error

Expand Down Expand Up @@ -239,12 +235,10 @@ module test_linalg_eigenvalues
if (allocated(error)) return
end subroutine test_eigvals_nondiagonal_B_${ri}$

#:endif
#:endfor

!> Simple complex matrix eigenvalues
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if ck!="xdp"
subroutine test_eig_complex_${ci}$(error)
type(error_type), allocatable, intent(out) :: error

Expand Down Expand Up @@ -309,8 +303,6 @@ module test_linalg_eigenvalues

lambda = eigvals(A, B, err=state)

print *, 'lambda = ',lambda

!> Expected eigenvalues
lres(1) = czero
lres(2) = 2*cone
Expand All @@ -324,10 +316,40 @@ module test_linalg_eigenvalues

end subroutine test_eigvals_generalized_complex_${ci}$

#:endif
! Generalized eigenvalues should not crash
subroutine test_issue_927_${ci}$(error)
type(error_type), allocatable, intent(out) :: error

${ct}$ :: A_Z(3,3),S_Z(3,3),vecs_r(3,3),eigs(3)
real(${ck}$) :: A_D(3,3),S_D(3,3)
type(linalg_state_type) :: state
integer :: i

! Set matrix
A_Z = reshape( [ [1, 6, 3], &
[9, 2, 1], &
[8, 3, 4] ], [3,3] )

S_Z = eye(3, mold=0.0_${ck}$)

A_D = real(A_Z)
S_D = real(S_Z)

call eig(A_D,S_D,eigs,right=vecs_r,err=state)
call check(error, state%ok(), 'test issue 927 (${ct}$): '//state%print())
if (allocated(error)) return

call eig(A_Z,S_Z,eigs,right=vecs_r,err=state) !Fails
call check(error, state%ok(), 'test issue 927 (${ct}$): '//state%print())
if (allocated(error)) return

end subroutine test_issue_927_${ci}$

#:endfor




end module test_linalg_eigenvalues

program test_eigenvalues
Expand Down

0 comments on commit 036b091

Please sign in to comment.