Skip to content

Commit cc129c6

Browse files
perazzperazz
authoredJun 30, 2024··
linalg: Eigenvalues and Eigenvectors (#816)
·
v0.7.0
2 parents d996e43 + d433869 commit cc129c6

11 files changed

+1287
-2
lines changed
 

‎doc/specs/stdlib_linalg.md

+176-1
Original file line numberDiff line numberDiff line change
@@ -687,6 +687,8 @@ Expert (`Pure`) interface:
687687

688688
`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
689689

690+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
691+
690692
### Return value
691693

692694
For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.
@@ -899,6 +901,180 @@ Exceptions trigger an `error stop`.
899901
{!example/linalg/example_determinant2.f90!}
900902
```
901903

904+
## `eig` - Eigenvalues and Eigenvectors of a Square Matrix
905+
906+
### Status
907+
908+
Experimental
909+
910+
### Description
911+
912+
This subroutine computes the solution to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank, `real` or `complex` matrix.
913+
914+
Result array `lambda` returns the eigenvalues of \( A \). The user can request eigenvectors to be returned: if provided, on output `left` will contain the left eigenvectors, `right` the right eigenvectors of \( A \).
915+
Both `left` and `right` are rank-2 arrays, where eigenvectors are stored as columns.
916+
The solver is based on LAPACK's `*GEEV` backends.
917+
918+
### Syntax
919+
920+
`call ` [[stdlib_linalg(module):eig(interface)]] `(a, lambda [, right] [,left] [,overwrite_a] [,err])`
921+
922+
### Arguments
923+
924+
`a` : `real` or `complex` square array containing the coefficient matrix. If `overwrite_a=.false.`, it is an `intent(in)` argument. Otherwise, it is an `intent(inout)` argument and is destroyed by the call.
925+
926+
`lambda`: Shall be a `complex` or `real` rank-1 array of the same kind as `a`, containing the eigenvalues, or their `real` component only. It is an `intent(out)` argument.
927+
928+
`right` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the right eigenvectors of `a`. It is an `intent(out)` argument.
929+
930+
`left` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the left eigenvectors of `a`. It is an `intent(out)` argument.
931+
932+
`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
933+
934+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
935+
936+
### Return value
937+
938+
Raises `LINALG_ERROR` if the calculation did not converge.
939+
Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
940+
Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts.
941+
If `err` is not present, exceptions trigger an `error stop`.
942+
943+
### Example
944+
945+
```fortran
946+
{!example/linalg/example_eig.f90!}
947+
```
948+
949+
## `eigh` - Eigenvalues and Eigenvectors of a Real symmetric or Complex Hermitian Square Matrix
950+
951+
### Status
952+
953+
Experimental
954+
955+
### Description
956+
957+
This subroutine computes the solution to the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \),
958+
where \( A \) is a square, full-rank, `real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix.
959+
960+
Result array `lambda` returns the `real` eigenvalues of \( A \). The user can request the orthogonal eigenvectors
961+
to be returned: on output `vectors` may contain the matrix of eigenvectors, returned as a column.
962+
963+
Normally, only the lower triangular part of \( A \) is accessed. On input, `logical` flag `upper_a`
964+
allows the user to request what triangular part of the matrix should be used.
965+
966+
The solver is based on LAPACK's `*SYEV` and `*HEEV` backends.
967+
968+
### Syntax
969+
970+
`call ` [[stdlib_linalg(module):eigh(interface)]] `(a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])`
971+
972+
### Arguments
973+
974+
`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call.
975+
976+
`lambda`: Shall be a `complex` rank-1 array of the same precision as `a`, containing the eigenvalues. It is an `intent(out)` argument.
977+
978+
`vectors` (optional): Shall be a rank-2 array of the same type, size and kind as `a`, containing the eigenvectors of `a`. It is an `intent(out)` argument.
979+
980+
`upper_a` (optional): Shall be an input `logical` flag. If `.true.`, the upper triangular part of `a` will be accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument.
981+
982+
`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
983+
984+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
985+
986+
### Return value
987+
988+
Raises `LINALG_ERROR` if the calculation did not converge.
989+
Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
990+
If `err` is not present, exceptions trigger an `error stop`.
991+
992+
### Example
993+
994+
```fortran
995+
{!example/linalg/example_eigh.f90!}
996+
```
997+
998+
## `eigvals` - Eigenvalues of a Square Matrix
999+
1000+
### Status
1001+
1002+
Experimental
1003+
1004+
### Description
1005+
1006+
This function returns the eigenvalues to matrix \( A \): a square, full-rank, `real` or `complex` matrix.
1007+
The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \).
1008+
1009+
Result array `lambda` is `complex`, and returns the eigenvalues of \( A \).
1010+
The solver is based on LAPACK's `*GEEV` backends.
1011+
1012+
### Syntax
1013+
1014+
`lambda = ` [[stdlib_linalg(module):eigvals(interface)]] `(a, [,err])`
1015+
1016+
### Arguments
1017+
1018+
`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument.
1019+
1020+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
1021+
1022+
### Return value
1023+
1024+
Returns a `complex` array containing the eigenvalues of `a`.
1025+
1026+
Raises `LINALG_ERROR` if the calculation did not converge.
1027+
Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
1028+
If `err` is not present, exceptions trigger an `error stop`.
1029+
1030+
1031+
### Example
1032+
1033+
```fortran
1034+
{!example/linalg/example_eigvals.f90!}
1035+
```
1036+
1037+
## `eigvalsh` - Eigenvalues of a Real Symmetric or Complex Hermitian Square Matrix
1038+
1039+
### Status
1040+
1041+
Experimental
1042+
1043+
### Description
1044+
1045+
This function returns the eigenvalues to matrix \( A \): a where \( A \) is a square, full-rank,
1046+
`real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix.
1047+
The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \).
1048+
1049+
Result array `lambda` is `real`, and returns the eigenvalues of \( A \).
1050+
The solver is based on LAPACK's `*SYEV` and `*HEEV` backends.
1051+
1052+
### Syntax
1053+
1054+
`lambda = ` [[stdlib_linalg(module):eigvalsh(interface)]] `(a, [, upper_a] [,err])`
1055+
1056+
### Arguments
1057+
1058+
`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument.
1059+
1060+
`upper_a` (optional): Shall be an input logical flag. If `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed (default). It is an `intent(in)` argument.
1061+
1062+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
1063+
1064+
### Return value
1065+
1066+
Returns a `real` array containing the eigenvalues of `a`.
1067+
1068+
Raises `LINALG_ERROR` if the calculation did not converge.
1069+
Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
1070+
If `err` is not present, exceptions trigger an `error stop`.
1071+
1072+
### Example
1073+
1074+
```fortran
1075+
{!example/linalg/example_eigvalsh.f90!}
1076+
```
1077+
9021078
## `svd` - Compute the singular value decomposition of a rank-2 array (matrix).
9031079

