From 0222737c64436e38dff559ef2247451772ff1271 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Sat, 16 Nov 2024 09:38:01 +0100
Subject: [PATCH 01/24] implement Schur

---
 src/stdlib_linalg.fypp       |  74 ++++++++++++++++++
 src/stdlib_linalg_schur.fypp | 140 +++++++++++++++++++++++++++++++++++
 2 files changed, 214 insertions(+)
 create mode 100644 src/stdlib_linalg_schur.fypp

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index ad36ad677..a3ec82a46 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -598,6 +598,80 @@ module stdlib_linalg
     #:endfor
   end interface qr_space
  
+  ! Schur decomposition of rank-2 array A
+  interface schur
+    !! version: experimental
+    !!
+    !! Computes the Schur decomposition of matrix \( A = Z T Z^H \).
+    !! ([Specification](../page/specs/stdlib_linalg.html#schur-compute-the-schur-decomposition-of-a-matrix))
+    !!
+    !!### Summary
+    !! Compute the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is
+    !! orthonormal/unitary and \( T \) is upper-triangular or quasi-upper-triangular. Matrix \( A \) has size `[m,m]`.
+    !!
+    !!### Description
+    !! 
+    !! This interface provides methods for computing the Schur decomposition of a matrix.
+    !! Supported data types include `real` and `complex`. If a pre-allocated workspace is provided, no internal 
+    !! memory allocations take place when using this interface.
+    !!
+    !! The output matrix \( T \) is upper-triangular for `complex` input matrices and quasi-upper-triangular
+    !! for `real` input matrices (with possible \( 2 \times 2 \) blocks on the diagonal).
+    !! The user can optionally request sorting of eigenvalues based on conditions, or a custom sorting function.
+    !!
+    !!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`). Sorting options 
+    !! are implemented using LAPACK's eigenvalue sorting mechanism.
+    !! 
+    #:for rk,rt,ri in RC_KINDS_TYPES
+      pure module subroutine stdlib_linalg_${ri}$_schur(a,t,z,sdim,output,overwrite_a,sort,storage,err)
+         !> Input matrix a[m,m]
+         ${rt}$, intent(inout), target :: a(:,:)
+         !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
+         ${rt}$, intent(out), contiguous, target :: t(:,:)
+         !> Unitary/orthonormal transformation matrix Z
+         ${rt}$, intent(out), contiguous, target :: z(:,:)
+         !> [optional] Number of eigenvalues satisfying the sort condition
+         integer(ilp), optional, intent(out) :: sdim
+         !> [optional] Output type: 'real' or 'complex'
+         character(*), optional, intent(in) :: output
+         !> [optional] Can A data be overwritten and destroyed?
+         logical(lk), optional, intent(in) :: overwrite_a
+         !> [optional] Sorting criterion: callable or predefined values ('lhp', 'rhp', 'iuc', 'ouc', or None)
+         class(*), optional, intent(in) :: sort
+         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
+         ${rt}$, intent(out), optional, target :: storage(:)
+         !> [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_${ri}$_schur
+    #:endfor
+  end interface schur
+
+  ! Return the working array space required by the Schur decomposition solver
+  interface schur_space
+    !! version: experimental
+    !!
+    !! Computes the working array space required by the Schur decomposition solver
+    !! ([Specification](../page/specs/stdlib_linalg.html#schur-space-compute-internal-working-space-requirements-for-the-schur-decomposition))
+    !!
+    !!### Description
+    !! 
+    !! This interface returns the size of the `real` or `complex` working storage required by the 
+    !! Schur decomposition solver. The working size only depends on the kind (`real` or `complex`) and size of
+    !! the matrix being decomposed. Storage size can be used to pre-allocate a working array in case several 
+    !! repeated Schur decompositions of same-size matrices are sought. If pre-allocated working arrays 
+    !! are provided, no internal allocations will take place during the decomposition.
+    !!     
+    #:for rk,rt,ri in RC_KINDS_TYPES
+      pure module subroutine get_schur_${ri}$_workspace(a,lwork,err)
+         !> Input matrix a[m,m]
+         ${rt}$, intent(in), target :: a(:,:)
+         !> Minimum workspace size for the decomposition operation
+         integer(ilp), intent(out) :: lwork
+         !> State return flag. Returns an error if the query failed
+         type(linalg_state_type), optional, intent(out) :: err
+      end subroutine get_schur_${ri}$_workspace
+    #:endfor
+  end interface schur_space
 
   interface det
     !! version: experimental 
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
new file mode 100644
index 000000000..0fb17b738
--- /dev/null
+++ b/src/stdlib_linalg_schur.fypp
@@ -0,0 +1,140 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_linalg) stdlib_linalg_schur
+    use stdlib_linalg_constants
+    use stdlib_linalg_lapack, only: gees
+    use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+        LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
+    implicit none(type,external)
+
+    character(*), parameter :: this = 'schur'
+
+    contains
+
+    elemental subroutine handle_gees_info(info, m, sort, err)
+        integer(ilp), intent(in) :: info, m
+        logical, intent(in) :: sort
+        type(linalg_state_type), intent(out) :: err
+
+        ! Process GEES output
+        select case (info)
+        case (0)
+            ! Success
+        case (-1)
+            err = linalg_state_type(this, LINALG_VALUE_ERROR, 'invalid matrix size m=', m)
+        case default
+            if (sort .and. info > 0) then
+                err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'sorting eigenvalues failed at index ', info)
+            else
+                err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
+            end if
+        end select
+    end subroutine handle_gees_info
+
+    #:for rk, rt, ri in RC_KINDS_TYPES
+
+    ! Schur decomposition subroutine
+    pure module subroutine stdlib_linalg_${ri}$_schur(a, t, z, lwork, overwrite_a, sort, err)
+        !> Input matrix a[m, m]
+        ${rt}$, intent(inout), target :: a(:,:)
+        !> Schur form of the matrix A
+        ${rt}$, intent(out), contiguous, target :: t(:,:)
+        !> Unitary/orthogonal matrix Z
+        ${rt}$, intent(out), contiguous, target :: z(:,:)
+        !> Workspace size (optional)
+        integer(ilp), optional, intent(inout) :: lwork
+        !> Overwrite input matrix A? (optional)
+        logical(lk), optional, intent(in) :: overwrite_a
+        !> Sorting eigenvalues? (optional)
+        logical(lk), optional, intent(in) :: sort
+        !> 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, lda, info, liwork
+        logical(lk) :: overwrite_a_
+        logical, pointer :: bwork(:)
+        integer(ilp), allocatable :: iwork(:)
+        ${rt}$, pointer :: amat(:,:), tau(:), work(:)
+        ${rt}$ :: rwork_dummy(1)  ! Dummy for real/complex cases
+        ${rt}$, allocatable :: tmat(:,:), zmat(:,:)
+        character :: jobz
+
+        ! Problem size
+        m = size(a, 1, kind=ilp)
+        lda = size(a, 1, kind=ilp)
+
+        ! Validate dimensions
+        if (size(a, 1, kind=ilp) /= size(a, 2, kind=ilp)) then
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix A must be square: a=', [size(a,1), size(a,2)])
+            call linalg_error_handling(err0, err)
+            return
+        end if
+
+        ! Set default values
+        overwrite_a_ = .false._lk
+        if (present(overwrite_a)) overwrite_a_ = overwrite_a
+
+        ! Job type based on sorting request
+        if (present(sort) .and. sort) then
+            jobz = 'S'  ! Compute and sort eigenvalues
+        else
+            jobz = 'N'  ! Compute Schur decomposition without sorting
+        end if
+
+        ! Prepare storage
+        allocate(tmat(m, m), source=0.0_${rk}$)
+        allocate(zmat(m, m), source=0.0_${rk}$)
+
+        if (overwrite_a_) then
+            amat => a
+        else
+            allocate(amat(m, m), source=a)
+        end if
+
+        ! Allocate workspace
+        liwork = -1
+        if (present(lwork)) then
+            allocate(work(lwork))
+        else
+            allocate(work(1))  ! Temporary workspace for querying size
+        end if
+
+        ! Workspace query
+        call #{if rt.startswith('complex')}# zgees #{else}# gees #{endif}# &
+            (jobz, 'N', nullify(bwork), m, amat, lda, tau, zmat, lda, work, liwork, iwork, info)
+        call handle_gees_info(info, m, .false._lk, err0)
+        if (err0%error()) then
+            call linalg_error_handling(err0, err)
+            return
+        end if
+
+        ! Update workspace size
+        if (.not.present(lwork)) then
+            liwork = ceiling(real(work(1), kind=${rk}$), kind=ilp)
+            deallocate(work)
+            allocate(work(liwork))
+        end if
+
+        ! Compute Schur decomposition
+        call #{if rt.startswith('complex')}# zgees #{else}# gees #{endif}# &
+            (jobz, 'N', nullify(bwork), m, amat, lda, tau, zmat, lda, work, liwork, iwork, info)
+        call handle_gees_info(info, m, present(sort) .and. sort, err0)
+        if (err0%error()) then
+            call linalg_error_handling(err0, err)
+            return
+        end if
+
+        ! Output results
+        t = amat
+        z = zmat
+
+        if (.not.overwrite_a_) deallocate(amat)
+        if (.not.present(lwork)) deallocate(work)
+    end subroutine stdlib_linalg_${ri}$_schur
+
+    #:endfor
+
+end submodule stdlib_linalg_schur
+

From 04d843f14388ca10230257c19a00ff8893941716 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 01:52:46 -0600
Subject: [PATCH 02/24] handle `*GEES` output

---
 src/stdlib_linalg_schur.fypp | 44 +++++++++++++++++++++++++++---------
 1 file changed, 33 insertions(+), 11 deletions(-)

diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 0fb17b738..683e96c55 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -11,24 +11,46 @@ submodule (stdlib_linalg) stdlib_linalg_schur
 
     contains
 
-    elemental subroutine handle_gees_info(info, m, sort, err)
-        integer(ilp), intent(in) :: info, m
-        logical, intent(in) :: sort
-        type(linalg_state_type), intent(out) :: err
+    !> Wrapper function to handle GEES error codes
+    elemental subroutine handle_gees_info(info, m, n, ldvs, err)
+        integer(ilp), intent(in) :: info, m, n, ldvs
+        type(linalg_state), intent(out) :: err
 
         ! Process GEES output
         select case (info)
-        case (0)
+        case (0_ilp)
             ! Success
-        case (-1)
-            err = linalg_state_type(this, LINALG_VALUE_ERROR, 'invalid matrix size m=', m)
-        case default
-            if (sort .and. info > 0) then
-                err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'sorting eigenvalues failed at index ', info)
+        case (-1_ilp)
+            ! Vector not wanted, but task is wrong
+            err = linalg_state(this, LINALG_INTERNAL_ERROR,'Invalid Schur vector task request')
+        case (-2_ilp)
+            ! Vector not wanted, but task is wrong
+            err = linalg_state(this, LINALG_INTERNAL_ERROR,'Invalid sorting task request')            
+        case (-4_ilp,-6_ilp)
+            ! Vector not wanted, but task is wrong
+            err = linalg_state(this, LINALG_VALUE_ERROR,'Invalid/non-square input matrix size:',[m,n])  
+        case (-11_ilp)
+            err = linalg_state(this, LINALG_VALUE_ERROR,'Schur vector matrix has insufficient size',[ldvs,n])
+        case (-13_ilp)
+            err = linalg_state(this, LINALG_INTERNAL_ERROR,'Insufficient working storage size')
+        case (1_ilp:)
+            
+            if (info==n+2) then 
+                err = linalg_state(this, LINALG_ERROR, 'Ill-conditioned problem: could not sort eigenvalues')
+            elseif (info==n+1) then 
+                err = linalg_state(this, LINALG_ERROR, 'Some selected eigenvalues lost property due to sorting')                
+            elseif (info==n) then 
+                err = linalg_state(this, LINALG_ERROR, 'Convergence failure: no converged eigenvalues')
             else
-                err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
+                err = linalg_state(this, LINALG_ERROR, 'Convergence failure; converged range is',[info,n])                            
             end if
+            
+        case default
+            
+            err = linalg_state(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
+
         end select
+        
     end subroutine handle_gees_info
 
     #:for rk, rt, ri in RC_KINDS_TYPES

From 908cab16aa22baea36aeae6b3d79dd4fefc57531 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 01:53:17 -0600
Subject: [PATCH 03/24] process `*GEES` tasks

---
 src/stdlib_linalg_schur.fypp | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 683e96c55..59a4804a8 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -8,9 +8,33 @@ submodule (stdlib_linalg) stdlib_linalg_schur
     implicit none(type,external)
 
     character(*), parameter :: this = 'schur'
+    
+    !> List of internal GEES tasks:     
+    !> No task request
+    character, parameter :: GEES_NOT            = 'N'
+    
+    !> Request Schur vectors to be computed
+    character, parameter :: GEES_WITH_VECTORS   = 'V'
+    
+    !> Request Schur vectors to be sorted
+    character, parameter :: GEES_SORTED_VECTORS = 'S' 
 
     contains
 
+    !> Wrapper function for Schur vectors request
+    elemental character function gees_vectors(wanted) 
+        !> Are Schur vectors wanted?
+        logical(lk), intent(in) :: wanted 
+        gees_vectors = merge(GEES_WITH_VECTORS,GEES_NOT,wanted)
+    end function gees_vectors
+    
+    !> Wrapper function for Schur vectors request
+    elemental character function gees_sort_eigs(sorted) 
+        !> Should the eigenvalues be sorted?
+        logical(lk), intent(in) :: sorted 
+        gees_sort_eigs = merge(GEES_SORTED_VECTORS,GEES_NOT,sorted)
+    end function gees_sort_eigs    
+
     !> Wrapper function to handle GEES error codes
     elemental subroutine handle_gees_info(info, m, n, ldvs, err)
         integer(ilp), intent(in) :: info, m, n, ldvs

From dc816880e46d0c9659fb651a4a967edb550e680c Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 01:55:42 -0600
Subject: [PATCH 04/24] workspace query `schur_space`

---
 src/stdlib_linalg.fypp       |  2 +-
 src/stdlib_linalg_schur.fypp | 50 ++++++++++++++++++++++++++++++++++--
 2 files changed, 49 insertions(+), 3 deletions(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index a3ec82a46..9b00040b8 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -662,7 +662,7 @@ module stdlib_linalg
     !! are provided, no internal allocations will take place during the decomposition.
     !!     
     #:for rk,rt,ri in RC_KINDS_TYPES
-      pure module subroutine get_schur_${ri}$_workspace(a,lwork,err)
+      module subroutine get_schur_${ri}$_workspace(a,lwork,err)
          !> Input matrix a[m,m]
          ${rt}$, intent(in), target :: a(:,:)
          !> Minimum workspace size for the decomposition operation
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 59a4804a8..11d3e9466 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -38,7 +38,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
     !> Wrapper function to handle GEES error codes
     elemental subroutine handle_gees_info(info, m, n, ldvs, err)
         integer(ilp), intent(in) :: info, m, n, ldvs
-        type(linalg_state), intent(out) :: err
+        type(linalg_state_type), intent(out) :: err
 
         ! Process GEES output
         select case (info)
@@ -77,7 +77,53 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         
     end subroutine handle_gees_info
 
-    #:for rk, rt, ri in RC_KINDS_TYPES
+    #:for rk, rt, ri in RC_KINDS_TYPES    
+    !> Workspace query
+    subroutine get_schur_${ri}$_workspace(a,lwork,err)
+        !> Input matrix a[m,m]
+        ${rt}$, intent(in), target :: a(:,:)
+        !> Minimum workspace size for the decomposition operation
+        integer(ilp), intent(out) :: lwork
+        !> State return flag. Returns an error if the query failed
+        type(linalg_state_type), optional, intent(out) :: err
+        
+        integer(ilp) :: m,n,sdim,info
+        character :: jobvs,sort
+        logical(lk) :: bwork_dummy(1)
+        ${rt}$, pointer :: amat(:,:)
+        real(${rk}$) :: rwork_dummy(1)
+        ${rt}$ :: wr_dummy(1),wi_dummy(1),vs_dummy(1,1),work_dummy(1)
+        type(linalg_state_type) :: err0
+        
+        !> Initialize problem
+        lwork = -1_ilp
+        m = size(a,1,kind=ilp)
+        n = size(a,2,kind=ilp)         
+        
+        !> Create a dummy intent(inout) argument
+        amat => a
+        
+        !> Select dummy task
+        jobvs = gees_vectors(.true.)
+        sort  = gees_sort_eigs(.false.)
+        sdim  = 0_ilp
+        
+        !> Get Schur workspace 
+        call gees(jobvs,sort,do_not_select,n,amat,m,sdim,wr_dummy,#{if rt.startswith('r')}#wi_dummy, #{endif}#&
+                  vs_dummy,m,work_dummy,lwork,#{if rt.startswith('c')}#rwork_dummy,#{endif}#bwork_dummy,info)
+        if (info==0) lwork = nint(real(work_dummy(1),kind=${rk}$),kind=ilp)
+        call handle_gees_info(info,m,n,m,err0)
+        call linalg_error_handling(err0,err)
+        
+        contains
+        
+            ! Interface to dummy select routine
+            pure logical(lk) function do_not_select(alpha#{if rt.startswith('r')}#r,alphai#{endif}#) 
+                ${rt}$, intent(in) :: alpha#{if rt.startswith('r')}#r,alphai#{endif}# 
+                do_not_select = .false.
+            end function do_not_select 
+        
+    end subroutine get_schur_${ri}$_workspace   
 
     ! Schur decomposition subroutine
     pure module subroutine stdlib_linalg_${ri}$_schur(a, t, z, lwork, overwrite_a, sort, err)

From 81e5583ea7138e9fa81276a288da0d396ebb80a4 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 01:56:59 -0600
Subject: [PATCH 05/24] schur: simplify interface, do not allocate where
 possible

---
 src/stdlib_linalg.fypp       |  16 +--
 src/stdlib_linalg_schur.fypp | 218 ++++++++++++++++++++++-------------
 2 files changed, 144 insertions(+), 90 deletions(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 9b00040b8..1e1a39a02 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -623,23 +623,17 @@ module stdlib_linalg
     !! are implemented using LAPACK's eigenvalue sorting mechanism.
     !! 
     #:for rk,rt,ri in RC_KINDS_TYPES
-      pure module subroutine stdlib_linalg_${ri}$_schur(a,t,z,sdim,output,overwrite_a,sort,storage,err)
+      module subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, storage, err)
          !> Input matrix a[m,m]
          ${rt}$, intent(inout), target :: a(:,:)
          !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
          ${rt}$, intent(out), contiguous, target :: t(:,:)
          !> Unitary/orthonormal transformation matrix Z
-         ${rt}$, intent(out), contiguous, target :: z(:,:)
-         !> [optional] Number of eigenvalues satisfying the sort condition
-         integer(ilp), optional, intent(out) :: sdim
-         !> [optional] Output type: 'real' or 'complex'
-         character(*), optional, intent(in) :: output
-         !> [optional] Can A data be overwritten and destroyed?
-         logical(lk), optional, intent(in) :: overwrite_a
-         !> [optional] Sorting criterion: callable or predefined values ('lhp', 'rhp', 'iuc', 'ouc', or None)
-         class(*), optional, intent(in) :: sort
+         ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
+         !> [optional] Output eigenvalues that appear on the diagonal of T
+         complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
          !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
-         ${rt}$, intent(out), optional, target :: storage(:)
+         ${rt}$, optional, intent(inout), target :: storage(:)
          !> [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_${ri}$_schur
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 11d3e9466..090137f5e 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -126,104 +126,164 @@ submodule (stdlib_linalg) stdlib_linalg_schur
     end subroutine get_schur_${ri}$_workspace   
 
     ! Schur decomposition subroutine
-    pure module subroutine stdlib_linalg_${ri}$_schur(a, t, z, lwork, overwrite_a, sort, err)
-        !> Input matrix a[m, m]
+    subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, storage, err)
+        !> Input matrix a[m,m]
         ${rt}$, intent(inout), target :: a(:,:)
-        !> Schur form of the matrix A
+        !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
         ${rt}$, intent(out), contiguous, target :: t(:,:)
-        !> Unitary/orthogonal matrix Z
-        ${rt}$, intent(out), contiguous, target :: z(:,:)
-        !> Workspace size (optional)
-        integer(ilp), optional, intent(inout) :: lwork
-        !> Overwrite input matrix A? (optional)
-        logical(lk), optional, intent(in) :: overwrite_a
-        !> Sorting eigenvalues? (optional)
-        logical(lk), optional, intent(in) :: sort
-        !> State return flag. On error if not requested, the code will stop
-        type(linalg_state_type), optional, intent(out) :: err
+        !> Unitary/orthonormal transformation matrix Z
+        ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
+        !> [optional] Output eigenvalues that appear on the diagonal of T
+        complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
+        !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
+        ${rt}$, optional, intent(inout), target :: storage(:)
+        !> [optional] State return flag. On error if not requested, the code will stop
+        type(linalg_state), optional, intent(out) :: err
 
         ! Local variables
-        type(linalg_state_type) :: err0
-        integer(ilp) :: m, lda, info, liwork
-        logical(lk) :: overwrite_a_
-        logical, pointer :: bwork(:)
-        integer(ilp), allocatable :: iwork(:)
-        ${rt}$, pointer :: amat(:,:), tau(:), work(:)
-        ${rt}$ :: rwork_dummy(1)  ! Dummy for real/complex cases
-        ${rt}$, allocatable :: tmat(:,:), zmat(:,:)
-        character :: jobz
-
+        integer(ilp) :: m,n,mt,nt,ldvs,nvs,lde,lwork,sdim,info
+        logical(lk), target :: bwork_dummy(1),local_eigs
+        logical(lk), pointer :: bwork(:)
+        real(${rk}$), allocatable :: rwork(:)
+        ${rt}$, target :: vs_dummy(1,1)
+        ${rt}$, pointer :: vs(:,:),work(:),eigs(:)#{if rt.startswith('r')}#,eigi(:)#{endif}#
+        character :: jobvs,sort
+        type(linalg_state) :: err0
+        
         ! Problem size
-        m = size(a, 1, kind=ilp)
-        lda = size(a, 1, kind=ilp)
-
+        m  = size(a, 1, kind=ilp)
+        n  = size(a, 2, kind=ilp)        
+        mt = size(t, 1, kind=ilp)
+        nt = size(t, 2, kind=ilp)
+        
         ! Validate dimensions
-        if (size(a, 1, kind=ilp) /= size(a, 2, kind=ilp)) then
-            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix A must be square: a=', [size(a,1), size(a,2)])
+        if (m/=n .or. m<=0 .or. n<=0) then
+            err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Matrix A must be square: size(a)=',[m,n])
             call linalg_error_handling(err0, err)
             return
-        end if
-
-        ! Set default values
-        overwrite_a_ = .false._lk
-        if (present(overwrite_a)) overwrite_a_ = overwrite_a
-
-        ! Job type based on sorting request
-        if (present(sort) .and. sort) then
-            jobz = 'S'  ! Compute and sort eigenvalues
+        end if    
+        if (mt/=nt .or. mt/=n .or. nt/=n) then
+            err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Matrix T must be square: size(T)=',[mt,nt], &
+                                                          'should be',[m,n])
+            call linalg_error_handling(err0, err)
+            return
+        end if                
+        
+        !> Copy data into the output array
+        t = a 
+        
+        !> SORTING: no sorting options are currently supported
+        sort  = gees_sort_eigs(.false.)
+        sdim  = 0_ilp     
+        
+        if (sort/=GEES_NOT) then 
+            
+           allocate(bwork(n),source=.false.)
+        
         else
-            jobz = 'N'  ! Compute Schur decomposition without sorting
-        end if
-
-        ! Prepare storage
-        allocate(tmat(m, m), source=0.0_${rk}$)
-        allocate(zmat(m, m), source=0.0_${rk}$)
-
-        if (overwrite_a_) then
-            amat => a
+            
+           bwork => bwork_dummy 
+            
+        end if        
+        
+        !> Schur vectors
+        jobvs = gees_vectors(present(z))        
+        if (present(z)) then 
+            vs => z
+            
+            ldvs = size(vs, 1, kind=ilp)
+            nvs  = size(vs, 2, kind=ilp)
+            
+            if (ldvs<n .or. nvs/=n) then
+                err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Schur vectors size=',[ldvs,nvs], &
+                                                              'should be n=',n)
+                goto 1
+            end if             
+            
         else
-            allocate(amat(m, m), source=a)
+            vs => vs_dummy
+            ldvs = size(vs, 1, kind=ilp)
+            nvs  = size(vs, 2, kind=ilp)
         end if
-
-        ! Allocate workspace
-        liwork = -1
-        if (present(lwork)) then
-            allocate(work(lwork))
+        
+        !> User or self-allocated storage
+        if (present(storage)) then 
+            
+            work => storage
+            lwork = size(work, 1, kind=ilp)
+            
         else
-            allocate(work(1))  ! Temporary workspace for querying size
+            
+            ! Query optimal workspace            
+            call get_schur_${ri}$_workspace(a,lwork,err0)
+            if (err0%error()) goto 1
+            allocate(work(lwork))
+            
         end if
-
-        ! Workspace query
-        call #{if rt.startswith('complex')}# zgees #{else}# gees #{endif}# &
-            (jobz, 'N', nullify(bwork), m, amat, lda, tau, zmat, lda, work, liwork, iwork, info)
-        call handle_gees_info(info, m, .false._lk, err0)
-        if (err0%error()) then
-            call linalg_error_handling(err0, err)
-            return
+        
+        !> User or self-allocated eigenvalue storage
+        if (present(eigvals)) then 
+            
+            lde = size(eigvals, 1, kind=ilp)
+            
+            #:if rt.startswith('c')
+            eigs => eigvals
+            local_eigs = .false.
+            #:else
+            allocate(eigs(n),eigi(n))
+            local_eigs = .true.
+            #:endif
+            
+        else
+            
+            allocate(eigs(n)#{if rt.startswith('r')}#,eigi(n)#{endif}#)
+            local_eigs = .true.
+            lde = n
+            
         end if
+        
+        #:if rt.startswith('c')
+        allocate(rwork(n))
+        #:endif
 
-        ! Update workspace size
-        if (.not.present(lwork)) then
-            liwork = ceiling(real(work(1), kind=${rk}$), kind=ilp)
-            deallocate(work)
-            allocate(work(liwork))
-        end if
+        if (lde<n) then
+            err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Insufficient eigenvalue array size=',lde, &
+                                                          'should be >=',n)
+            goto 2
+        end if  
 
         ! Compute Schur decomposition
-        call #{if rt.startswith('complex')}# zgees #{else}# gees #{endif}# &
-            (jobz, 'N', nullify(bwork), m, amat, lda, tau, zmat, lda, work, liwork, iwork, info)
-        call handle_gees_info(info, m, present(sort) .and. sort, err0)
-        if (err0%error()) then
-            call linalg_error_handling(err0, err)
-            return
-        end if
+        call gees(jobvs,sort,eig_select,nt,t,mt,sdim,eigs,#{if rt.startswith('r')}#eigi,#{endif}# &
+                  vs,ldvs,work,lwork,#{if rt.startswith('c')}#rwork,#{endif}#bwork,info)
+        call handle_gees_info(info,m,n,m,err0)
 
-        ! Output results
-        t = amat
-        z = zmat
+      2 eigenvalue_output: if (local_eigs) then 
+           #:if rt.startswith('r')
+           ! Build complex eigenvalues
+           eigvals = cmplx(eigs,eigi,kind=${rk}$)
+           deallocate(eigs,eigi)
+           #:else
+           deallocate(eigs)
+           #:endif
+        endif eigenvalue_output
+        if (.not.present(storage)) deallocate(work)
+      1 if (sort/=GEES_NOT) deallocate(bwork)
+        call linalg_error_handling(err0,err)
+        
+        contains
+        
+            ! Dummy select routine: currently, no sorting options are offered
+            pure logical(lk) function eig_select(alpha#{if rt.startswith('r')}#r,alphai#{endif}#)             
+                #:if rt.startswith('r')
+                ${rt}$, intent(in) :: alphar,alphai
+                complex(${rk}$)    :: alpha
+                alpha = cmplx(alphar,alphai,kind=${rk}$)
+                #:else
+                ${rt}$, intent(in) :: alpha
+                #:endif
+                eig_select = .false.
+            end function eig_select 
 
-        if (.not.overwrite_a_) deallocate(amat)
-        if (.not.present(lwork)) deallocate(work)
     end subroutine stdlib_linalg_${ri}$_schur
 
     #:endfor

From 95a4900fe987f885e723029ee9d6427afbe078c9 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 02:16:42 -0600
Subject: [PATCH 06/24] `overwrite_a` option

---
 src/stdlib_linalg.fypp       |  6 ++++--
 src/stdlib_linalg_schur.fypp | 35 ++++++++++++++++++++++++++++-------
 2 files changed, 32 insertions(+), 9 deletions(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 1e1a39a02..479eb61d1 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -623,7 +623,7 @@ module stdlib_linalg
     !! are implemented using LAPACK's eigenvalue sorting mechanism.
     !! 
     #:for rk,rt,ri in RC_KINDS_TYPES
-      module subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, storage, err)
+      module subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, overwrite_a, storage, err)
          !> Input matrix a[m,m]
          ${rt}$, intent(inout), target :: a(:,:)
          !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
@@ -632,7 +632,9 @@ module stdlib_linalg
          ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
          !> [optional] Output eigenvalues that appear on the diagonal of T
          complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
-         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
+         !> [optional] Can A data be overwritten and destroyed?
+         logical(lk), optional, intent(in) :: overwrite_a                 
+         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space         
          ${rt}$, optional, intent(inout), target :: storage(:)
          !> [optional] State return flag. On error if not requested, the code will stop
          type(linalg_state_type), optional, intent(out) :: err
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 090137f5e..e36b5351e 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -79,7 +79,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
 
     #:for rk, rt, ri in RC_KINDS_TYPES    
     !> Workspace query
-    subroutine get_schur_${ri}$_workspace(a,lwork,err)
+    module subroutine get_schur_${ri}$_workspace(a,lwork,err)
         !> Input matrix a[m,m]
         ${rt}$, intent(in), target :: a(:,:)
         !> Minimum workspace size for the decomposition operation
@@ -126,7 +126,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
     end subroutine get_schur_${ri}$_workspace   
 
     ! Schur decomposition subroutine
-    subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, storage, err)
+    module subroutine stdlib_linalg_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
         !> Input matrix a[m,m]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
@@ -137,11 +137,14 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
         ${rt}$, optional, intent(inout), target :: storage(:)
+        !> [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), optional, intent(out) :: err
 
         ! Local variables
         integer(ilp) :: m,n,mt,nt,ldvs,nvs,lde,lwork,sdim,info
+        logical(lk) :: overwrite_a_
         logical(lk), target :: bwork_dummy(1),local_eigs
         logical(lk), pointer :: bwork(:)
         real(${rk}$), allocatable :: rwork(:)
@@ -172,6 +175,13 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         !> Copy data into the output array
         t = a 
         
+        ! Can A be overwritten? By default, do not overwrite
+        if (present(overwrite_a)) then 
+            overwrite_a_ = overwrite_a .and. n>=2    
+        else
+            overwrite_a_ = .false._lk
+        endif
+        
         !> SORTING: no sorting options are currently supported
         sort  = gees_sort_eigs(.false.)
         sdim  = 0_ilp     
@@ -230,13 +240,26 @@ submodule (stdlib_linalg) stdlib_linalg_schur
             eigs => eigvals
             local_eigs = .false.
             #:else
-            allocate(eigs(n),eigi(n))
+            ! use A storage if possible
+            if (overwrite_a_) then                 
+                eigs => a(:,1)
+                eigi => a(:,2)
+            else
+                allocate(eigs(n),eigi(n))
+            end if            
             local_eigs = .true.
             #:endif
             
         else
             
-            allocate(eigs(n)#{if rt.startswith('r')}#,eigi(n)#{endif}#)
+            ! Use A storage if possible
+            if (overwrite_a_) then                 
+                eigs => a(:,1)
+                eigi => a(:,2)
+            else
+                allocate(eigs(n)#{if rt.startswith('r')}#,eigi(n)#{endif}#)
+            end if               
+            
             local_eigs = .true.
             lde = n
             
@@ -261,10 +284,8 @@ submodule (stdlib_linalg) stdlib_linalg_schur
            #:if rt.startswith('r')
            ! Build complex eigenvalues
            eigvals = cmplx(eigs,eigi,kind=${rk}$)
-           deallocate(eigs,eigi)
-           #:else
-           deallocate(eigs)
            #:endif
+           if (.not.overwrite_a_) deallocate(eigs#{if rt.startswith('r')}#,eigi#{endif}#)
         endif eigenvalue_output
         if (.not.present(storage)) deallocate(work)
       1 if (sort/=GEES_NOT) deallocate(bwork)

From ead68dde3e7866cab9e95fe324538ba68af21c13 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 02:58:29 -0600
Subject: [PATCH 07/24] fix, export interface

---
 src/CMakeLists.txt           |  1 +
 src/stdlib_linalg.fypp       |  2 ++
 src/stdlib_linalg_schur.fypp | 66 ++++++++++++++++--------------------
 3 files changed, 32 insertions(+), 37 deletions(-)

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ff9f39417..d6f521bc2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -36,6 +36,7 @@ set(fppFiles
     stdlib_linalg_state.fypp
     stdlib_linalg_svd.fypp 
     stdlib_linalg_cholesky.fypp
+    stdlib_linalg_schur.fypp
     stdlib_optval.fypp
     stdlib_selection.fypp
     stdlib_sorting.fypp
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 479eb61d1..f896eff57 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -44,6 +44,8 @@ module stdlib_linalg
   public :: cross_product
   public :: qr
   public :: qr_space
+  public :: schur
+  public :: schur_space
   public :: is_square
   public :: is_diagonal
   public :: is_symmetric
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index e36b5351e..26ca20850 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -46,32 +46,32 @@ submodule (stdlib_linalg) stdlib_linalg_schur
             ! Success
         case (-1_ilp)
             ! Vector not wanted, but task is wrong
-            err = linalg_state(this, LINALG_INTERNAL_ERROR,'Invalid Schur vector task request')
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid Schur vector task request')
         case (-2_ilp)
             ! Vector not wanted, but task is wrong
-            err = linalg_state(this, LINALG_INTERNAL_ERROR,'Invalid sorting task request')            
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid sorting task request')            
         case (-4_ilp,-6_ilp)
             ! Vector not wanted, but task is wrong
-            err = linalg_state(this, LINALG_VALUE_ERROR,'Invalid/non-square input matrix size:',[m,n])  
+            err = linalg_state_type(this, LINALG_VALUE_ERROR,'Invalid/non-square input matrix size:',[m,n])  
         case (-11_ilp)
-            err = linalg_state(this, LINALG_VALUE_ERROR,'Schur vector matrix has insufficient size',[ldvs,n])
+            err = linalg_state_type(this, LINALG_VALUE_ERROR,'Schur vector matrix has insufficient size',[ldvs,n])
         case (-13_ilp)
-            err = linalg_state(this, LINALG_INTERNAL_ERROR,'Insufficient working storage size')
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Insufficient working storage size')
         case (1_ilp:)
             
             if (info==n+2) then 
-                err = linalg_state(this, LINALG_ERROR, 'Ill-conditioned problem: could not sort eigenvalues')
+                err = linalg_state_type(this, LINALG_ERROR, 'Ill-conditioned problem: could not sort eigenvalues')
             elseif (info==n+1) then 
-                err = linalg_state(this, LINALG_ERROR, 'Some selected eigenvalues lost property due to sorting')                
+                err = linalg_state_type(this, LINALG_ERROR, 'Some selected eigenvalues lost property due to sorting')                
             elseif (info==n) then 
-                err = linalg_state(this, LINALG_ERROR, 'Convergence failure: no converged eigenvalues')
+                err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure: no converged eigenvalues')
             else
-                err = linalg_state(this, LINALG_ERROR, 'Convergence failure; converged range is',[info,n])                            
+                err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure; converged range is',[info,n])                            
             end if
             
         case default
             
-            err = linalg_state(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
 
         end select
         
@@ -140,7 +140,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         !> [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), optional, intent(out) :: err
+        type(linalg_state_type), optional, intent(out) :: err
 
         ! Local variables
         integer(ilp) :: m,n,mt,nt,ldvs,nvs,lde,lwork,sdim,info
@@ -151,7 +151,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         ${rt}$, target :: vs_dummy(1,1)
         ${rt}$, pointer :: vs(:,:),work(:),eigs(:)#{if rt.startswith('r')}#,eigi(:)#{endif}#
         character :: jobvs,sort
-        type(linalg_state) :: err0
+        type(linalg_state_type) :: err0
         
         ! Problem size
         m  = size(a, 1, kind=ilp)
@@ -161,12 +161,12 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         
         ! Validate dimensions
         if (m/=n .or. m<=0 .or. n<=0) then
-            err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Matrix A must be square: size(a)=',[m,n])
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix A must be square: size(a)=',[m,n])
             call linalg_error_handling(err0, err)
             return
         end if    
         if (mt/=nt .or. mt/=n .or. nt/=n) then
-            err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Matrix T must be square: size(T)=',[mt,nt], &
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix T must be square: size(T)=',[mt,nt], &
                                                           'should be',[m,n])
             call linalg_error_handling(err0, err)
             return
@@ -205,7 +205,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
             nvs  = size(vs, 2, kind=ilp)
             
             if (ldvs<n .or. nvs/=n) then
-                err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Schur vectors size=',[ldvs,nvs], &
+                err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Schur vectors size=',[ldvs,nvs], &
                                                               'should be n=',n)
                 goto 1
             end if             
@@ -232,45 +232,37 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         end if
         
         !> User or self-allocated eigenvalue storage
-        if (present(eigvals)) then 
-            
-            lde = size(eigvals, 1, kind=ilp)
-            
+        if (present(eigvals)) then             
+            lde = size(eigvals, 1, kind=ilp)            
             #:if rt.startswith('c')
             eigs => eigvals
             local_eigs = .false.
             #:else
-            ! use A storage if possible
-            if (overwrite_a_) then                 
-                eigs => a(:,1)
-                eigi => a(:,2)
-            else
-                allocate(eigs(n),eigi(n))
-            end if            
             local_eigs = .true.
-            #:endif
-            
-        else
-            
+            #:endif            
+        else            
+            local_eigs = .true.
+            lde = n            
+        end if
+
+        if (local_eigs) then 
             ! Use A storage if possible
             if (overwrite_a_) then                 
                 eigs => a(:,1)
+                #:if rt.startswith('r') 
                 eigi => a(:,2)
+                #:endif
             else
                 allocate(eigs(n)#{if rt.startswith('r')}#,eigi(n)#{endif}#)
-            end if               
-            
-            local_eigs = .true.
-            lde = n
-            
-        end if
+            end if  
+        endif
         
         #:if rt.startswith('c')
         allocate(rwork(n))
         #:endif
 
         if (lde<n) then
-            err0 = linalg_state(this, LINALG_VALUE_ERROR, 'Insufficient eigenvalue array size=',lde, &
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Insufficient eigenvalue array size=',lde, &
                                                           'should be >=',n)
             goto 2
         end if  

From 0551a66e64a40f26d3b1b6659e540fd90fb0179c Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 02:58:37 -0600
Subject: [PATCH 08/24] add examples

---
 example/linalg/CMakeLists.txt            |  3 +++
 example/linalg/example_schur_complex.f90 | 33 +++++++++++++++++++++++
 example/linalg/example_schur_eigvals.f90 | 34 ++++++++++++++++++++++++
 example/linalg/example_schur_real.f90    | 33 +++++++++++++++++++++++
 4 files changed, 103 insertions(+)
 create mode 100644 example/linalg/example_schur_complex.f90
 create mode 100644 example/linalg/example_schur_eigvals.f90
 create mode 100644 example/linalg/example_schur_real.f90

diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt
index fc3f823d9..01fe71936 100644
--- a/example/linalg/CMakeLists.txt
+++ b/example/linalg/CMakeLists.txt
@@ -24,6 +24,9 @@ ADD_EXAMPLE(eigvalsh)
 ADD_EXAMPLE(trace)
 ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)
+ADD_EXAMPLE(schur_real)
+ADD_EXAMPLE(schur_complex)
+ADD_EXAMPLE(schur_eigvals)
 ADD_EXAMPLE(blas_gemv)
 ADD_EXAMPLE(lapack_getrf)
 ADD_EXAMPLE(lstsq1)
diff --git a/example/linalg/example_schur_complex.f90 b/example/linalg/example_schur_complex.f90
new file mode 100644
index 000000000..741f40bb8
--- /dev/null
+++ b/example/linalg/example_schur_complex.f90
@@ -0,0 +1,33 @@
+! This example demonstrates the Schur decomposition for a complex-valued matrix.
+program example_schur_complex
+  use stdlib_linalg, only: schur
+  use stdlib_linalg_constants, only: dp
+  implicit none
+  complex(dp), allocatable :: A(:,:), T(:,:), Z(:,:)
+  integer :: n
+
+  ! Initialize a complex-valued square matrix
+  n = 3
+  allocate(A(n,n), T(n,n), Z(n,n))
+  A = reshape([ &
+       (1.0_dp, 2.0_dp), (3.0_dp, -1.0_dp), (4.0_dp, 1.0_dp), &
+       (0.0_dp, -1.0_dp), (2.0_dp, 0.0_dp), (1.0_dp, -2.0_dp), &
+       (2.0_dp, 3.0_dp), (1.0_dp, 1.0_dp), (0.0_dp, -1.0_dp) ], shape=[n,n])
+
+  ! Compute the Schur decomposition: A = Z T Z^H
+  call schur(A, T, Z)
+
+  ! Output results
+  print *, "Original Matrix A:"
+  print *, A
+  print *, "Schur Form Matrix T:"
+  print *, T
+  print *, "Unitary Matrix Z:"
+  print *, Z
+
+  ! Test factorization: Z*T*Z^H = A
+  print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, conjg(transpose(Z)))) - A))
+
+  deallocate(A, T, Z)
+end program example_schur_complex
+
diff --git a/example/linalg/example_schur_eigvals.f90 b/example/linalg/example_schur_eigvals.f90
new file mode 100644
index 000000000..829165b22
--- /dev/null
+++ b/example/linalg/example_schur_eigvals.f90
@@ -0,0 +1,34 @@
+! This example includes eigenvalue computation in addition to 
+! the Schur decomposition for a randomly generated matrix.
+program example_schur_eigenvalues
+  use stdlib_linalg, only: schur
+  use stdlib_linalg_constants, only: dp
+  implicit none
+  real(dp), allocatable :: A(:,:), T(:,:), Z(:,:)
+  complex(dp), allocatable :: eigenvalues(:)
+  integer :: n
+
+  ! Create a random real-valued square matrix
+  n = 5
+  allocate(A(n,n), T(n,n), Z(n,n), eigenvalues(n))
+  call random_number(A)
+
+  ! Compute the Schur decomposition and eigenvalues
+  call schur(A, T, Z, eigenvalues)
+
+  ! Output results
+  print *, "Random Matrix A:"
+  print *, A
+  print *, "Schur Form Matrix T:"
+  print *, T
+  print *, "Orthogonal Matrix Z:"
+  print *, Z
+  print *, "Eigenvalues:"
+  print *, eigenvalues
+
+  ! Test factorization: Z*T*Z^T = A
+  print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
+
+  deallocate(A, T, Z, eigenvalues)
+end program example_schur_eigenvalues
+
diff --git a/example/linalg/example_schur_real.f90 b/example/linalg/example_schur_real.f90
new file mode 100644
index 000000000..140ee5ce6
--- /dev/null
+++ b/example/linalg/example_schur_real.f90
@@ -0,0 +1,33 @@
+! This example computes the Schur decomposition of a real-valued square matrix.
+program example_schur_real
+  use stdlib_linalg, only: schur
+  use stdlib_linalg_constants, only: dp
+  implicit none
+  real(dp), allocatable :: A(:,:), T(:,:), Z(:,:)
+  integer :: n
+
+  ! Initialize a real-valued square matrix
+  n = 3
+  allocate(A(n,n), T(n,n), Z(n,n))
+  A = reshape([ &
+       0.0_dp, 2.0_dp, 2.0_dp, &
+       0.0_dp, 1.0_dp, 2.0_dp, &
+       1.0_dp, 0.0_dp, 1.0_dp], shape=[n,n])
+
+  ! Compute the Schur decomposition: A = Z T Z^T
+  call schur(A, T, Z)
+
+  ! Output results
+  print *, "Original Matrix A:"
+  print *, A
+  print *, "Schur Form Matrix T:"
+  print *, T
+  print *, "Orthogonal Matrix Z:"
+  print *, Z
+
+  ! Test factorization: Z*T*Z^T = A
+  print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
+
+  deallocate(A, T, Z)
+end program example_schur_real
+

From 14f771b97bf40c346c7596168f04243788abea4b Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 04:07:52 -0600
Subject: [PATCH 09/24] fix: only return eigenvalues if present

---
 src/stdlib_linalg_schur.fypp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 26ca20850..4695b782e 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -275,7 +275,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
       2 eigenvalue_output: if (local_eigs) then 
            #:if rt.startswith('r')
            ! Build complex eigenvalues
-           eigvals = cmplx(eigs,eigi,kind=${rk}$)
+           if (present(eigvals)) eigvals = cmplx(eigs,eigi,kind=${rk}$)
            #:endif
            if (.not.overwrite_a_) deallocate(eigs#{if rt.startswith('r')}#,eigi#{endif}#)
         endif eigenvalue_output

From 80a1df19214929bfcc17d7c1c1f9189b10a3a2ca Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 07:07:43 -0600
Subject: [PATCH 10/24] add tests

---
 test/linalg/CMakeLists.txt         |   4 +-
 test/linalg/test_linalg_schur.fypp | 193 +++++++++++++++++++++++++++++
 2 files changed, 196 insertions(+), 1 deletion(-)
 create mode 100644 test/linalg/test_linalg_schur.fypp

diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt
index ad41e5bb0..ecf645f74 100644
--- a/test/linalg/CMakeLists.txt
+++ b/test/linalg/CMakeLists.txt
@@ -10,6 +10,7 @@ set(
   "test_linalg_norm.fypp"
   "test_linalg_determinant.fypp"
   "test_linalg_qr.fypp"
+  "test_linalg_schur.fypp"
   "test_linalg_svd.fypp"
   "test_linalg_matrix_property_checks.fypp"
   "test_linalg_sparse.fypp"
@@ -26,6 +27,7 @@ ADDTEST(linalg_norm)
 ADDTEST(linalg_solve)
 ADDTEST(linalg_lstsq)
 ADDTEST(linalg_qr)
+ADDTEST(linalg_schur)
 ADDTEST(linalg_svd)
 ADDTEST(blas_lapack)
-ADDTEST(linalg_sparse)
\ No newline at end of file
+ADDTEST(linalg_sparse)
diff --git a/test/linalg/test_linalg_schur.fypp b/test/linalg/test_linalg_schur.fypp
new file mode 100644
index 000000000..0e3f63e5b
--- /dev/null
+++ b/test/linalg/test_linalg_schur.fypp
@@ -0,0 +1,193 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+! Test Schur decomposition
+module test_linalg_schur
+    use testdrive, only: error_type, check, new_unittest, unittest_type
+    use stdlib_linalg_constants
+    use stdlib_linalg_state, only: LINALG_VALUE_ERROR,linalg_state_type
+    use stdlib_linalg, only: schur,schur_space
+    use ieee_arithmetic, only: ieee_value,ieee_quiet_nan
+
+    implicit none (type,external)
+    
+    public :: test_schur_decomposition
+
+    contains
+
+    !> schur decomposition tests
+    subroutine test_schur_decomposition(tests)
+        !> Collection of tests
+        type(unittest_type), allocatable, intent(out) :: tests(:)
+        
+        allocate(tests(0))
+        
+        #:for rk,rt,ri in RC_KINDS_TYPES
+        tests = [tests,new_unittest("schur_api_${ri}$",test_schur_api_${ri}$), &
+                       new_unittest("schur_random_${ri}$",test_schur_random_${ri}$)]
+        #:endfor
+
+    end subroutine test_schur_decomposition
+
+    !> schur decomposition of a random matrix
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    subroutine test_schur_api_${ri}$(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        integer(ilp), parameter :: n   = 15_ilp
+        integer(ilp) :: lwork
+        real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$))
+        complex(${rk}$) :: eigs(n)
+        ${rt}$, dimension(n,n) :: a,t,z
+        ${rt}$, allocatable :: storage(:)
+        #:if 'complex' in rt
+        real(${rk}$) :: rea(n,n),ima(n,n)
+        #:endif
+        type(linalg_state_type) :: state
+        
+        #:if 'complex' in rt
+        call random_number(rea)
+        call random_number(ima)
+        a = cmplx(rea,ima,kind=${rk}$)
+        #:else
+        call random_number(a)
+        #:endif
+        
+        ! Test simple API
+        call schur(a,t,err=state)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return      
+        
+        ! Test output transformation matrix
+        call schur(a,t,z,err=state)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return      
+
+        ! Test output eigenvalues
+        call schur(a,t,eigvals=eigs,err=state)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return      
+        
+        ! Test storage query
+        call schur_space(a,lwork,err=state)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return      
+        
+        ! Test with user-defined storage
+        allocate(storage(lwork))        
+        call schur(a,t,eigvals=eigs,storage=storage,err=state)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return      
+        
+    end subroutine test_schur_api_${ri}$    
+    
+    subroutine test_schur_random_${ri}$(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        integer(ilp), parameter :: n = 3_ilp
+        real(${rk}$), parameter :: rtol = 1.0e-4_${rk}$
+        real(${rk}$), parameter :: eps  = sqrt(epsilon(0.0_${rk}$))
+        integer(ilp) :: lwork
+        ${rt}$, allocatable :: storage(:)
+        ${rt}$, dimension(n,n) :: a,t,z,aorig
+        #:if 'complex' in rt
+        real(${rk}$), dimension(n,n) :: a_re,a_im
+        #:endif        
+        type(linalg_state_type) :: state
+
+        #:if 'complex' in rt
+        call random_number(a_re)
+        call random_number(a_im)
+        a = cmplx(a_re,a_im,kind=${rk}$)
+        #:else
+        call random_number(a)
+        #:endif
+        aorig = a
+
+        ! 1) Run schur (standard)
+        call schur(a,t,z,err=state)
+        
+        ! Check return code
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return            
+        
+        ! Check solution
+        call check(error, all(schur_error(a,z,t)<=max(rtol*abs(a),eps)), &
+                          'converged solution (${rt}$)')
+        if (allocated(error)) return             
+
+        ! 2) Run schur (overwrite A)
+        call schur(a,t,z,overwrite_a=.true.,err=state)
+
+        ! Check return code
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return            
+        
+        ! Check solution
+        call check(error, all(schur_error(aorig,z,t)<=max(rtol*abs(aorig),eps)), &
+                          'converged solution (${rt}$ - overwrite A)')
+        if (allocated(error)) return            
+        
+        ! 3) Use working storage
+        a = aorig   
+        call schur_space(a,lwork,err=state)
+         
+        ! Check return code
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return
+        
+        allocate(storage(lwork))         
+        call schur(a,t,z,storage=storage,err=state)
+        
+        ! Check return code
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return            
+        
+        ! Check solution
+        call check(error, all(schur_error(a,z,t)<=max(rtol*abs(a),eps)), &
+                          'converged solution (${rt}$ - external storage)')
+        if (allocated(error)) return            
+
+    contains
+    
+        pure function schur_error(a,z,t) result(err)
+            ${rt}$, intent(in), dimension(:,:) :: a,z,t
+            real(${rk}$), dimension(size(a,1),size(a,2)) :: err
+            
+            #:if rt.startswith('real')
+            err = abs(matmul(matmul(z,t),transpose(z)) - a)
+            #:else
+            err = abs(matmul(matmul(z,t),conjg(transpose(z))) - a)
+            #:endif                
+        end function schur_error
+        
+    end subroutine test_schur_random_${ri}$
+
+    #:endfor
+
+end module test_linalg_schur
+
+program test_schur
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_linalg_schur, only : test_schur_decomposition
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+
+    testsuites = [ &
+        new_testsuite("linalg_schur", test_schur_decomposition) &
+        ]
+
+    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_schur

From 99ca1063c2f04f38eaba5c4e8aa74e95e831d5d2 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Wed, 20 Nov 2024 07:24:55 -0600
Subject: [PATCH 11/24] documentation

---
 doc/specs/stdlib_linalg.md | 75 ++++++++++++++++++++++++++++++++++++++
 src/stdlib_linalg.fypp     |  4 +-
 2 files changed, 76 insertions(+), 3 deletions(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 2f03f0fad..83b8da3f5 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -979,6 +979,81 @@ This subroutine computes the internal working space requirements for the QR fact
 {!example/linalg/example_qr_space.f90!}
 ```
 
+## `schur` - Compute the Schur decomposition of a matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is unitary (orthonormal) and \( T \) is upper-triangular (for complex matrices) or quasi-upper-triangular (for real matrices, with possible \( 2 \times 2 \) blocks on the diagonal). Matrix \( A \) has size `[n,n]`.
+
+The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \) is the Schur form, and matrix \( Z \) is the unitary transformation matrix such that \( A = Z T Z^H \). If requested, the eigenvalues of \( T \) can also be returned as a `complex` array of size `[n]`.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):schur(interface)]] `(a, t, z, [, eigvals] [, overwrite_a] [, storage] [, err])`
+
+### Arguments
+
+- `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(inout)` argument if `overwrite_a = .true.`; otherwise, it is an `intent(in)` argument.
+
+- `t`: Shall be a rank-2 array of the same kind as `a`, containing the Schur form \( T \) of the matrix. It is an `intent(out)` argument and should have a shape of `[n,n]`.
+
+- `z`: Shall be a rank-2 array of the same kind as `a`, containing the unitary matrix \( Z \). It is an `intent(out)` argument and is optional. If provided, it should have the shape `[n,n]`.
+
+- `eigvals` (optional): Shall be a rank-1 `complex` array containing the eigenvalues of \( A \) (the diagonal elements of \( T \)). The array must be of size `[n]`. If not provided, the eigenvalues are not returned.
+
+- `overwrite_a` (optional): Shall be a `logical` flag (default: `.false.`). If `.true.`, the input matrix `a` will be overwritten and destroyed upon return, avoiding internal data allocation. It is an `intent(in)` argument.
+
+- `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to [[stdlib_linalg(module):schur_space(interface)]]. It is an `intent(inout)` argument.
+
+- `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. If not provided, exceptions trigger an `error stop`.
+
+### Return value
+
+Returns the Schur decomposition matrices into the \( T \) and \( Z \) arguments. If `eigvals` is provided, it will also return the eigenvalues of the matrix \( A \).
+
+Raises `LINALG_VALUE_ERROR` if any of the matrices have invalid or unsuitable size for the decomposition.  
+Raises `LINALG_ERROR` on insufficient user storage space.  
+If the state argument `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_schur_eigvals.f90!}
+```
+
+---
+
+## `schur_space` - Compute internal working space requirements for the Schur decomposition
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the internal working space requirements for the Schur decomposition, [[stdlib_linalg(module):schur(interface)]].
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):schur_space(interface)]] `(a, lwork, [, err])`
+
+### Arguments
+
+- `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(in)` argument.
+
+- `lwork`: Shall be an `integer` scalar that returns the minimum array size required for the working storage in [[stdlib_linalg(module):schur(interface)]] to decompose `a`.
+
+- `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
+
+### Return value
+
+Returns the required working storage size `lwork` for the Schur decomposition. This value can be used to pre-allocate a workspace array in case multiple Schur decompositions of the same matrix size are needed. If pre-allocated working arrays are provided, no internal memory allocations will take place during the decomposition.
+
+
 ## `eig` - Eigenvalues and Eigenvectors of a Square Matrix
 
 ### Status
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index f896eff57..cbeba9cbf 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -619,10 +619,8 @@ module stdlib_linalg
     !!
     !! The output matrix \( T \) is upper-triangular for `complex` input matrices and quasi-upper-triangular
     !! for `real` input matrices (with possible \( 2 \times 2 \) blocks on the diagonal).
-    !! The user can optionally request sorting of eigenvalues based on conditions, or a custom sorting function.
     !!
-    !!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`). Sorting options 
-    !! are implemented using LAPACK's eigenvalue sorting mechanism.
+    !!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`). 
     !! 
     #:for rk,rt,ri in RC_KINDS_TYPES
       module subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, overwrite_a, storage, err)

From ec3c6c3b539bead964daa14f6f0b0401df9945a7 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 10:49:46 +0100
Subject: [PATCH 12/24] docs: `Z` is an optional argument

---
 doc/specs/stdlib_linalg.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 83b8da3f5..091c14e0e 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -993,7 +993,7 @@ The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \)
 
 ### Syntax
 
-`call ` [[stdlib_linalg(module):schur(interface)]] `(a, t, z, [, eigvals] [, overwrite_a] [, storage] [, err])`
+`call ` [[stdlib_linalg(module):schur(interface)]] `(a, t [, z,] [, eigvals] [, overwrite_a] [, storage] [, err])`
 
 ### Arguments
 

From 0ee05fc76f29137a667e3bf3821768c6eb72f5a9 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 10:51:56 +0100
Subject: [PATCH 13/24] make storage `intent(out)`

---
 doc/specs/stdlib_linalg.md   | 2 +-
 src/stdlib_linalg_schur.fypp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 091c14e0e..ef69fa5cc 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -1007,7 +1007,7 @@ The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \)
 
 - `overwrite_a` (optional): Shall be a `logical` flag (default: `.false.`). If `.true.`, the input matrix `a` will be overwritten and destroyed upon return, avoiding internal data allocation. It is an `intent(in)` argument.
 
-- `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to [[stdlib_linalg(module):schur_space(interface)]]. It is an `intent(inout)` argument.
+- `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to [[stdlib_linalg(module):schur_space(interface)]]. It is an `intent(out)` argument.
 
 - `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. If not provided, exceptions trigger an `error stop`.
 
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 4695b782e..440e0d28d 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -136,7 +136,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         !> [optional] Output eigenvalues that appear on the diagonal of T
         complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
