linalg eig: generalized eigenvalue problem (#909)

Author: perazz

doc/specs/stdlib_linalg.md

@@ -1091,20 +1091,30 @@ Stable
 
 ### 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.
+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, or to the generalized eigenproblem \( A \cdot \bar{v} - \lambda \cdot B \cdot \bar{v} \), 
+where \( B \) is a square matrix with the same type, kind and size as \( A \).
 
 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.
+The solver is based on LAPACK's `*GEEV` (standard eigenproblem) and `*GGEV` (generalized eigenproblem) backends.
 
 ### Syntax
 
+For the standard eigenproblem: 
+
 `call ` [[stdlib_linalg(module):eig(interface)]] `(a, lambda [, right] [,left] [,overwrite_a] [,err])`
 
+For the generalized eigenproblem: 
+
+`call ` [[stdlib_linalg(module):eig(interface)]] `(a, b, lambda [, right] [, left] [, overwrite_a] [, overwrite_b] [, 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. 
 
+`b`: `real` or `complex` square array containing the second coefficient matrix. If `overwrite_b=.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.
@@ -1113,6 +1123,8 @@ The solver is based on LAPACK's `*GEEV` backends.
 
 `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.
 
+`overwrite_b` (optional): Shall be an input logical flag. If `.true.`, input matrix `b` 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
@@ -1185,28 +1197,41 @@ Stable
 
 ### 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} \).
+This function computes the eigenvalues for either a standard or generalized eigenproblem:
 
-Result array `lambda` is `complex`, and returns the eigenvalues of \( A \). 
-The solver is based on LAPACK's `*GEEV` backends.
+- **Standard eigenproblem**: \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank `real` or `complex` matrix.
+- **Generalized eigenproblem**: \( A \cdot \bar{v} - \lambda \cdot B \cdot \bar{v} \), where \( B \) is a square matrix with the same type and kind as \( A \).
+
+The eigenvalues are stored in the result array `lambda`, which is `complex` (even for real input matrices).  
+The solver uses LAPACK's `*GEEV` and `*GGEV` backends for the standard and generalized problems, respectively.
 
 ### Syntax
 
-`lambda = ` [[stdlib_linalg(module):eigvals(interface)]] `(a, [,err])`
+For the standard eigenproblem:
+
+`lambda = ` [[stdlib_linalg(module):eigvals(interface)]] `(a [, err])`
+
+For the generalized eigenproblem:
+
+`lambda = ` [[stdlib_linalg(module):eigvals(interface)]] `(a, b [, err])`
 
 ### Arguments
 
-`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument.
+`a`:  
+Shall be 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.
+`b` (optional):  
+Shall be a `real` or `complex` square array containing the second coefficient matrix for the generalized problem. It is an `intent(in)` argument.
 
-### Return value
+`err` (optional):  
+Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
 
-Returns a `complex` array containing the eigenvalues of `a`. 
+### Return Value
 
-Raises `LINALG_ERROR` if the calculation did not converge.
-Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
+Returns a `complex` rank-1 array containing the eigenvalues of the problem.  
+
+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

example/linalg/CMakeLists.txt

@@ -21,8 +21,10 @@ ADD_EXAMPLE(pseudoinverse)
 ADD_EXAMPLE(outer_product)
 ADD_EXAMPLE(eig)
 ADD_EXAMPLE(eigh)
+ADD_EXAMPLE(eig_generalized)
 ADD_EXAMPLE(eigvals)
 ADD_EXAMPLE(eigvalsh)
+ADD_EXAMPLE(eigvals_generalized)
 ADD_EXAMPLE(trace)
 ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)

example/linalg/example_eig_generalized.f90

@@ -0,0 +1,30 @@
+! Eigendecomposition of a real square matrix for the generalized eigenproblem
+program example_eig_generalized
+  use stdlib_linalg, only: eig, eye
+  implicit none
+
+  integer :: i
+  real, allocatable :: A(:,:), B(:,:)
+  complex, allocatable :: lambda(:), vectors(:,:)
+
+  ! Matrices for the generalized eigenproblem: A * v = lambda * B * v
+  ! NB Fortran is column-major -> transpose input
+  A = transpose(reshape([ [2, 2, 4], &
+                          [1, 3, 5], &
+                          [2, 3, 4] ], [3,3]))
+
+  B = eye(3)
+
+  ! Allocate eigenvalues and right eigenvectors
+  allocate(lambda(3), vectors(3,3))
+
+  ! Get eigenvalues and right eigenvectors for the generalized problem
+  call eig(A, B, lambda, right=vectors)
+
+  do i = 1, 3
+     print *, 'Eigenvalue  ', i, ': ', lambda(i)
+     print *, 'Eigenvector ', i, ': ', vectors(:,i)
+  end do
+
+end program example_eig_generalized
+

example/linalg/example_eigvals_generalized.f90