9041080
### Status
@@ -989,4 +1165,3 @@ Exceptions trigger an `error stop`, unless argument `err` is present.
9891165
```fortran
9901166
{!example/linalg/example_svdvals.f90!}
9911167
```
992-

‎example/linalg/CMakeLists.txt

+4
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,10 @@ ADD_EXAMPLE(is_square)
1313
ADD_EXAMPLE(is_symmetric)
1414
ADD_EXAMPLE(is_triangular)
1515
ADD_EXAMPLE(outer_product)
16+
ADD_EXAMPLE(eig)
17+
ADD_EXAMPLE(eigh)
18+
ADD_EXAMPLE(eigvals)
19+
ADD_EXAMPLE(eigvalsh)
1620
ADD_EXAMPLE(trace)
1721
ADD_EXAMPLE(state1)
1822
ADD_EXAMPLE(state2)

‎example/linalg/example_eig.f90

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
! Eigendecomposition of a real square matrix
2+
program example_eig
3+
use stdlib_linalg, only: eig
4+
implicit none
5+
6+
integer :: i
7+
real, allocatable :: A(:,:)
8+
complex, allocatable :: lambda(:),vectors(:,:)
9+
10+
! Decomposition of this square matrix
11+
! NB Fortran is column-major -> transpose input
12+
A = transpose(reshape( [ [2, 2, 4], &
13+
[1, 3, 5], &
14+
[2, 3, 4] ], [3,3] ))
15+
16+
! Get eigenvalues and right eigenvectors
17+
allocate(lambda(3),vectors(3,3))
18+
19+
call eig(A, lambda, right=vectors)
20+
21+
do i=1,3
22+
print *, 'eigenvalue ',i,': ',lambda(i)
23+
print *, 'eigenvector ',i,': ',vectors(:,i)
24+
end do
25+
26+
end program example_eig

