la_schur.fypp

#:include "common.fypp"
module la_schur
    use la_constants
    use la_lapack, only: gees
    use la_state_type, only: la_state, LINALG_ERROR, LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
    implicit none(type,external)
    private
    
    public :: schur
    public :: schur_space

    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' 
    
  ! 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/la_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 ALL_KINDS_TYPES
      module procedure la_${ri}$_schur
      module procedure la_real_eig_${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/la_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 ALL_KINDS_TYPES
      module procedure get_schur_${ri}$_workspace
    #:endfor
  end interface schur_space
    
    
    
    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
        type(la_state), intent(out) :: err

        ! Process GEES output
        select case (info)
        case (0_ilp)
            ! Success
        case (-1_ilp)
            ! Vector not wanted, but task is wrong
            err = la_state(this, LINALG_INTERNAL_ERROR,'Invalid Schur vector task request')
        case (-2_ilp)
            ! Vector not wanted, but task is wrong
            err = la_state(this, LINALG_INTERNAL_ERROR,'Invalid sorting task request')            
        case (-4_ilp,-6_ilp)
            ! Vector not wanted, but task is wrong
            err = la_state(this, LINALG_VALUE_ERROR,'Invalid/non-square input matrix size:',[m,n])  
        case (-11_ilp)
            err = la_state(this, LINALG_VALUE_ERROR,'Schur vector matrix has insufficient size',[ldvs,n])
        case (-13_ilp)
            err = la_state(this, LINALG_INTERNAL_ERROR,'Insufficient working storage size')
        case (1_ilp:)
            
            if (info==n+2) then 
                err = la_state(this, LINALG_ERROR, 'Ill-conditioned problem: could not sort eigenvalues')
            elseif (info==n+1) then 
                err = la_state(this, LINALG_ERROR, 'Some selected eigenvalues lost property due to sorting')                
            elseif (info==n) then 
                err = la_state(this, LINALG_ERROR, 'Convergence failure: no converged eigenvalues')
            else
                err = la_state(this, LINALG_ERROR, 'Convergence failure; converged range is',[info,n])                            
            end if
            
        case default
            
            err = la_state(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)

        end select
        
    end subroutine handle_gees_info
    
    
    #:for rk,rt,ri in ALL_KINDS_TYPES
    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(la_state), 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(la_state) :: 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 err0%handle(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
    subroutine la_${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(:,:)
        !> [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] 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(la_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(:)
        ${rt}$, target :: vs_dummy(1,1)
        ${rt}$, pointer :: vs(:,:),work(:),eigs(:)#{if rt.startswith('r')}#,eigi(:)#{endif}#
        character :: jobvs,sort
        type(la_state) :: err0
        
        ! Problem size
        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 (m/=n .or. m<=0 .or. n<=0) then
            err0 = la_state(this, LINALG_VALUE_ERROR, 'Matrix A must be square: size(a)=',[m,n])
            call err0%handle(err)
            return
        end if    
        if (mt/=nt .or. mt/=n .or. nt/=n) then
            err0 = la_state(this, LINALG_VALUE_ERROR, 'Matrix T must be square: size(T)=',[mt,nt], &
                                                          'should be',[m,n])
            call err0%handle(err)
            return
        end if                
        
        !> Copy data into the output array
        t = a 
        
        ! Can A be overwritten? By default, do not overwrite
        overwrite_a_ = .false._lk
        if (present(overwrite_a)) overwrite_a_ = overwrite_a .and. n>=2    
        
        !> 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 = la_state(this, LINALG_VALUE_ERROR, 'Schur vectors size=',[ldvs,nvs], &
                                                              'should be n=',n)
                call err0%handle(err)
                return
            end if             
            
        else
            vs => vs_dummy
            ldvs = size(vs, 1, kind=ilp)
            nvs  = size(vs, 2, kind=ilp)
        end if
        
        !> User or self-allocated storage
        if (present(storage)) then 
            
            work => storage
            lwork = size(work, 1, kind=ilp)
            
        else
            
            ! Query optimal workspace            
            call get_schur_${ri}$_workspace(a,lwork,err0)
            
            if (err0%error()) then 
                call err0%handle(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             
            lde = size(eigvals, 1, kind=ilp)            
            #:if rt.startswith('c')
            eigs => eigvals
            local_eigs = .false.
            #:else
            local_eigs = .true.
            #: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  
        endif
        
        #:if rt.startswith('c')
        allocate(rwork(n))
        #:endif

        if (lde<n) then
            
            err0 = la_state(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  

        eigenvalue_output: if (local_eigs) then 
           #:if rt.startswith('r')
           ! Build complex eigenvalues
           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
        if (.not.present(storage)) deallocate(work)
        if (sort/=GEES_NOT) deallocate(bwork)
        call err0%handle(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 

    end subroutine la_${ri}$_schur

    ! Schur decomposition subroutine: real eigenvalue interface
    module subroutine la_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(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(la_state), optional, intent(out) :: err
        
        type(la_state) :: 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 la_${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 = la_state(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 err0%handle(err)        
        
    end subroutine la_real_eig_${ri}$_schur

    #:endfor

end module la_schur