From e6aa0f50e95eefe96cd5e6ade2c50e42fa158cdc Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 5 Jun 2024 09:19:15 +0200
Subject: [PATCH] Option to output `real` eigenvalues only

---
 doc/specs/stdlib_linalg.md         |  3 ++-
 src/stdlib_linalg.fypp             | 21 +++++++++++++++
 src/stdlib_linalg_eigenvalues.fypp | 42 +++++++++++++++++++++++++++++-
 3 files changed, 64 insertions(+), 2 deletions(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index cc3ee2778..aa50aab51 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -923,7 +923,7 @@ The solver is based on LAPACK's `*GEEV` backends.
 
 `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` rank-1 array of the same kind as `a`, containing the eigenvalues. It is an `intent(out)` argument. 
+`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.
 
@@ -937,6 +937,7 @@ The solver is based on LAPACK's `*GEEV` backends.
 
 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
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 90b920395..2e8d9ef8a 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -597,6 +597,27 @@ module stdlib_linalg
     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
diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp
index ae7ac482c..37f40ff8e 100644
--- a/src/stdlib_linalg_eigenvalues.fypp
+++ b/src/stdlib_linalg_eigenvalues.fypp
@@ -315,6 +315,47 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
 
      end subroutine stdlib_linalg_eig_${ri}$
      
+     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}$     
+     
      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]
@@ -535,5 +576,4 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
      #:endif
      #:endfor
 
-
 end submodule stdlib_linalg_eigenvalues