diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index a8fd633c7..4e8fe3b19 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -687,6 +687,8 @@ Expert (`Pure`) interface: `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. +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + ### Return value 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`. {!example/linalg/example_determinant2.f90!} ``` +## `eig` - Eigenvalues and Eigenvectors of a Square Matrix + +### Status + +Experimental + +### Description + +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. + +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 \). +Both `left` and `right` are rank-2 arrays, where eigenvectors are stored as columns. +The solver is based on LAPACK's `*GEEV` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):eig(interface)]] `(a, lambda [, right] [,left] [,overwrite_a] [,err])` + +### Arguments + +`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. + +`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. + +`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. + +`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. + +`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. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Raises `LINALG_ERROR` if the calculation did not converge. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_eig.f90!} +``` + +## `eigh` - Eigenvalues and Eigenvectors of a Real symmetric or Complex Hermitian Square Matrix + +### Status + +Experimental + +### Description + +This subroutine computes the solution to the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), +where \( A \) is a square, full-rank, `real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix. + +Result array `lambda` returns the `real` eigenvalues of \( A \). The user can request the orthogonal eigenvectors +to be returned: on output `vectors` may contain the matrix of eigenvectors, returned as a column. + +Normally, only the lower triangular part of \( A \) is accessed. On input, `logical` flag `upper_a` +allows the user to request what triangular part of the matrix should be used. + +The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):eigh(interface)]] `(a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])` + +### Arguments + +`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. + +`lambda`: Shall be a `complex` rank-1 array of the same precision as `a`, containing the eigenvalues. It is an `intent(out)` argument. + +`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. + +`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. + +`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. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Raises `LINALG_ERROR` if the calculation did not converge. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_eigh.f90!} +``` + +## `eigvals` - Eigenvalues of a Square Matrix + +### Status + +Experimental + +### Description + +This function returns the eigenvalues to matrix \( A \): a square, full-rank, `real` or `complex` matrix. +The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \). + +Result array `lambda` is `complex`, and returns the eigenvalues of \( A \). +The solver is based on LAPACK's `*GEEV` backends. + +### Syntax + +`lambda = ` [[stdlib_linalg(module):eigvals(interface)]] `(a, [,err])` + +### Arguments + +`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns a `complex` array containing the eigenvalues of `a`. + +Raises `LINALG_ERROR` if the calculation did not converge. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + + +### Example + +```fortran +{!example/linalg/example_eigvals.f90!} +``` + +## `eigvalsh` - Eigenvalues of a Real Symmetric or Complex Hermitian Square Matrix + +### Status + +Experimental + +### Description + +This function returns the eigenvalues to matrix \( A \): a where \( A \) is a square, full-rank, +`real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix. +The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \). + +Result array `lambda` is `real`, and returns the eigenvalues of \( A \). +The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. + +### Syntax + +`lambda = ` [[stdlib_linalg(module):eigvalsh(interface)]] `(a, [, upper_a] [,err])` + +### Arguments + +`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. + +`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. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns a `real` array containing the eigenvalues of `a`. + +Raises `LINALG_ERROR` if the calculation did not converge. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_eigvalsh.f90!} +``` + ## `svd` - Compute the singular value decomposition of a rank-2 array (matrix). ### Status @@ -989,4 +1165,3 @@ Exceptions trigger an `error stop`, unless argument `err` is present. ```fortran {!example/linalg/example_svdvals.f90!} ``` - diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index e9ad22a03..3606f0a7e 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -13,6 +13,10 @@ ADD_EXAMPLE(is_square) ADD_EXAMPLE(is_symmetric) ADD_EXAMPLE(is_triangular) ADD_EXAMPLE(outer_product) +ADD_EXAMPLE(eig) +ADD_EXAMPLE(eigh) +ADD_EXAMPLE(eigvals) +ADD_EXAMPLE(eigvalsh) ADD_EXAMPLE(trace) ADD_EXAMPLE(state1) ADD_EXAMPLE(state2) diff --git a/example/linalg/example_eig.f90 b/example/linalg/example_eig.f90 new file mode 100644 index 000000000..ca7e6728b --- /dev/null +++ b/example/linalg/example_eig.f90 @@ -0,0 +1,26 @@ +! Eigendecomposition of a real square matrix +program example_eig + use stdlib_linalg, only: eig + implicit none + + integer :: i + real, allocatable :: A(:,:) + complex, allocatable :: lambda(:),vectors(:,:) + + ! Decomposition of this square matrix + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 2, 4], & + [1, 3, 5], & + [2, 3, 4] ], [3,3] )) + + ! Get eigenvalues and right eigenvectors + allocate(lambda(3),vectors(3,3)) + + call eig(A, lambda, right=vectors) + + do i=1,3 + print *, 'eigenvalue ',i,': ',lambda(i) + print *, 'eigenvector ',i,': ',vectors(:,i) + end do + +end program example_eig diff --git a/example/linalg/example_eigh.f90 b/example/linalg/example_eigh.f90 new file mode 100644 index 000000000..7f19946ae --- /dev/null +++ b/example/linalg/example_eigh.f90 @@ -0,0 +1,38 @@ +! Eigendecomposition of a real symmetric matrix +program example_eigh + use stdlib_linalg, only: eigh + implicit none + + integer :: i + real, allocatable :: A(:,:),lambda(:),vectors(:,:) + complex, allocatable :: cA(:,:),cvectors(:,:) + + ! Decomposition of this symmetric matrix + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 1, 4], & + [1, 3, 5], & + [4, 5, 4] ], [3,3] )) + + ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors + allocate(lambda(3),vectors(3,3)) + call eigh(A, lambda, vectors=vectors) + + print *, 'Real matrix' + do i=1,3 + print *, 'eigenvalue ',i,': ',lambda(i) + print *, 'eigenvector ',i,': ',vectors(:,i) + end do + + ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors + cA = A + + allocate(cvectors(3,3)) + call eigh(cA, lambda, vectors=cvectors) + + print *, 'Complex matrix' + do i=1,3 + print *, 'eigenvalue ',i,': ',lambda(i) + print *, 'eigenvector ',i,': ',cvectors(:,i) + end do + +end program example_eigh diff --git a/example/linalg/example_eigvals.f90 b/example/linalg/example_eigvals.f90 new file mode 100644 index 000000000..7a8727e75 --- /dev/null +++ b/example/linalg/example_eigvals.f90 @@ -0,0 +1,24 @@ +! Eigenvalues of a general real / complex matrix +program example_eigvals + use stdlib_linalg, only: eigvals + implicit none + + integer :: i + real, allocatable :: A(:,:),lambda(:) + complex, allocatable :: cA(:,:),clambda(:) + + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 8, 4], & + [1, 3, 5], & + [9, 5,-2] ], [3,3] )) + + ! Note: real symmetric matrix + lambda = eigvals(A) + print *, 'Real matrix eigenvalues: ',lambda + + ! Complex general matrix + cA = cmplx(A, -2*A) + clambda = eigvals(cA) + print *, 'Complex matrix eigenvalues: ',clambda + +end program example_eigvals diff --git a/example/linalg/example_eigvalsh.f90 b/example/linalg/example_eigvalsh.f90 new file mode 100644 index 000000000..de672b0df --- /dev/null +++ b/example/linalg/example_eigvalsh.f90 @@ -0,0 +1,25 @@ +! Eigenvalues of a real symmetric / complex hermitian matrix +program example_eigvalsh + use stdlib_linalg, only: eigvalsh + implicit none + + integer :: i + real, allocatable :: A(:,:),lambda(:) + complex, allocatable :: cA(:,:) + + ! Decomposition of this symmetric matrix + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 1, 4], & + [1, 3, 5], & + [4, 5, 4] ], [3,3] )) + + ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors + lambda = eigvalsh(A) + print *, 'Symmetric matrix eigenvalues: ',lambda + + ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors + cA = A + lambda = eigvalsh(cA) + print *, 'Hermitian matrix eigenvalues: ',lambda + +end program example_eigvalsh diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5e27927b1..fa1831e98 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,7 +27,8 @@ set(fppFiles stdlib_linalg_outer_product.fypp stdlib_linalg_kronecker.fypp stdlib_linalg_cross_product.fypp - stdlib_linalg_solve.fypp + stdlib_linalg_eigenvalues.fypp + stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp stdlib_linalg_state.fypp stdlib_linalg_svd.fypp diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 10cd16650..aa4365b97 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -19,6 +19,10 @@ module stdlib_linalg public :: det public :: operator(.det.) public :: diag + public :: eig + public :: eigh + public :: eigvals + public :: eigvalsh public :: eye public :: lstsq public :: lstsq_space @@ -554,6 +558,210 @@ module stdlib_linalg #:endfor end interface + ! Eigendecomposition of a square matrix: eigenvalues, and optionally eigenvectors + interface eig + !! version: experimental + !! + !! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for square matrix \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#eig-eigenvalues-and-eigenvectors-of-a-square-matrix)) + !! + !!### Summary + !! Subroutine interface for computing eigenvalues and eigenvectors of a square matrix. + !! + !!### Description + !! + !! This interface provides methods for computing the eigenvalues, and optionally eigenvectors, + !! of a general square matrix. Supported data types include `real` and `complex`, and no assumption is + !! made on the matrix structure. The user may request either left, right, or both + !! eigenvectors to be returned. They are returned as columns of a square matrix with the same size as `A`. + !! Preallocated space for both eigenvalues `lambda` and the eigenvector matrices must be user-provided. + !! + !!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, + !! and optionally right or left eigenvectors. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of eigenvalues + complex(${rk}$), intent(out) :: lambda(:) + !> The columns of RIGHT contain the right eigenvectors of A + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> The columns of LEFT contain the left eigenvectors of A + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_eig_${ri}$ + #:endif + #:endfor + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, + !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had + !! non-trivial imaginary parts. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of real eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of RIGHT contain the right eigenvectors of A + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> The columns of LEFT contain the left eigenvectors of A + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_real_eig_${ri}$ + #:endif + #:endfor + end interface eig + + ! Eigenvalues of a square matrix + interface eigvals + !! version: experimental + !! + !! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for square matrix \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#eigvals-eigenvalues-of-a-square-matrix)) + !! + !!### Summary + !! Function interface for computing the eigenvalues of a square matrix. + !! + !!### Description + !! + !! This interface provides functions for returning the eigenvalues of a general square matrix. + !! Supported data types include `real` and `complex`, and no assumption is made on the matrix structure. + !! An `error stop` is thrown in case of failure; otherwise, error information can be returned + !! as an optional `type(linalg_state_type)` output flag. + !! + !!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda) + !! Return an array of eigenvalues of matrix A. + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out) :: err + !> Array of singular values + complex(${rk}$), allocatable :: lambda(:) + end function stdlib_linalg_eigvals_${ri}$ + + module function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda) + !! Return an array of eigenvalues of matrix A. + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Array of singular values + complex(${rk}$), allocatable :: lambda(:) + end function stdlib_linalg_eigvals_noerr_${ri}$ + #:endif + #:endfor + end interface eigvals + + ! Eigendecomposition of a real symmetric or complex hermitian matrix + interface eigh + !! version: experimental + !! + !! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for a real symmetric + !! \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix. + !! ([Specification](../page/specs/stdlib_linalg.html#eigh-eigenvalues-and-eigenvectors-of-a-real-symmetric-or-complex-hermitian-square-matrix)) + !! + !!### Summary + !! Subroutine interface for computing eigenvalues and eigenvectors of a real symmetric or complex Hermitian square matrix. + !! + !!### Description + !! + !! This interface provides methods for computing the eigenvalues, and optionally eigenvectors, + !! of a real symmetric or complex Hermitian square matrix. Supported data types include `real` and `complex`. + !! The matrix must be symmetric (if `real`) or Hermitian (if `complex`). Only the lower or upper + !! half of the matrix is accessed, and the user can select which using the optional `upper_a` + !! flag (default: use lower half). The vectors are orthogonal, and may be returned as columns of an optional + !! matrix `vectors` with the same kind and size as `A`. + !! Preallocated space for both eigenvalues `lambda` and the eigenvector matrix must be user-provided. + !! + !!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) + !! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` + !! of eigenvalues, and optionally right or left eigenvectors. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of vectors contain the orthonormal eigenvectors of A + ${rt}$, optional, intent(out), target :: vectors(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Should the upper/lower half of A be used? Default: lower + logical(lk), optional, intent(in) :: upper_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_eigh_${ri}$ + #:endif + #:endfor + end interface eigh + + ! Eigenvalues of a real symmetric or complex hermitian matrix + interface eigvalsh + !! version: experimental + !! + !! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for a real + !! symmetric \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix. + !! ([Specification](../page/specs/stdlib_linalg.html#eigvalsh-eigenvalues-of-a-real-symmetric-or-complex-hermitian-square-matrix)) + !! + !!### Summary + !! Function interface for computing the eigenvalues of a real symmetric or complex hermitian square matrix. + !! + !!### Description + !! + !! This interface provides functions for returning the eigenvalues of a real symmetric or complex Hermitian + !! square matrix. Supported data types include `real` and `complex`. The matrix must be symmetric + !! (if `real`) or Hermitian (if `complex`). Only the lower or upper half of the matrix is accessed, + !! and the user can select which using the optional `upper_a` flag (default: use lower half). + !! An `error stop` is thrown in case of failure; otherwise, error information can be returned + !! as an optional `type(linalg_state_type)` output flag. + !! + !!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) + !! Return an array of eigenvalues of real symmetric / complex hermitian A + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Should the upper/lower half of A be used? Default: lower + logical(lk), optional, intent(in) :: upper_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out) :: err + !> Array of singular values + real(${rk}$), allocatable :: lambda(:) + end function stdlib_linalg_eigvalsh_${ri}$ + + module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda) + !! Return an array of eigenvalues of real symmetric / complex hermitian A + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Should the upper/lower half of A be used? Default: lower + logical(lk), optional, intent(in) :: upper_a + !> Array of singular values + real(${rk}$), allocatable :: lambda(:) + end function stdlib_linalg_eigvalsh_noerr_${ri}$ + #:endif + #:endfor + end interface eigvalsh + + ! Singular value decomposition interface svd !! version: experimental diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp new file mode 100644 index 000000000..2c8ee8783 --- /dev/null +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -0,0 +1,572 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +submodule (stdlib_linalg) stdlib_linalg_eigenvalues +!! Compute eigenvalues and eigenvectors + use stdlib_linalg_constants + use stdlib_linalg_lapack, only: geev, heev, syev + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS + implicit none(type,external) + + character(*), parameter :: this = 'eigenvalues' + + contains + + !> Request for eigenvector calculation + elemental character function eigenvectors_task(required) + logical(lk), intent(in) :: required + eigenvectors_task = merge('V','N',required) + end function eigenvectors_task + + !> Request for symmetry side (default: lower) + elemental character function symmetric_triangle_task(upper) + logical(lk), optional, intent(in) :: upper + symmetric_triangle_task = 'L' + if (present(upper)) symmetric_triangle_task = merge('U','L',upper) + end function symmetric_triangle_task + + !> Process GEEV output flags + elemental subroutine handle_geev_info(err,info,m,n) + !> Error handler + type(linalg_state_type), intent(inout) :: err + !> GEEV return flag + integer(ilp), intent(in) :: info + !> Input matrix size + integer(ilp), intent(in) :: m,n + + select case (info) + case (0) + ! Success! + err%state = LINALG_SUCCESS + case (-1) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.') + case (-2) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.') + case (-5,-3) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + case (-9) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.') + case (-11) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.') + case (-13) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.') + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.') + end select + + end subroutine handle_geev_info + + !> Process SYEV/HEEV output flags + elemental subroutine handle_heev_info(err,info,m,n) + !> Error handler + type(linalg_state_type), intent(inout) :: err + !> SYEV/HEEV return flag + integer(ilp), intent(in) :: info + !> Input matrix size + integer(ilp), intent(in) :: m,n + + select case (info) + case (0) + ! Success! + err%state = LINALG_SUCCESS + case (-1) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.') + case (-2) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.') + case (-5,-3) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + case (-8) + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.') + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.') + end select + + end subroutine handle_heev_info + + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + + module function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda) + !! Return an array of eigenvalues of matrix A. + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [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 :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eig_${ri}$(amat,lambda,overwrite_a=.false.,err=err) + + end function stdlib_linalg_eigvals_${ri}$ + + module function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda) + !! Return an array of eigenvalues of matrix A. + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Array of eigenvalues + complex(${rk}$), allocatable :: lambda(:) + + !> Create + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eig_${ri}$(amat,lambda,overwrite_a=.false.) + + end function stdlib_linalg_eigvals_noerr_${ri}$ + + module subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, + !! and optionally right or left eigenvectors. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of eigenvalues + complex(${rk}$), intent(out) :: lambda(:) + !> [optional] RIGHT eigenvectors of A (as columns) + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> [optional] LEFT eigenvectors of A (as columns) + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + !> Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: m,n,lda,ldu,ldv,info,k,lwork,lrwork,neig + logical(lk) :: copy_a + character :: task_u,task_v + ${rt}$, target :: work_dummy(1),u_dummy(1,1),v_dummy(1,1) + ${rt}$, allocatable :: work(:) + real(${rk}$), allocatable :: rwork(:) + ${rt}$, pointer :: amat(:,:),lreal(:),limag(:),umat(:,:),vmat(:,:) + + !> Matrix size + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + neig = size(lambda,kind=ilp) + lda = m + + if (k<=0 .or. m/=n) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& + 'invalid or matrix size a=',[m,n],', must be nonempty square.') + call linalg_error_handling(err0,err) + return + elseif (neig right + #:else + ! For a real matrix, GEEV returns real arrays. + ! Allocate temporary reals, will be converted to complex vectors at the end. + allocate(vmat(n,n)) + #:endif + + if (size(vmat,2,kind=ilp) left + #:else + ! For a real matrix, GEEV returns real arrays. + ! Allocate temporary reals, will be converted to complex vectors at the end. + allocate(umat(n,n)) + #:endif + + if (size(umat,2,kind=ilp) Prepare working storage + lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) + allocate(work(lwork)) + + !> Compute eigensystem + call geev(task_u,task_v,n,amat,lda,& + #{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# & + umat,ldu,vmat,ldv,& + work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info) + call handle_geev_info(err0,info,m,n) + + endif + + ! Finalize storage and process output flag + #:if not rt.startswith('complex') + lambda(:n) = cmplx(lreal(:n),limag(:n),kind=${rk}$) + + ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, + ! geev returns reals as: + ! u(j) = VL(:,j) + i*VL(:,j+1) and + ! u(j+1) = VL(:,j) - i*VL(:,j+1). + ! Convert these to complex numbers here. + if (present(right)) call assign_real_eigenvectors_${rk}$(n,lambda,vmat,right) + if (present(left)) call assign_real_eigenvectors_${rk}$(n,lambda,umat,left) + #:endif + + endif get_geev + + if (copy_a) deallocate(amat) + #:if not rt.startswith('complex') + if (present(right)) deallocate(vmat) + if (present(left)) deallocate(umat) + #:endif + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_eig_${ri}$ + + module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) + !! Return an array of eigenvalues of real symmetric / complex hermitian A + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Should the upper/lower half of A be used? Default: lower + logical(lk), optional, intent(in) :: upper_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out) :: err + !> Array of eigenvalues + real(${rk}$), allocatable :: lambda(:) + + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.,err=err) + + end function stdlib_linalg_eigvalsh_${ri}$ + + module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda) + !! Return an array of eigenvalues of real symmetric / complex hermitian A + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Should the upper/lower half of A be used? Default: use lower + logical(lk), optional, intent(in) :: upper_a + !> Array of eigenvalues + real(${rk}$), allocatable :: lambda(:) + + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.) + + end function stdlib_linalg_eigvalsh_noerr_${ri}$ + + module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) + !! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` + !! of eigenvalues, and optionally right or left eigenvectors. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of vectors contain the orthonormal eigenvectors of A + ${rt}$, optional, intent(out), target :: vectors(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Should the upper/lower half of A be used? Default: use lower + logical(lk), optional, intent(in) :: upper_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + !> Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: m,n,lda,info,k,lwork,neig + logical(lk) :: copy_a + character :: triangle,task + ${rt}$, target :: work_dummy(1) + ${rt}$, allocatable :: work(:) + real(${rk}$), allocatable :: rwork(:) + ${rt}$, pointer :: amat(:,:) + + !> Matrix size + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + neig = size(lambda,kind=ilp) + + if (k<=0 .or. m/=n) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & + ', must be non-empty square.') + call linalg_error_handling(err0,err) + return + elseif (neig= n=',n) + call linalg_error_handling(err0,err) + return + endif + + ! Check if input A can be overwritten + copy_a = .true._lk + if (present(vectors)) then + ! No need to copy A anyways + copy_a = .false. + elseif (present(overwrite_a)) then + copy_a = .not.overwrite_a + endif + + ! Should we use the upper or lower half of the matrix? + triangle = symmetric_triangle_task(upper_a) + + ! Request for eigenvectors + task = eigenvectors_task(present(vectors)) + + if (present(vectors)) then + + ! Check size + if (any(shape(vectors,kind=ilp) Prepare working storage + lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) + allocate(work(lwork)) + + !> Compute eigensystem + #:if rt.startswith('complex') + call heev(task,triangle,n,amat,lda,lambda,work,lwork,rwork,info) + #:else + call syev(task,triangle,n,amat,lda,lambda,work,lwork,info) + #:endif + call handle_heev_info(err0,info,m,n) + + endif + + ! Finalize storage and process output flag + if (copy_a) deallocate(amat) + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_eigh_${ri}$ + + #:endif + #:endfor + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + pure subroutine assign_real_eigenvectors_${rk}$(n,lambda,lmat,out_mat) + !! GEEV for real matrices returns complex eigenvalues in real arrays, where two consecutive + !! reals at [j,j+1] locations represent the real and imaginary parts of two complex conjugate + !! eigenvalues. Convert them to complex here, following the GEEV logic. + !> Problem size + integer(ilp), intent(in) :: n + !> Array of eigenvalues + complex(${rk}$), intent(in) :: lambda(:) + !> Real matrix as returned by geev + ${rt}$, intent(in) :: lmat(:,:) + !> Complex matrix as returned by eig + complex(${rk}$), intent(out) :: out_mat(:,:) + + integer(ilp) :: i,j + + ! Copy matrix + do concurrent(i=1:n,j=1:n) + out_mat(i,j) = lmat(i,j) + end do + + ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, + ! geev returns them as reals as: + ! u(j) = VL(:,j) + i*VL(:,j+1) and + ! u(j+1) = VL(:,j) - i*VL(:,j+1). + ! Convert these to complex numbers here. + do j=1,n-1 + if (lambda(j)==conjg(lambda(j+1))) then + out_mat(:, j) = cmplx(lmat(:,j), lmat(:,j+1),kind=${rk}$) + out_mat(:,j+1) = cmplx(lmat(:,j),-lmat(:,j+1),kind=${rk}$) + endif + end do + + end subroutine assign_real_eigenvectors_${rk}$ + + module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, + !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had + !! non-trivial imaginary parts. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of real eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of RIGHT contain the right eigenvectors of A + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> The columns of LEFT contain the left eigenvectors of A + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + type(linalg_state_type) :: err0 + integer(ilp) :: n + complex(${rk}$), allocatable :: clambda(:) + real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$) + real(${rk}$), parameter :: atol = tiny(0.0_${rk}$) + + n = size(lambda,dim=1,kind=ilp) + allocate(clambda(n)) + + call stdlib_linalg_eig_${ri}$(a,clambda,right,left,overwrite_a,err0) + + ! Check that no eigenvalues have meaningful imaginary part + if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR, & + 'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(clambda))) + endif + + ! Return real components only + lambda(:n) = real(clambda,kind=${rk}$) + + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_real_eig_${ri}$ + + #:endif + #:endfor + +end submodule stdlib_linalg_eigenvalues diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index eacfff53c..cff60532d 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -2,6 +2,7 @@ set( fppFiles "test_linalg.fypp" "test_blas_lapack.fypp" + "test_linalg_eigenvalues.fypp" "test_linalg_solve.fypp" "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" @@ -12,6 +13,7 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_determinant) +ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) diff --git a/test/linalg/test_linalg_eigenvalues.fypp b/test/linalg/test_linalg_eigenvalues.fypp new file mode 100644 index 000000000..ecec0be3b --- /dev/null +++ b/test/linalg/test_linalg_eigenvalues.fypp @@ -0,0 +1,210 @@ +#:include "common.fypp" +! Test eigenvalues and eigendecompositions +module test_linalg_eigenvalues + use stdlib_linalg_constants + use stdlib_linalg_state + use stdlib_linalg, only: eig, eigh, eigvals, eigvalsh, diag + use testdrive, only: error_type, check, new_unittest, unittest_type + + implicit none (type,external) + private + + public :: test_eig_eigh + + contains + + !> SVD tests + subroutine test_eig_eigh(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + 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_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}$)] + #:endif + #: 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 + + !> Reference solution + real(${rk}$), parameter :: zero = 0.0_${rk}$ + real(${rk}$), parameter :: two = 2.0_${rk}$ + real(${rk}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${rk}$ + real(${rk}$), parameter :: tol = sqrt(epsilon(zero)) + + !> Local variables + type(linalg_state_type) :: state + ${rt}$ :: A(3,3),B(2,2) + complex(${rk}$) :: lambda(3),Bvec(2,2),Bres(2,2) + + !> Matrix with real eigenvalues + A = reshape([1,0,0, & + 0,2,0, & + 0,0,3],[3,3]) + + call eig(A,lambda,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(aimag(lambda)==zero.and.real(lambda,kind=${rk}$)==[1,2,3]),'expected results') + if (allocated(error)) return + + !> Matrix with complex eigenvalues + B = transpose(reshape([1, -1, & + 1, 1],[2,2])) + + !> Expected right eigenvectors + Bres(1,1:2) = sqrt2o2 + Bres(2,1) = cmplx(zero,-sqrt2o2,kind=${rk}$) + Bres(2,2) = cmplx(zero,+sqrt2o2,kind=${rk}$) + + call eig(B,lambda,right=Bvec,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(Bres-Bvec)<=tol),'expected results') + if (allocated(error)) return + + end subroutine test_eig_real_${ri}$ + + ! Symmetric matrix eigenvalues + subroutine test_eigh_real_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + !> Reference solution + real(${rk}$), parameter :: zero = 0.0_${rk}$ + real(${rk}$), parameter :: tol = sqrt(epsilon(zero)) + real(${rk}$), parameter :: A(4,4) = reshape([6,3,1,5, & + 3,0,5,1, & + 1,5,6,2, & + 5,1,2,2],[4,4]) + + !> Local variables + real(${rk}$) :: Amat(4,4),lambda(4),vect(4,4),Av(4,4),lv(4,4) + type(linalg_state_type) :: state + + Amat = A + + call eigh(Amat,lambda,vect,err=state) + + Av = matmul(A,vect) + lv = matmul(vect,diag(lambda)) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(Av-lv)<=tol*abs(Av)),'expected results') + if (allocated(error)) return + + !> Test functional versions: no state interface + lambda = eigvalsh(Amat) + + !> State interface + lambda = eigvalsh(Amat,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + !> Functional version, lower A + Amat = A + lambda = eigvalsh(Amat,upper_a=.false.,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + end subroutine test_eigh_real_${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 + + !> Reference solution + real(${ck}$), parameter :: zero = 0.0_${ck}$ + real(${ck}$), parameter :: two = 2.0_${ck}$ + real(${ck}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${ck}$ + real(${ck}$), parameter :: tol = sqrt(epsilon(zero)) + ${ct}$, parameter :: cone = (1.0_${ck}$,0.0_${ck}$) + ${ct}$, parameter :: cimg = (0.0_${ck}$,1.0_${ck}$) + ${ct}$, parameter :: czero = (0.0_${ck}$,0.0_${ck}$) + + !> Local vaciables + type(linalg_state_type) :: state + ${ct}$ :: A(2,2),lambda(2),Avec(2,2),Ares(2,2),lres(2) + + !> Matcix with real eigenvalues + A = transpose(reshape([ cone, cimg, & + -cimg, cone], [2,2])) + + call eig(A,lambda,right=Avec,err=state) + + !> Expected eigenvalues and eigenvectors + lres(1) = two + lres(2) = zero + + !> Eigenvectors may vary: do not use for error + Ares(1,1) = cmplx(zero,sqrt2o2,kind=${ck}$) + Ares(1,2) = cmplx(sqrt2o2,zero,kind=${ck}$) + Ares(2,1) = cmplx(sqrt2o2,zero,kind=${ck}$) + Ares(2,2) = cmplx(zero,sqrt2o2,kind=${ck}$) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(lambda-lres)<=tol), 'results match expected') + if (allocated(error)) return + + end subroutine test_eig_complex_${ci}$ + + #:endif + #:endfor + + +end module test_linalg_eigenvalues + +program test_eigenvalues + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_eigenvalues, only : test_eig_eigh + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_eigenvalues", test_eig_eigh) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_eigenvalues