Skip to content

Fix generalized eig rwork size #929

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Jan 30, 2025
Merged

Fix generalized eig rwork size #929

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e84f91a
fix `ggev` rwork size
perazz Jan 30, 2025
08becd8
add `ggev` tests
perazz Jan 30, 2025
a91ae04
lift `xdp` restriction
perazz Jan 30, 2025
8a2c6e8
remove `dimension` attributes
perazz Jan 30, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 10 additions & 10 deletions src/stdlib_linalg_eigenvalues.fypp
View file
Open in desktop
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
Open in desktop
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
Loading