‎example/linalg/example_eigh.f90

+38
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
! Eigendecomposition of a real symmetric matrix
2+
program example_eigh
3+
use stdlib_linalg, only: eigh
4+
implicit none
5+
6+
integer :: i
7+
real, allocatable :: A(:,:),lambda(:),vectors(:,:)
8+
complex, allocatable :: cA(:,:),cvectors(:,:)
9+
10+
! Decomposition of this symmetric matrix
11+
! NB Fortran is column-major -> transpose input
12+
A = transpose(reshape( [ [2, 1, 4], &
13+
[1, 3, 5], &
14+
[4, 5, 4] ], [3,3] ))
15+
16+
! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors
17+
allocate(lambda(3),vectors(3,3))
18+
call eigh(A, lambda, vectors=vectors)
19+
20+
print *, 'Real matrix'
21+
do i=1,3
22+
print *, 'eigenvalue ',i,': ',lambda(i)
23+
print *, 'eigenvector ',i,': ',vectors(:,i)
24+
end do
25+
26+
! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors
27+
cA = A
28+
29+
allocate(cvectors(3,3))
30+
call eigh(cA, lambda, vectors=cvectors)
31+
32+
print *, 'Complex matrix'
33+
do i=1,3
34+
print *, 'eigenvalue ',i,': ',lambda(i)
35+
print *, 'eigenvector ',i,': ',cvectors(:,i)
36+
end do
37+
38+
end program example_eigh

‎example/linalg/example_eigvals.f90

+24
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
! Eigenvalues of a general real / complex matrix
2+
program example_eigvals
3+
use stdlib_linalg, only: eigvals
4+
implicit none
5+
6+
integer :: i
7+
real, allocatable :: A(:,:),lambda(:)
8+
complex, allocatable :: cA(:,:),clambda(:)
9+
10+
! NB Fortran is column-major -> transpose input
11+
A = transpose(reshape( [ [2, 8, 4], &
12+
[1, 3, 5], &
13+
[9, 5,-2] ], [3,3] ))
14+
15+
! Note: real symmetric matrix
16+
lambda = eigvals(A)
17+
print *, 'Real matrix eigenvalues: ',lambda
18+
19+
! Complex general matrix
20+
cA = cmplx(A, -2*A)
21+
clambda = eigvals(cA)
22+
print *, 'Complex matrix eigenvalues: ',clambda
23+
24+
end program example_eigvals

‎example/linalg/example_eigvalsh.f90

+25
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
! Eigenvalues of a real symmetric / complex hermitian matrix
2+
program example_eigvalsh
3+
use stdlib_linalg, only: eigvalsh
4+
implicit none
5+
6+
integer :: i
7+
real, allocatable :: A(:,:),lambda(:)
8+
complex, allocatable :: cA(:,:)
9+
10+
! Decomposition of this symmetric matrix
11+
! NB Fortran is column-major -> transpose input
12+
A = transpose(reshape( [ [2, 1, 4], &
13+
[1, 3, 5], &
14+
[4, 5, 4] ], [3,3] ))
15+
16+
! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors
17+
lambda = eigvalsh(A)
18+
print *, 'Symmetric matrix eigenvalues: ',lambda
19+
20+
! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors
21+
cA = A
22+
lambda = eigvalsh(cA)
23+
print *, 'Hermitian matrix eigenvalues: ',lambda
24+
25+
end program example_eigvalsh

‎src/CMakeLists.txt

