linalg: Eigenvalues and Eigenvectors (fortran-lang#816)

Author: perazz

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!}
 ```
-

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)

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

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

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

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

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