@@ -0,0 +1,26 @@
+! Eigenvalues of a general real/complex matrix for the generalized eigenproblem
+program example_eigvals_generalized
+  use stdlib_linalg, only: eigvals, eye
+  implicit none
+
+  real, allocatable :: A(:,:), B(:,:), lambda(:)
+  complex, allocatable :: cA(:,:), cB(:,:), clambda(:)
+
+  ! NB Fortran is column-major -> transpose input
+  A = transpose(reshape([ [2, 8, 4], &
+                          [1, 3, 5], &
+                          [9, 5,-2] ], [3,3]))
+
+  B = eye(3)
+
+  ! Real generalized eigenproblem
+  lambda = eigvals(A, B)
+  print *, 'Real generalized matrix eigenvalues: ', lambda
+
+  ! Complex generalized eigenproblem
+  cA = cmplx(A, -2*A)
+  cB = cmplx(B, 0.5*B)
+  clambda = eigvals(cA, cB)
+  print *, 'Complex generalized matrix eigenvalues: ', clambda
+
+end program example_eigvals_generalized

src/stdlib_linalg.fypp

@@ -5,6 +5,9 @@
 #:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]]
 #:set RHS_EMPTY  = [emptyranksuffix(r) for r in [1,2]]
 #:set ALL_RHS    = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY))
+#:set EIG_PROBLEM  = ["standard", "generalized"]
+#:set EIG_FUNCTION = ["geev","ggev"]
+#:set EIG_PROBLEM_LIST = list(zip(EIG_PROBLEM, EIG_FUNCTION))
 module stdlib_linalg
   !!Provides a support for various linear algebra procedures
   !! ([Specification](../page/specs/stdlib_linalg.html))
@@ -1084,12 +1087,17 @@ module stdlib_linalg
      !!@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)
+    #:for ep,ei in EIG_PROBLEM_LIST
+    module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
+                                                      overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,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(:,:)
+         #:if ei=='ggev'
+         !> Generalized problem matrix B[n,n]
+         ${rt}$, intent(inout), target :: b(:,:)
+         #:endif         
          !> Array of eigenvalues
          complex(${rk}$), intent(out) :: lambda(:)
          !> The columns of RIGHT contain the right eigenvectors of A
@@ -1098,19 +1106,25 @@ module stdlib_linalg
          complex(${rk}$), optional, intent(out), target :: left(:,:)
          !> [optional] Can A data be overwritten and destroyed?
          logical(lk), optional, intent(in) :: overwrite_a
+         #:if ei=='ggev'
+         !> [optional] Can B data be overwritten and destroyed? (default: no)
+         logical(lk), optional, intent(in) :: overwrite_b
+         #:endif         
          !> [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)
+    end subroutine stdlib_linalg_eig_${ep}$_${ri}$
+
+    module subroutine stdlib_linalg_real_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
+                                                           overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,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(:,:)
+         #:if ei=='ggev'
+         !> Generalized problem matrix B[n,n]
+         ${rt}$, intent(inout), target :: b(:,:)  
+         #:endif
          !> Array of real eigenvalues
          real(${rk}$), intent(out) :: lambda(:)
          !> The columns of RIGHT contain the right eigenvectors of A
@@ -1119,11 +1133,15 @@ module stdlib_linalg
          complex(${rk}$), optional, intent(out), target :: left(:,:)
          !> [optional] Can A data be overwritten and destroyed?
          logical(lk), optional, intent(in) :: overwrite_a
+         #:if ei=='ggev'
+         !> [optional] Can B data be overwritten and destroyed? (default: no)
+         logical(lk), optional, intent(in) :: overwrite_b
+         #:endif         
          !> [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
+    end subroutine stdlib_linalg_real_eig_${ep}$_${ri}$
     #:endfor    
+    #:endfor
   end interface eig
 
   ! Eigenvalues of a square matrix
@@ -1147,25 +1165,33 @@ module stdlib_linalg
      !!@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)
+    #:for ep,ei in EIG_PROBLEM_LIST
+    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), target :: a(:,:)
+         ${rt}$, intent(in), dimension(:,:), target :: a 
+         #:if ei=='ggev'
+         !> Generalized problem matrix B[n,n]
+         ${rt}$, intent(inout), dimension(:,:), 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 singular values
          complex(${rk}$), allocatable :: lambda(:)        
-    end function stdlib_linalg_eigvals_${ri}$
+    end function stdlib_linalg_eigvals_${ep}$_${ri}$
 
-    module function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda)
+    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), target :: a(:,:)
+         ${rt}$, intent(in), dimension(:,:), target :: a 
+         #:if ei=='ggev'
+         !> Generalized problem matrix B[n,n]
+         ${rt}$, intent(inout), dimension(:,:), target :: b         
+         #:endif                  
          !> Array of singular values
          complex(${rk}$), allocatable :: lambda(:)
-    end function stdlib_linalg_eigvals_noerr_${ri}$
-    #:endif
+    end function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$
+    #:endfor
     #:endfor
   end interface eigvals
 
@@ -1194,7 +1220,6 @@ module stdlib_linalg
      !!@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.        
@@ -1211,7 +1236,6 @@ module stdlib_linalg
          !> [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
 
@@ -1239,7 +1263,6 @@ module stdlib_linalg
      !!@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]
@@ -1261,7 +1284,6 @@ module stdlib_linalg
          !> Array of singular values
          real(${rk}$), allocatable :: lambda(:)
     end function stdlib_linalg_eigvalsh_noerr_${ri}$
-    #:endif
     #:endfor                
   end interface eigvalsh