+2-1
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,8 @@ set(fppFiles
2727
stdlib_linalg_outer_product.fypp
2828
stdlib_linalg_kronecker.fypp
2929
stdlib_linalg_cross_product.fypp
30-
stdlib_linalg_solve.fypp
30+
stdlib_linalg_eigenvalues.fypp
31+
stdlib_linalg_solve.fypp
3132
stdlib_linalg_determinant.fypp
3233
stdlib_linalg_state.fypp
3334
stdlib_linalg_svd.fypp

‎src/stdlib_linalg.fypp

+208
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,10 @@ module stdlib_linalg
1919
public :: det
2020
public :: operator(.det.)
2121
public :: diag
22+
public :: eig
23+
public :: eigh
24+
public :: eigvals
25+
public :: eigvalsh
2226
public :: eye
2327
public :: lstsq
2428
public :: lstsq_space
@@ -554,6 +558,210 @@ module stdlib_linalg
554558
#:endfor
555559
end interface
556560

561+
! Eigendecomposition of a square matrix: eigenvalues, and optionally eigenvectors
562+
interface eig
563+
!! version: experimental
564+
!!
565+
!! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for square matrix \( A \).
566+
!! ([Specification](../page/specs/stdlib_linalg.html#eig-eigenvalues-and-eigenvectors-of-a-square-matrix))
567+
!!
568+
!!### Summary
569+
!! Subroutine interface for computing eigenvalues and eigenvectors of a square matrix.
570+
!!
571+
!!### Description
572+
!!
573+
!! This interface provides methods for computing the eigenvalues, and optionally eigenvectors,
574+
!! of a general square matrix. Supported data types include `real` and `complex`, and no assumption is
575+
!! made on the matrix structure. The user may request either left, right, or both
576+
!! eigenvectors to be returned. They are returned as columns of a square matrix with the same size as `A`.
577+
!! Preallocated space for both eigenvalues `lambda` and the eigenvector matrices must be user-provided.
578+
!!
579+
!!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`.
580+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
581+
!!
582+
#:for rk,rt,ri in RC_KINDS_TYPES
583+
#:if rk!="xdp"
584+
module subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err)
585+
!! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues,
586+
!! and optionally right or left eigenvectors.
587+
!> Input matrix A[m,n]
588+
${rt}$, intent(inout), target :: a(:,:)
589+
!> Array of eigenvalues
590+
complex(${rk}$), intent(out) :: lambda(:)
591+
!> The columns of RIGHT contain the right eigenvectors of A
592+
complex(${rk}$), optional, intent(out), target :: right(:,:)
593+
!> The columns of LEFT contain the left eigenvectors of A
594+
complex(${rk}$), optional, intent(out), target :: left(:,:)
595+
!> [optional] Can A data be overwritten and destroyed?
596+
logical(lk), optional, intent(in) :: overwrite_a
597+
!> [optional] state return flag. On error if not requested, the code will stop
598+
type(linalg_state_type), optional, intent(out) :: err
599+
end subroutine stdlib_linalg_eig_${ri}$
600+
#:endif
601+
#:endfor
602+
#:for rk,rt,ri in REAL_KINDS_TYPES
603+
#:if rk!="xdp"
604+
module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err)
605+
!! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues,
606+
!! and optionally right or left eigenvectors. Returns an error if the eigenvalues had
607+
!! non-trivial imaginary parts.
608+
!> Input matrix A[m,n]
609+
${rt}$, intent(inout), target :: a(:,:)
610+
!> Array of real eigenvalues
611+
real(${rk}$), intent(out) :: lambda(:)
612+
!> The columns of RIGHT contain the right eigenvectors of A
613+
complex(${rk}$), optional, intent(out), target :: right(:,:)
614+
!> The columns of LEFT contain the left eigenvectors of A
615+
complex(${rk}$), optional, intent(out), target :: left(:,:)
616+
!> [optional] Can A data be overwritten and destroyed?
617+
logical(lk), optional, intent(in) :: overwrite_a
618+
!> [optional] state return flag. On error if not requested, the code will stop
619+
type(linalg_state_type), optional, intent(out) :: err
620+
end subroutine stdlib_linalg_real_eig_${ri}$
621+
#:endif
622+
#:endfor
623+
end interface eig
624+
625+
! Eigenvalues of a square matrix
626+
interface eigvals
627+
!! version: experimental
628+
!!
629+
!! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for square matrix \( A \).
630+
!! ([Specification](../page/specs/stdlib_linalg.html#eigvals-eigenvalues-of-a-square-matrix))
631+
!!
632+
!!### Summary
633+
!! Function interface for computing the eigenvalues of a square matrix.
634+
!!
635+
!!### Description
636+
!!
637+
!! This interface provides functions for returning the eigenvalues of a general square matrix.
638+
!! Supported data types include `real` and `complex`, and no assumption is made on the matrix structure.
639+
!! An `error stop` is thrown in case of failure; otherwise, error information can be returned
640+
!! as an optional `type(linalg_state_type)` output flag.
641+
!!
642+
!!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`.
643+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
644+
!!
645+
#:for rk,rt,ri in RC_KINDS_TYPES
646+
#:if rk!="xdp"
647+
module function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda)
648+
!! Return an array of eigenvalues of matrix A.
649+
!> Input matrix A[m,n]
650+
${rt}$, intent(in), target :: a(:,:)
651+
!> [optional] state return flag. On error if not requested, the code will stop
652+
type(linalg_state_type), intent(out) :: err
653+
!> Array of singular values
654+
complex(${rk}$), allocatable :: lambda(:)
655+
end function stdlib_linalg_eigvals_${ri}$
656+
657+
module function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda)
658+
!! Return an array of eigenvalues of matrix A.
659+
!> Input matrix A[m,n]
660+
${rt}$, intent(in), target :: a(:,:)
661+
!> Array of singular values
662+
complex(${rk}$), allocatable :: lambda(:)
663+
end function stdlib_linalg_eigvals_noerr_${ri}$
664+
#:endif
665+
#:endfor
666+
end interface eigvals
667+
668+
! Eigendecomposition of a real symmetric or complex hermitian matrix
669+
interface eigh
670+
!! version: experimental
671+
!!
672+
!! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for a real symmetric
673+
!! \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix.
674+
!! ([Specification](../page/specs/stdlib_linalg.html#eigh-eigenvalues-and-eigenvectors-of-a-real-symmetric-or-complex-hermitian-square-matrix))
675+
!!
676+
!!### Summary
677+
!! Subroutine interface for computing eigenvalues and eigenvectors of a real symmetric or complex Hermitian square matrix.
678+
!!
679+
!!### Description
680+
!!
681+
!! This interface provides methods for computing the eigenvalues, and optionally eigenvectors,
682+
!! of a real symmetric or complex Hermitian square matrix. Supported data types include `real` and `complex`.
683+
!! The matrix must be symmetric (if `real`) or Hermitian (if `complex`). Only the lower or upper
684+
!! half of the matrix is accessed, and the user can select which using the optional `upper_a`
685+
!! flag (default: use lower half). The vectors are orthogonal, and may be returned as columns of an optional
686+
!! matrix `vectors` with the same kind and size as `A`.
687+
!! Preallocated space for both eigenvalues `lambda` and the eigenvector matrix must be user-provided.
688+
!!
689+
!!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`.
690+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
691+
!!
692+
#:for rk,rt,ri in RC_KINDS_TYPES
693+
#:if rk!="xdp"
694+
module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err)
695+
!! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda`
696+
!! of eigenvalues, and optionally right or left eigenvectors.
697+
!> Input matrix A[m,n]
698+
${rt}$, intent(inout), target :: a(:,:)
699+
!> Array of eigenvalues
700+
real(${rk}$), intent(out) :: lambda(:)
701+
!> The columns of vectors contain the orthonormal eigenvectors of A
702+
${rt}$, optional, intent(out), target :: vectors(:,:)
703+
!> [optional] Can A data be overwritten and destroyed?
704+
logical(lk), optional, intent(in) :: overwrite_a
705+
!> [optional] Should the upper/lower half of A be used? Default: lower
706+
logical(lk), optional, intent(in) :: upper_a
707+
!> [optional] state return flag. On error if not requested, the code will stop
708+
type(linalg_state_type), optional, intent(out) :: err
709+
end subroutine stdlib_linalg_eigh_${ri}$
710+
#:endif
711+
#:endfor
712+
end interface eigh
713+
714+
! Eigenvalues of a real symmetric or complex hermitian matrix
715+
interface eigvalsh
716+
!! version: experimental
717+
!!
718+
!! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for a real
719+
!! symmetric \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix.
720+
!! ([Specification](../page/specs/stdlib_linalg.html#eigvalsh-eigenvalues-of-a-real-symmetric-or-complex-hermitian-square-matrix))
721+
!!
722+
!!### Summary
723+
!! Function interface for computing the eigenvalues of a real symmetric or complex hermitian square matrix.
724+
!!
725+
!!### Description
726+
!!
727+
!! This interface provides functions for returning the eigenvalues of a real symmetric or complex Hermitian
728+
!! square matrix. Supported data types include `real` and `complex`. The matrix must be symmetric
729+
!! (if `real`) or Hermitian (if `complex`). Only the lower or upper half of the matrix is accessed,
730+
!! and the user can select which using the optional `upper_a` flag (default: use lower half).
731+
!! An `error stop` is thrown in case of failure; otherwise, error information can be returned
732+
!! as an optional `type(linalg_state_type)` output flag.
733+
!!
734+
!!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`.
735+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
736+
!!
737+
#:for rk,rt,ri in RC_KINDS_TYPES
738+
#:if rk!="xdp"
739+
module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda)
740+
!! Return an array of eigenvalues of real symmetric / complex hermitian A
741+
!> Input matrix A[m,n]
742+
${rt}$, intent(in), target :: a(:,:)
743+
!> [optional] Should the upper/lower half of A be used? Default: lower
744+
logical(lk), optional, intent(in) :: upper_a
745+
!> [optional] state return flag. On error if not requested, the code will stop
746+
type(linalg_state_type), intent(out) :: err
747+
!> Array of singular values
748+
real(${rk}$), allocatable :: lambda(:)
749+
end function stdlib_linalg_eigvalsh_${ri}$
750+
751+
module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda)
752+
!! Return an array of eigenvalues of real symmetric / complex hermitian A
753+
!> Input matrix A[m,n]
754+
${rt}$, intent(in), target :: a(:,:)
755+
!> [optional] Should the upper/lower half of A be used? Default: lower
756+
logical(lk), optional, intent(in) :: upper_a
757+
!> Array of singular values
758+
real(${rk}$), allocatable :: lambda(:)
759+
end function stdlib_linalg_eigvalsh_noerr_${ri}$
760+
#:endif
761+
#:endfor
762+
end interface eigvalsh
763+
764+
557765
! Singular value decomposition
558766
interface svd
559767
!! version: experimental

‎src/stdlib_linalg_eigenvalues.fypp

+572
Large diffs are not rendered by default.

‎test/linalg/CMakeLists.txt

+2
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@ set(
22
fppFiles
33
"test_linalg.fypp"
44
"test_blas_lapack.fypp"
5+
"test_linalg_eigenvalues.fypp"
56
"test_linalg_solve.fypp"
67
"test_linalg_lstsq.fypp"
78
"test_linalg_determinant.fypp"
@@ -12,6 +13,7 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
1213

1314
ADDTEST(linalg)
1415
ADDTEST(linalg_determinant)
16+
ADDTEST(linalg_eigenvalues)
1517
ADDTEST(linalg_matrix_property_checks)
1618
ADDTEST(linalg_solve)
1719
ADDTEST(linalg_lstsq)

‎test/linalg/test_linalg_eigenvalues.fypp

+210
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,210 @@
1+
#:include "common.fypp"
2+
! Test eigenvalues and eigendecompositions
3+
module test_linalg_eigenvalues
4+
use stdlib_linalg_constants
5+
use stdlib_linalg_state
6+
use stdlib_linalg, only: eig, eigh, eigvals, eigvalsh, diag
7+
use testdrive, only: error_type, check, new_unittest, unittest_type
8+
9+
implicit none (type,external)
10+
private
11+
12+
public :: test_eig_eigh
13+
14+
contains
15+
16+
!> SVD tests
17+
subroutine test_eig_eigh(tests)
18+
!> Collection of tests
19+
type(unittest_type), allocatable, intent(out) :: tests(:)
20+
21+
allocate(tests(0))
22+
23+
#:for rk,rt,ri in REAL_KINDS_TYPES
24+
#:if rk!="xdp"
25+
tests = [tests,new_unittest("test_eig_real_${ri}$",test_eig_real_${ri}$), &
26+
new_unittest("test_eigh_real_${ri}$",test_eigh_real_${ri}$)]
27+
#:endif
28+
#: endfor
29+
30+
#:for ck,ct,ci in CMPLX_KINDS_TYPES
31+
#:if ck!="xdp"
32+
tests = [tests,new_unittest("test_eig_complex_${ci}$",test_eig_complex_${ci}$)]
33+
#:endif
34+
#: endfor
35+
36+
end subroutine test_eig_eigh
37+
38+
!> Simple real matrix eigenvalues
39+
#:for rk,rt,ri in REAL_KINDS_TYPES
40+
#:if rk!="xdp"
41+
subroutine test_eig_real_${ri}$(error)
42+
type(error_type), allocatable, intent(out) :: error
43+
44+
!> Reference solution
45+
real(${rk}$), parameter :: zero = 0.0_${rk}$
46+
real(${rk}$), parameter :: two = 2.0_${rk}$
47+
real(${rk}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${rk}$
48+
real(${rk}$), parameter :: tol = sqrt(epsilon(zero))
49+
50+
!> Local variables
51+
type(linalg_state_type) :: state
52+
${rt}$ :: A(3,3),B(2,2)
53+
complex(${rk}$) :: lambda(3),Bvec(2,2),Bres(2,2)
54+
55+
!> Matrix with real eigenvalues
56+
A = reshape([1,0,0, &
57+
0,2,0, &
58+
0,0,3],[3,3])
59+
60+
call eig(A,lambda,err=state)
61+
62+
call check(error,state%ok(),state%print())
63+
if (allocated(error)) return
64+
65+
call check(error, all(aimag(lambda)==zero.and.real(lambda,kind=${rk}$)==[1,2,3]),'expected results')
66+
if (allocated(error)) return
67+
68+
!> Matrix with complex eigenvalues
69+
B = transpose(reshape([1, -1, &
70+
1, 1],[2,2]))
71+
72+
!> Expected right eigenvectors
73+
Bres(1,1:2) = sqrt2o2
74+
Bres(2,1) = cmplx(zero,-sqrt2o2,kind=${rk}$)
75+
Bres(2,2) = cmplx(zero,+sqrt2o2,kind=${rk}$)
76+
77+
call eig(B,lambda,right=Bvec,err=state)
78+
79+
call check(error,state%ok(),state%print())
80+
if (allocated(error)) return
81+
82+
call check(error, all(abs(Bres-Bvec)<=tol),'expected results')
83+
if (allocated(error)) return
84+
85+
end subroutine test_eig_real_${ri}$
86+
87+
! Symmetric matrix eigenvalues
88+
subroutine test_eigh_real_${ri}$(error)
89+
type(error_type), allocatable, intent(out) :: error
90+
91+
!> Reference solution
92+
real(${rk}$), parameter :: zero = 0.0_${rk}$
93+
real(${rk}$), parameter :: tol = sqrt(epsilon(zero))
94+
real(${rk}$), parameter :: A(4,4) = reshape([6,3,1,5, &
95+
3,0,5,1, &
96+
1,5,6,2, &
97+
5,1,2,2],[4,4])
98+
99+
!> Local variables
100+
real(${rk}$) :: Amat(4,4),lambda(4),vect(4,4),Av(4,4),lv(4,4)
101+
type(linalg_state_type) :: state
102+
103+
Amat = A
104+
105+
call eigh(Amat,lambda,vect,err=state)
106+
107+
Av = matmul(A,vect)
108+
lv = matmul(vect,diag(lambda))
109+
110+
call check(error,state%ok(),state%print())
111+
if (allocated(error)) return
112+
113+
call check(error, all(abs(Av-lv)<=tol*abs(Av)),'expected results')
114+
if (allocated(error)) return
115+
116+
!> Test functional versions: no state interface
117+
lambda = eigvalsh(Amat)
118+
119+
!> State interface
120+
lambda = eigvalsh(Amat,err=state)
121+
122+
call check(error,state%ok(),state%print())
123+
if (allocated(error)) return
124+
125+
!> Functional version, lower A
126+
Amat = A
127+
lambda = eigvalsh(Amat,upper_a=.false.,err=state)
128+
129+
call check(error,state%ok(),state%print())
130+
if (allocated(error)) return
131+
132+
end subroutine test_eigh_real_${ri}$
133+
134+
#:endif
135+
#:endfor
136+
137+
!> Simple complex matrix eigenvalues
138+
#:for ck,ct,ci in CMPLX_KINDS_TYPES
139+
#:if ck!="xdp"
140+
subroutine test_eig_complex_${ci}$(error)
141+
type(error_type), allocatable, intent(out) :: error
142+
143+
!> Reference solution
144+
real(${ck}$), parameter :: zero = 0.0_${ck}$
145+
real(${ck}$), parameter :: two = 2.0_${ck}$
146+
real(${ck}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${ck}$
147+
real(${ck}$), parameter :: tol = sqrt(epsilon(zero))
148+
${ct}$, parameter :: cone = (1.0_${ck}$,0.0_${ck}$)
149+
${ct}$, parameter :: cimg = (0.0_${ck}$,1.0_${ck}$)
150+
${ct}$, parameter :: czero = (0.0_${ck}$,0.0_${ck}$)
151+
152+
!> Local vaciables
153+
type(linalg_state_type) :: state
154+
${ct}$ :: A(2,2),lambda(2),Avec(2,2),Ares(2,2),lres(2)
155+
156+
!> Matcix with real eigenvalues
157+
A = transpose(reshape([ cone, cimg, &
158+
-cimg, cone], [2,2]))
159+
160+
call eig(A,lambda,right=Avec,err=state)
161+
162+
!> Expected eigenvalues and eigenvectors
163+
lres(1) = two
164+
lres(2) = zero
165+
166+
!> Eigenvectors may vary: do not use for error
167+
Ares(1,1) = cmplx(zero,sqrt2o2,kind=${ck}$)
168+
Ares(1,2) = cmplx(sqrt2o2,zero,kind=${ck}$)
169+
Ares(2,1) = cmplx(sqrt2o2,zero,kind=${ck}$)
170+
Ares(2,2) = cmplx(zero,sqrt2o2,kind=${ck}$)
171+
172+
call check(error,state%ok(),state%print())
173+
if (allocated(error)) return
174+
175+
call check(error, all(abs(lambda-lres)<=tol), 'results match expected')
176+
if (allocated(error)) return
177+
178+
end subroutine test_eig_complex_${ci}$
179+
180+
#:endif
181+
#:endfor
182+
183+
184+
end module test_linalg_eigenvalues
185+
186+
program test_eigenvalues
187+
use, intrinsic :: iso_fortran_env, only : error_unit
188+
use testdrive, only : run_testsuite, new_testsuite, testsuite_type
189+
use test_linalg_eigenvalues, only : test_eig_eigh
190+
implicit none
191+
integer :: stat, is
192+
type(testsuite_type), allocatable :: testsuites(:)
193+
character(len=*), parameter :: fmt = '("#", *(1x, a))'
194+
195+
stat = 0
196+
197+
testsuites = [ &
198+
new_testsuite("linalg_eigenvalues", test_eig_eigh) &
199+
]
200+
201+
do is = 1, size(testsuites)
202+
write(error_unit, fmt) "Testing:", testsuites(is)%name
203+
call run_testsuite(testsuites(is)%collect, error_unit, stat)
204+
end do
205+
206+
if (stat > 0) then
207+
write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
208+
error stop
209+
end if
210+
end program test_eigenvalues

0 commit comments

Comments
 (0)
Please sign in to comment.