-        ${rt}$, optional, intent(inout), target :: storage(:)
+        ${rt}$, optional, intent(out), target :: storage(:)
         !> [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

From 19e2edeecb66d467e2947706830269f7e670f884 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 10:52:09 +0100
Subject: [PATCH 14/24] Update doc/specs/stdlib_linalg.md

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
---
 doc/specs/stdlib_linalg.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 83b8da3f5..b6ad20fc3 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -1045,7 +1045,7 @@ This subroutine computes the internal working space requirements for the Schur d
 
 - `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(in)` argument.
 
-- `lwork`: Shall be an `integer` scalar that returns the minimum array size required for the working storage in [[stdlib_linalg(module):schur(interface)]] to decompose `a`.
+- `lwork`: Shall be an `integer` scalar that returns the minimum array size required for the working storage in [[stdlib_linalg(module):schur(interface)]] to decompose `a`. It is an `intent(out)` argument.
 
 - `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
 

From 21fa91f28caa8fc7a073d6fb944899897762b289 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 10:54:17 +0100
Subject: [PATCH 15/24] tidy up complex example

---
 example/linalg/example_schur_complex.f90 | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/example/linalg/example_schur_complex.f90 b/example/linalg/example_schur_complex.f90
index 741f40bb8..dc4f4e913 100644
--- a/example/linalg/example_schur_complex.f90
+++ b/example/linalg/example_schur_complex.f90
@@ -3,16 +3,14 @@ program example_schur_complex
   use stdlib_linalg, only: schur
   use stdlib_linalg_constants, only: dp
   implicit none
-  complex(dp), allocatable :: A(:,:), T(:,:), Z(:,:)
-  integer :: n
+  
+  integer, parameter :: n = 3
+  complex(dp), dimension(n,n) :: A, T, Z
 
   ! Initialize a complex-valued square matrix
-  n = 3
-  allocate(A(n,n), T(n,n), Z(n,n))
-  A = reshape([ &
-       (1.0_dp, 2.0_dp), (3.0_dp, -1.0_dp), (4.0_dp, 1.0_dp), &
-       (0.0_dp, -1.0_dp), (2.0_dp, 0.0_dp), (1.0_dp, -2.0_dp), &
-       (2.0_dp, 3.0_dp), (1.0_dp, 1.0_dp), (0.0_dp, -1.0_dp) ], shape=[n,n])
+  A = reshape([ (1, 2), (3,-1), (4, 1), &
+                (0,-1), (2, 0), (1,-2), &
+                (2, 3), (1, 1), (0,-1) ], shape=[n,n])
 
   ! Compute the Schur decomposition: A = Z T Z^H
   call schur(A, T, Z)
@@ -28,6 +26,5 @@ program example_schur_complex
   ! Test factorization: Z*T*Z^H = A
   print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, conjg(transpose(Z)))) - A))
 
-  deallocate(A, T, Z)
 end program example_schur_complex
 

From 7d8ca7d6f34e7c5a769df9eafefbaa0d7cc705fa Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 10:55:13 +0100
Subject: [PATCH 16/24] tidy up eigenvalues example

---
 example/linalg/example_schur_eigvals.f90 | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/example/linalg/example_schur_eigvals.f90 b/example/linalg/example_schur_eigvals.f90
index 829165b22..c62887439 100644
--- a/example/linalg/example_schur_eigvals.f90
+++ b/example/linalg/example_schur_eigvals.f90
@@ -4,13 +4,12 @@ program example_schur_eigenvalues
   use stdlib_linalg, only: schur
   use stdlib_linalg_constants, only: dp
   implicit none
-  real(dp), allocatable :: A(:,:), T(:,:), Z(:,:)
-  complex(dp), allocatable :: eigenvalues(:)
-  integer :: n
-
+  
+  integer, parameter :: n = 5  
+  real(dp), dimension(n,n) :: A, T, Z
+  complex(dp), dimension(n) :: eigenvalues
+  
   ! Create a random real-valued square matrix
-  n = 5
-  allocate(A(n,n), T(n,n), Z(n,n), eigenvalues(n))
   call random_number(A)
 
   ! Compute the Schur decomposition and eigenvalues
@@ -29,6 +28,5 @@ program example_schur_eigenvalues
   ! Test factorization: Z*T*Z^T = A
   print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
 
-  deallocate(A, T, Z, eigenvalues)
 end program example_schur_eigenvalues
 

From d594573f5dbfff6c0392e55ec6f780686f466ee8 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 10:56:07 +0100
Subject: [PATCH 17/24] tidy up schur real example

---
 example/linalg/example_schur_real.f90 | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/example/linalg/example_schur_real.f90 b/example/linalg/example_schur_real.f90
index 140ee5ce6..c46c730db 100644
--- a/example/linalg/example_schur_real.f90
+++ b/example/linalg/example_schur_real.f90
@@ -3,16 +3,13 @@ program example_schur_real
   use stdlib_linalg, only: schur
   use stdlib_linalg_constants, only: dp
   implicit none
-  real(dp), allocatable :: A(:,:), T(:,:), Z(:,:)
-  integer :: n
+  integer, parameter :: n = 3
+  real(dp), dimension(n,n) :: A, T, Z
 
   ! Initialize a real-valued square matrix
-  n = 3
-  allocate(A(n,n), T(n,n), Z(n,n))
-  A = reshape([ &
-       0.0_dp, 2.0_dp, 2.0_dp, &
-       0.0_dp, 1.0_dp, 2.0_dp, &
-       1.0_dp, 0.0_dp, 1.0_dp], shape=[n,n])
+  A = reshape([ 0, 2, 2, &
+                0, 1, 2, &
+                1, 0, 1], shape=[n,n])
 
   ! Compute the Schur decomposition: A = Z T Z^T
   call schur(A, T, Z)
@@ -28,6 +25,5 @@ program example_schur_real
   ! Test factorization: Z*T*Z^T = A
   print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
 
-  deallocate(A, T, Z)
 end program example_schur_real
 

From 8caa4dc877de9c0a5102952e53f951248bc06786 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 11:05:51 +0100
Subject: [PATCH 18/24] remove `goto`s

---
 src/stdlib_linalg_schur.fypp | 67 +++++++++++++++++++++---------------
 1 file changed, 39 insertions(+), 28 deletions(-)

diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 440e0d28d..c1b0e7eac 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -182,20 +182,6 @@ submodule (stdlib_linalg) stdlib_linalg_schur
             overwrite_a_ = .false._lk
         endif
         
-        !> SORTING: no sorting options are currently supported
-        sort  = gees_sort_eigs(.false.)
-        sdim  = 0_ilp     
-        
-        if (sort/=GEES_NOT) then 
-            
-           allocate(bwork(n),source=.false.)
-        
-        else
-            
-           bwork => bwork_dummy 
-            
-        end if        
-        
         !> Schur vectors
         jobvs = gees_vectors(present(z))        
         if (present(z)) then 
@@ -207,7 +193,8 @@ submodule (stdlib_linalg) stdlib_linalg_schur
             if (ldvs<n .or. nvs/=n) then
                 err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Schur vectors size=',[ldvs,nvs], &
                                                               'should be n=',n)
-                goto 1
+                call linalg_error_handling(err0, err)
+                return
             end if             
             
         else
@@ -226,10 +213,29 @@ submodule (stdlib_linalg) stdlib_linalg_schur
             
             ! Query optimal workspace            
             call get_schur_${ri}$_workspace(a,lwork,err0)
-            if (err0%error()) goto 1
-            allocate(work(lwork))
             
-        end if
+            if (err0%error()) then 
+                call linalg_error_handling(err0, err)
+                return
+            else
+                allocate(work(lwork))
+            end if
+            
+        end if        
+        
+        !> SORTING: no sorting options are currently supported
+        sort  = gees_sort_eigs(.false.)
+        sdim  = 0_ilp     
+        
+        if (sort/=GEES_NOT) then 
+            
+           allocate(bwork(n),source=.false.)
+        
+        else
+            
+           bwork => bwork_dummy 
+            
+        end if            
         
         !> User or self-allocated eigenvalue storage
         if (present(eigvals)) then             
@@ -262,17 +268,22 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         #:endif
 
         if (lde<n) then
-            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Insufficient eigenvalue array size=',lde, &
-                                                          'should be >=',n)
-            goto 2
+            
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, &
+                                           'Insufficient eigenvalue array size=',lde, &
+                                           'should be >=',n)
+        
+        else
+            
+            ! Compute Schur decomposition
+            call gees(jobvs,sort,eig_select,nt,t,mt,sdim,eigs,#{if rt.startswith('r')}#eigi,#{endif}# &
+                      vs,ldvs,work,lwork,#{if rt.startswith('c')}#rwork,#{endif}#bwork,info)
+            call handle_gees_info(info,m,n,m,err0)            
+            
+        
         end if  
 
-        ! Compute Schur decomposition
-        call gees(jobvs,sort,eig_select,nt,t,mt,sdim,eigs,#{if rt.startswith('r')}#eigi,#{endif}# &
-                  vs,ldvs,work,lwork,#{if rt.startswith('c')}#rwork,#{endif}#bwork,info)
-        call handle_gees_info(info,m,n,m,err0)
-
-      2 eigenvalue_output: if (local_eigs) then 
+        eigenvalue_output: if (local_eigs) then 
            #:if rt.startswith('r')
            ! Build complex eigenvalues
            if (present(eigvals)) eigvals = cmplx(eigs,eigi,kind=${rk}$)
@@ -280,7 +291,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
            if (.not.overwrite_a_) deallocate(eigs#{if rt.startswith('r')}#,eigi#{endif}#)
         endif eigenvalue_output
         if (.not.present(storage)) deallocate(work)
-      1 if (sort/=GEES_NOT) deallocate(bwork)
+        if (sort/=GEES_NOT) deallocate(bwork)
         call linalg_error_handling(err0,err)
         
         contains

From 170ed548138f5c4b6bd117ce260a2493e84fbfe6 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Thu, 12 Dec 2024 11:06:45 +0100
Subject: [PATCH 19/24] fix interface

---
 src/stdlib_linalg.fypp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index cbeba9cbf..6522e65b9 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -635,7 +635,7 @@ module stdlib_linalg
          !> [optional] Can A data be overwritten and destroyed?
          logical(lk), optional, intent(in) :: overwrite_a                 
          !> [optional] Provide pre-allocated workspace, size to be checked with schur_space         
-         ${rt}$, optional, intent(inout), target :: storage(:)
+         ${rt}$, optional, intent(out), target :: storage(:)
          !> [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_${ri}$_schur

From daf7a5e23b39ea6b521d63ba0029af752619a35d Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Mon, 16 Dec 2024 16:34:11 +0100
Subject: [PATCH 20/24] schur: `real` eigenvalues option

---
 src/stdlib_linalg.fypp       | 18 ++++++++++++++++
 src/stdlib_linalg_schur.fypp | 42 ++++++++++++++++++++++++++++++++++++
 2 files changed, 60 insertions(+)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 7ca844181..0dabff4ad 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -640,6 +640,24 @@ 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_${ri}$_schur
+      
+      ! Schur decomposition subroutine: real eigenvalue interface
+      module subroutine stdlib_linalg_real_eig_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
+         !> Input matrix a[m,m]
+         ${rt}$, intent(inout), target :: a(:,:)
+         !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
+         ${rt}$, intent(out), contiguous, target :: t(:,:)
+         !> Unitary/orthonormal transformation matrix Z
+         ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
+         !> Output real eigenvalues that appear on the diagonal of T
+         real(${rk}$), intent(out), contiguous, target :: eigvals(:)
+         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
+         ${rt}$, optional, intent(out), target :: storage(:)
+         !> [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}$_schur 
     #:endfor
   end interface schur
 
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index c1b0e7eac..a8180adc6 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -310,6 +310,48 @@ submodule (stdlib_linalg) stdlib_linalg_schur
 
     end subroutine stdlib_linalg_${ri}$_schur
 
+    ! Schur decomposition subroutine: real eigenvalue interface
+    module subroutine stdlib_linalg_real_eig_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
+        !> Input matrix a[m,m]
+        ${rt}$, intent(inout), target :: a(:,:)
+        !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
+        ${rt}$, intent(out), contiguous, target :: t(:,:)
+        !> Unitary/orthonormal transformation matrix Z
+        ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
+        !> Output eigenvalues that appear on the diagonal of T
+        real(${rk}$), intent(out), contiguous, target :: eigvals(:)
+        !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
+        ${rt}$, optional, intent(out), target :: storage(:)
+        !> [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 :: ceigvals(:)
+        real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$)
+        real(${rk}$), parameter :: atol = tiny(0.0_${rk}$)
+          
+        n = size(eigvals,dim=1,kind=ilp)
+        allocate(ceigvals(n))
+          
+        !> Compute Schur decomposition with complex eigenvalues
+        call stdlib_linalg_${ri}$_schur(a,t,z,ceigvals,overwrite_a,storage,err0)
+          
+        ! Check that no eigenvalues have meaningful imaginary part
+        if (err0%ok() .and. any(aimag(ceigvals)>atol+rtol*abs(abs(ceigvals)))) then 
+           err0 = linalg_state_type(this,LINALG_VALUE_ERROR, &
+                              'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(ceigvals)))
+        endif
+          
+        ! Return real components only
+        eigvals(:n) = real(ceigvals,kind=${rk}$)
+          
+        call linalg_error_handling(err0,err)        
+        
+    end subroutine stdlib_linalg_real_eig_${ri}$_schur
+
     #:endfor
 
 end submodule stdlib_linalg_schur

From 589d94d6b8b3cd135109a42fdf28c549ea58acc8 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Mon, 16 Dec 2024 16:34:20 +0100
Subject: [PATCH 21/24] test `real` eigenvalue option

---
 test/linalg/test_linalg_schur.fypp | 47 +++++++++++++++++++++++++++++-
 1 file changed, 46 insertions(+), 1 deletion(-)

diff --git a/test/linalg/test_linalg_schur.fypp b/test/linalg/test_linalg_schur.fypp
index 0e3f63e5b..139fdf360 100644
--- a/test/linalg/test_linalg_schur.fypp
+++ b/test/linalg/test_linalg_schur.fypp
@@ -23,7 +23,8 @@ module test_linalg_schur
         
         #:for rk,rt,ri in RC_KINDS_TYPES
         tests = [tests,new_unittest("schur_api_${ri}$",test_schur_api_${ri}$), &
-                       new_unittest("schur_random_${ri}$",test_schur_random_${ri}$)]
+                       new_unittest("schur_random_${ri}$",test_schur_random_${ri}$), &
+                       new_unittest("schur_symmetric_${ri}$",test_schur_symmetric_${ri}$)]
         #:endfor
 
     end subroutine test_schur_decomposition
@@ -162,6 +163,50 @@ module test_linalg_schur
         
     end subroutine test_schur_random_${ri}$
 
+    !> Test symmetric matrix (real eigenvalues)
+    subroutine test_schur_symmetric_${ri}$(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        integer(ilp), parameter :: n = 3_ilp
+        real(${rk}$), parameter :: rtol = 1.0e-4_${rk}$
+        real(${rk}$), parameter :: eps  = sqrt(epsilon(0.0_${rk}$))
+        real(${rk}$) :: reigs(n)
+        ${rt}$, dimension(n,n) :: a, t, z
+        type(linalg_state_type) :: state
+
+        ! Define a symmetric 3x3 matrix with real eigenvalues
+        a = reshape([ 3, 1, 0, &
+                      1, 3, 1, &
+                      0, 1, 3], shape=[n, n])
+
+        ! Return real eigenvalues (Should trigger an error if they have an imaginary part)
+        call schur(a, t, z, eigvals=reigs, err=state)
+        
+        ! Check return code
+        call check(error, state%ok(), state%print())
+        if (allocated(error)) return
+
+        ! Check solution
+        call check(error, all(schur_error(a, z, t) <= max(rtol * abs(a), eps)), &
+                          'converged solution (real symmetric, real eigs)')
+        if (allocated(error)) return
+
+    contains
+    
+        pure function schur_error(a,z,t) result(err)
+            ${rt}$, intent(in), dimension(:,:) :: a,z,t
+            real(${rk}$), dimension(size(a,1),size(a,2)) :: err
+            
+            #:if rt.startswith('real')
+            err = abs(matmul(matmul(z,t),transpose(z)) - a)
+            #:else
+            err = abs(matmul(matmul(z,t),conjg(transpose(z))) - a)
+            #:endif                
+        end function schur_error        
+
+    end subroutine test_schur_symmetric_${ri}$
+
+
     #:endfor
 
 end module test_linalg_schur

From 46a7ec2e188957d32463484f14060e9c172c1310 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Mon, 16 Dec 2024 16:38:12 +0100
Subject: [PATCH 22/24] document `real` eigenvalues case

---
 doc/specs/stdlib_linalg.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 94c9bf848..86cdbf8a6 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -1003,7 +1003,7 @@ The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \)
 
 - `z`: Shall be a rank-2 array of the same kind as `a`, containing the unitary matrix \( Z \). It is an `intent(out)` argument and is optional. If provided, it should have the shape `[n,n]`.
 
-- `eigvals` (optional): Shall be a rank-1 `complex` array containing the eigenvalues of \( A \) (the diagonal elements of \( T \)). The array must be of size `[n]`. If not provided, the eigenvalues are not returned.
+- `eigvals` (optional): Shall be a rank-1 `complex` or `real` array of the same kind as `a`, containing the eigenvalues of \( A \) (the diagonal elements of \( T \)), or their `real` component only. The array must be of size `[n]`. If not provided, the eigenvalues are not returned. It is an `intent(out)` argument.
 
 - `overwrite_a` (optional): Shall be a `logical` flag (default: `.false.`). If `.true.`, the input matrix `a` will be overwritten and destroyed upon return, avoiding internal data allocation. It is an `intent(in)` argument.
 
@@ -1016,6 +1016,7 @@ The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \)
 Returns the Schur decomposition matrices into the \( T \) and \( Z \) arguments. If `eigvals` is provided, it will also return the eigenvalues of the matrix \( A \).
 
 Raises `LINALG_VALUE_ERROR` if any of the matrices have invalid or unsuitable size for the decomposition.  
+Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts.
 Raises `LINALG_ERROR` on insufficient user storage space.  
 If the state argument `err` is not present, exceptions trigger an `error stop`.
 

From d5e4019d407834767e8359ace69bac2c3a553f88 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Fri, 27 Dec 2024 12:27:15 +0100
Subject: [PATCH 23/24] style fix

---
 src/stdlib_linalg_schur.fypp | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index a8180adc6..28095f0eb 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -176,11 +176,8 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         t = a 
         
         ! Can A be overwritten? By default, do not overwrite
-        if (present(overwrite_a)) then 
-            overwrite_a_ = overwrite_a .and. n>=2    
-        else
-            overwrite_a_ = .false._lk
-        endif
+        overwrite_a_ = .false._lk
+        if (present(overwrite_a)) overwrite_a_ = overwrite_a .and. n>=2    
         
         !> Schur vectors
         jobvs = gees_vectors(present(z))        

From f9a31cd1f2bb2566b9f08df8a22deafb394928cb Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Fri, 27 Dec 2024 12:28:30 +0100
Subject: [PATCH 24/24] `storage`: make `intent(inout)`

---
 doc/specs/stdlib_linalg.md   | 2 +-
 src/stdlib_linalg.fypp       | 4 ++--
 src/stdlib_linalg_schur.fypp | 4 ++--
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 35f3eeed1..ead216f2d 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -1035,7 +1035,7 @@ The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \)
 
 - `overwrite_a` (optional): Shall be a `logical` flag (default: `.false.`). If `.true.`, the input matrix `a` will be overwritten and destroyed upon return, avoiding internal data allocation. It is an `intent(in)` argument.
 
-- `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to [[stdlib_linalg(module):schur_space(interface)]]. It is an `intent(out)` argument.
+- `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to [[stdlib_linalg(module):schur_space(interface)]]. It is an `intent(inout)` argument.
 
 - `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. If not provided, exceptions trigger an `error stop`.
 
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 269c08fa8..2f4419a22 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -675,7 +675,7 @@ module stdlib_linalg
          !> [optional] Can A data be overwritten and destroyed?
          logical(lk), optional, intent(in) :: overwrite_a                 
          !> [optional] Provide pre-allocated workspace, size to be checked with schur_space         
-         ${rt}$, optional, intent(out), target :: storage(:)
+         ${rt}$, optional, intent(inout), target :: storage(:)
          !> [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_${ri}$_schur
@@ -691,7 +691,7 @@ module stdlib_linalg
          !> Output real eigenvalues that appear on the diagonal of T
          real(${rk}$), intent(out), contiguous, target :: eigvals(:)
          !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
-         ${rt}$, optional, intent(out), target :: storage(:)
+         ${rt}$, optional, intent(inout), target :: storage(:)
          !> [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
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 28095f0eb..3e477a319 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -136,7 +136,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         !> [optional] Output eigenvalues that appear on the diagonal of T
         complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
-        ${rt}$, optional, intent(out), target :: storage(:)
+        ${rt}$, optional, intent(inout), target :: storage(:)
         !> [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
@@ -318,7 +318,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         !> Output eigenvalues that appear on the diagonal of T
         real(${rk}$), intent(out), contiguous, target :: eigvals(:)
         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
-        ${rt}$, optional, intent(out), target :: storage(:)
+        ${rt}$, optional, intent(inout), target :: storage(:)
         !> [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