#:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES #:set EIG_PROBLEM = ["standard", "generalized"] #:set EIG_FUNCTION = ["geev","ggev"] #:set EIG_PROBLEM_LIST = list(zip(EIG_PROBLEM, EIG_FUNCTION)) submodule (stdlib_linalg) stdlib_linalg_eigenvalues !! Compute eigenvalues and eigenvectors use stdlib_linalg_constants use stdlib_linalg_lapack, only: geev, ggev, heev, syev use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_positive_inf, ieee_quiet_nan implicit none(type,external) character(*), parameter :: this = 'eigenvalues' !> Utility function: Scale generalized eigenvalue interface scale_general_eig #:for rk,rt,ri in RC_KINDS_TYPES module procedure scale_general_eig_${ri}$ #:endfor end interface scale_general_eig contains !> Request for eigenvector calculation elemental character function eigenvectors_task(required) logical(lk), intent(in) :: required eigenvectors_task = merge('V','N',required) end function eigenvectors_task !> Request for symmetry side (default: lower) elemental character function symmetric_triangle_task(upper) logical(lk), optional, intent(in) :: upper symmetric_triangle_task = 'L' if (present(upper)) symmetric_triangle_task = merge('U','L',upper) end function symmetric_triangle_task !> Process GEEV output flags pure subroutine handle_geev_info(err,info,shapea) !> Error handler type(linalg_state_type), intent(inout) :: err !> GEEV return flag integer(ilp), intent(in) :: info !> Input matrix size integer(ilp), intent(in) :: shapea(2) select case (info) case (0) ! Success! err%state = LINALG_SUCCESS case (-1) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.') case (-2) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.') case (-5,-3) err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea) case (-9) err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.') case (-11) err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.') case (-13) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.') case (1:) err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.') case default err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.') end select end subroutine handle_geev_info !> Process GGEV output flags pure subroutine handle_ggev_info(err,info,shapea,shapeb) !> Error handler type(linalg_state_type), intent(inout) :: err !> GEEV return flag integer(ilp), intent(in) :: info !> Input matrix size integer(ilp), intent(in) :: shapea(2),shapeb(2) select case (info) case (0) ! Success! err%state = LINALG_SUCCESS case (-1) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.') case (-2) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.') case (-5,-3) err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea) case (-7) err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: b=',shapeb) case (-12) err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.') case (-14) err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.') case (-16) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.') case (1:) err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.') case default err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by ggev.') end select end subroutine handle_ggev_info !> Process SYEV/HEEV output flags elemental subroutine handle_heev_info(err,info,m,n) !> Error handler type(linalg_state_type), intent(inout) :: err !> SYEV/HEEV return flag integer(ilp), intent(in) :: info !> Input matrix size integer(ilp), intent(in) :: m,n select case (info) case (0) ! Success! err%state = LINALG_SUCCESS case (-1) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.') case (-2) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.') case (-5,-3) err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) case (-8) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.') case (1:) err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.') case default err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.') end select end subroutine handle_heev_info #:for rk,rt,ri in RC_KINDS_TYPES #:for ep,ei in EIG_PROBLEM_LIST module function stdlib_linalg_eigvals_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,err) result(lambda) !! Return an array of eigenvalues of matrix A. !> Input matrix A[m,n] ${rt}$, intent(in), dimension(:,:), target :: a #:if ei=='ggev' !> Generalized problem matrix B[n,n] ${rt}$, intent(inout), dimension(:,:), target :: b #:endif !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), intent(out) :: err !> Array of eigenvalues complex(${rk}$), allocatable :: lambda(:) !> Create ${rt}$, pointer, dimension(:,:) :: amat#{if ei=='ggev'}#, bmat #{endif}# integer(ilp) :: m,n,k !> Create an internal pointer so the intent of A won't affect the next call amat => a #{if ei=='ggev'}#bmat => b#{endif}# m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) k = min(m,n) !> Allocate return storage allocate(lambda(k)) !> Compute eigenvalues only call stdlib_linalg_eig_${ep}$_${ri}$(amat#{if ei=='ggev'}#,bmat#{endif}#,lambda,err=err) end function stdlib_linalg_eigvals_${ep}$_${ri}$ module function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#) result(lambda) !! Return an array of eigenvalues of matrix A. !> Input matrix A[m,n] ${rt}$, intent(in), dimension(:,:), target :: a #:if ei=='ggev' !> Generalized problem matrix B[n,n] ${rt}$, intent(inout), dimension(:,:), target :: b #:endif !> Array of eigenvalues complex(${rk}$), allocatable :: lambda(:) !> Create ${rt}$, pointer, dimension(:,:) :: amat#{if ei=='ggev'}#, bmat #{endif}# integer(ilp) :: m,n,k !> Create an internal pointer so the intent of A won't affect the next call amat => a #{if ei=='ggev'}#bmat => b#{endif}# m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) k = min(m,n) !> Allocate return storage allocate(lambda(k)) !> Compute eigenvalues only call stdlib_linalg_eig_${ep}$_${ri}$(amat#{if ei=='ggev'}#,bmat#{endif}#,lambda,overwrite_a=.false.) end function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$ module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left,& overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err) !! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, !! and optionally right or left eigenvectors. !> Input matrix A[m,n] ${rt}$, intent(inout), dimension(:,:), target :: a #:if ei=='ggev' !> Generalized problem matrix B[n,n] ${rt}$, intent(inout), dimension(:,:), target :: b #:endif !> Array of eigenvalues complex(${rk}$), intent(out) :: lambda(:) !> [optional] RIGHT eigenvectors of A (as columns) complex(${rk}$), optional, intent(out), target :: right(:,:) !> [optional] LEFT eigenvectors of A (as columns) complex(${rk}$), optional, intent(out), target :: left(:,:) !> [optional] Can A data be overwritten and destroyed? (default: no) logical(lk), optional, intent(in) :: overwrite_a #:if ei=='ggev' !> [optional] Can B data be overwritten and destroyed? (default: no) logical(lk), optional, intent(in) :: overwrite_b #:endif !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err !> Local variables type(linalg_state_type) :: err0 integer(ilp) :: m,n,lda,ldu,ldv,info,k,lwork,neig#{if ei=='ggev'}#,ldb,nb#{endif}# logical(lk) :: copy_a#{if ei=='ggev'}#,copy_b#{endif}# character :: task_u,task_v ${rt}$, target :: work_dummy(1),u_dummy(1,1),v_dummy(1,1) ${rt}$, allocatable :: work(:) ${rt}$, dimension(:,:), pointer :: amat,umat,vmat#{if ei=='ggev'}#,bmat#{endif}# #:if rt.startswith('complex') real(${rk}$), allocatable :: rwork(:) #:else ${rt}$, pointer :: lreal(:),limag(:) #:endif #:if ei=='ggev' ${rt}$, allocatable :: beta(:) #:endif !> Matrix size m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) k = min(m,n) neig = size(lambda,kind=ilp) lda = m if (k<=0 .or. m/=n) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& 'invalid or matrix size a=',[m,n],', must be nonempty square.') call linalg_error_handling(err0,err) return elseif (neig<k) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& 'eigenvalue array has insufficient size:',& ' lambda=',neig,', n=',n) call linalg_error_handling(err0,err) return endif #:if ep=='generalized' ldb = size(b,1,kind=ilp) nb = size(b,2,kind=ilp) if (ldb/=n .or. nb/=n) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& 'invalid or matrix size b=',[ldb,nb],', must be same as a=',[m,n]) call linalg_error_handling(err0,err) return end if #:endif ! Can A be overwritten? By default, do not overwrite copy_a = .true._lk if (present(overwrite_a)) copy_a = .not.overwrite_a ! Initialize a matrix temporary if (copy_a) then allocate(amat(m,n),source=a) else amat => a endif #:if ep=='generalized' ! Can B be overwritten? By default, do not overwrite copy_b = .true._lk if (present(overwrite_b)) copy_b = .not.overwrite_b ! Initialize a matrix temporary if (copy_b) then allocate(bmat,source=b) else bmat => b endif allocate(beta(n)) #:endif ! Decide if U, V eigenvectors task_u = eigenvectors_task(present(left)) task_v = eigenvectors_task(present(right)) if (present(right)) then #:if rt.startswith('complex') ! For a complex matrix, GEEV returns complex arrays. ! Point directly to output. vmat => right #:else ! For a real matrix, GEEV returns real arrays. ! Allocate temporary reals, will be converted to complex vectors at the end. allocate(vmat(n,n)) #:endif if (size(vmat,2,kind=ilp)<n) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& 'right eigenvector matrix has insufficient size: ',& shape(vmat),', with n=',n) endif else vmat => v_dummy endif if (present(left)) then #:if rt.startswith('complex') ! For a complex matrix, GEEV returns complex arrays. ! Point directly to output. umat => left #:else ! For a real matrix, GEEV returns real arrays. ! Allocate temporary reals, will be converted to complex vectors at the end. allocate(umat(n,n)) #:endif if (size(umat,2,kind=ilp)<n) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& 'left eigenvector matrix has insufficient size: ',& shape(umat),', with n=',n) endif else umat => u_dummy endif get_${ei}$: if (err0%ok()) then ldu = size(umat,1,kind=ilp) ldv = size(vmat,1,kind=ilp) ! Compute workspace size #:if rt.startswith('complex') allocate(rwork(2*n)) #:else allocate(lreal(n),limag(n)) #:endif lwork = -1_ilp call ${ei}$(task_u,task_v,n,amat,lda,& #:if ep=='generalized' bmat,ldb, & #:endif #{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# & #:if ep=='generalized' beta, & #:endif umat,ldu,vmat,ldv,& work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info) call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#) ! Compute eigenvalues if (info==0) then !> Prepare working storage lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) allocate(work(lwork)) !> Compute eigensystem call ${ei}$(task_u,task_v,n,amat,lda,& #:if ep=='generalized' bmat,ldb, & #:endif #{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# & #:if ep=='generalized' beta, & #:endif umat,ldu,vmat,ldv,& work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info) call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#) endif ! Finalize storage and process output flag #:if not rt.startswith('complex') lambda(:n) = cmplx(lreal(:n),limag(:n),kind=${rk}$) #:endif #:if ep=='generalized' ! Scale generalized eigenvalues lambda(:n) = scale_general_eig(lambda(:n),beta) #:endif #:if not rt.startswith('complex') ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, ! ${ei}$ returns reals as: ! u(j) = VL(:,j) + i*VL(:,j+1) and ! u(j+1) = VL(:,j) - i*VL(:,j+1). ! Convert these to complex numbers here. if (present(right)) call assign_real_eigenvectors_${rk}$(n,lambda,vmat,right) if (present(left)) call assign_real_eigenvectors_${rk}$(n,lambda,umat,left) #:endif endif get_${ei}$ if (copy_a) deallocate(amat) #:if ep=='generalized' if (copy_b) deallocate(bmat) #:endif #:if not rt.startswith('complex') if (present(right)) deallocate(vmat) if (present(left)) deallocate(umat) #:endif call linalg_error_handling(err0,err) end subroutine stdlib_linalg_eig_${ep}$_${ri}$ #:endfor module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] Should the upper/lower half of A be used? Default: lower logical(lk), optional, intent(in) :: upper_a !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), intent(out) :: err !> Array of eigenvalues real(${rk}$), allocatable :: lambda(:) ${rt}$, pointer :: amat(:,:) integer(ilp) :: m,n,k !> Create an internal pointer so the intent of A won't affect the next call amat => a m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) k = min(m,n) !> Allocate return storage allocate(lambda(k)) !> Compute eigenvalues only call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,err=err) end function stdlib_linalg_eigvalsh_${ri}$ module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda) !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] Should the upper/lower half of A be used? Default: use lower logical(lk), optional, intent(in) :: upper_a !> Array of eigenvalues real(${rk}$), allocatable :: lambda(:) ${rt}$, pointer :: amat(:,:) integer(ilp) :: m,n,k !> Create an internal pointer so the intent of A won't affect the next call amat => a m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) k = min(m,n) !> Allocate return storage allocate(lambda(k)) !> Compute eigenvalues only call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.) end function stdlib_linalg_eigvalsh_noerr_${ri}$ module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) !! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` !! of eigenvalues, and optionally right or left eigenvectors. !> Input matrix A[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Array of eigenvalues real(${rk}$), intent(out) :: lambda(:) !> The columns of vectors contain the orthonormal eigenvectors of A ${rt}$, optional, intent(out), target :: vectors(:,:) !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a !> [optional] Should the upper/lower half of A be used? Default: use lower logical(lk), optional, intent(in) :: upper_a !> [optional] 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,n,lda,info,k,lwork,neig logical(lk) :: copy_a character :: triangle,task ${rt}$, target :: work_dummy(1) ${rt}$, allocatable :: work(:) #:if rt.startswith('complex') real(${rk}$), allocatable :: rwork(:) #:endif ${rt}$, pointer :: amat(:,:) !> Matrix size m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) k = min(m,n) neig = size(lambda,kind=ilp) if (k<=0 .or. m/=n) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & ', must be non-empty square.') call linalg_error_handling(err0,err) return elseif (neig<k) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& ' lambda=',neig,' must be >= n=',n) call linalg_error_handling(err0,err) return endif ! Check if input A can be overwritten copy_a = .true._lk if (present(vectors)) then ! No need to copy A anyways copy_a = .false. elseif (present(overwrite_a)) then copy_a = .not.overwrite_a endif ! Should we use the upper or lower half of the matrix? triangle = symmetric_triangle_task(upper_a) ! Request for eigenvectors task = eigenvectors_task(present(vectors)) if (present(vectors)) then ! Check size if (any(shape(vectors,kind=ilp)<n)) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& 'eigenvector matrix has insufficient size: ',& shape(vectors),', with n=',n) call linalg_error_handling(err0,err) return endif ! The input matrix will be overwritten by the vectors. ! So, use this one as storage for syev/heev amat => vectors ! Copy data in amat(:n,:n) = a(:n,:n) elseif (copy_a) then ! Initialize a matrix temporary allocate(amat(m,n),source=a) else ! Overwrite A amat => a endif lda = size(amat,1,kind=ilp) ! Request workspace size lwork = -1_ilp #:if rt.startswith('complex') allocate(rwork(max(1,3*n-2))) call heev(task,triangle,n,amat,lda,lambda,work_dummy,lwork,rwork,info) #:else call syev(task,triangle,n,amat,lda,lambda,work_dummy,lwork,info) #:endif call handle_heev_info(err0,info,m,n) ! Compute eigenvalues if (info==0) then !> Prepare working storage lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) allocate(work(lwork)) !> Compute eigensystem #:if rt.startswith('complex') call heev(task,triangle,n,amat,lda,lambda,work,lwork,rwork,info) #:else call syev(task,triangle,n,amat,lda,lambda,work,lwork,info) #:endif call handle_heev_info(err0,info,m,n) endif ! Finalize storage and process output flag if (copy_a) deallocate(amat) call linalg_error_handling(err0,err) end subroutine stdlib_linalg_eigh_${ri}$ #:endfor #:for rk,rt,ri in REAL_KINDS_TYPES pure subroutine assign_real_eigenvectors_${rk}$(n,lambda,lmat,out_mat) !! GEEV for real matrices returns complex eigenvalues in real arrays, where two consecutive !! reals at [j,j+1] locations represent the real and imaginary parts of two complex conjugate !! eigenvalues. Convert them to complex here, following the GEEV logic. !> Problem size integer(ilp), intent(in) :: n !> Array of eigenvalues complex(${rk}$), intent(in) :: lambda(:) !> Real matrix as returned by geev ${rt}$, intent(in) :: lmat(:,:) !> Complex matrix as returned by eig complex(${rk}$), intent(out) :: out_mat(:,:) integer(ilp) :: i,j ! Copy matrix do concurrent(i=1:n,j=1:n) out_mat(i,j) = lmat(i,j) end do ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, ! geev returns them as reals as: ! u(j) = VL(:,j) + i*VL(:,j+1) and ! u(j+1) = VL(:,j) - i*VL(:,j+1). ! Convert these to complex numbers here. do j=1,n-1 if (lambda(j)==conjg(lambda(j+1))) then out_mat(:, j) = cmplx(lmat(:,j), lmat(:,j+1),kind=${rk}$) out_mat(:,j+1) = cmplx(lmat(:,j),-lmat(:,j+1),kind=${rk}$) endif end do end subroutine assign_real_eigenvectors_${rk}$ #:for ep,ei in EIG_PROBLEM_LIST module subroutine stdlib_linalg_real_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, & overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err) !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had !! non-trivial imaginary parts. !> Input matrix A[m,n] ${rt}$, intent(inout), target :: a(:,:) #:if ei=='ggev' !> Generalized problem matrix B[n,n] ${rt}$, intent(inout), target :: b(:,:) #:endif !> Array of real eigenvalues real(${rk}$), intent(out) :: lambda(:) !> The columns of RIGHT contain the right eigenvectors of A complex(${rk}$), optional, intent(out), target :: right(:,:) !> The columns of LEFT contain the left eigenvectors of A complex(${rk}$), optional, intent(out), target :: left(:,:) !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a #:if ei=='ggev' !> [optional] Can B data be overwritten and destroyed? (default: no) logical(lk), optional, intent(in) :: overwrite_b #:endif !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err type(linalg_state_type) :: err0 integer(ilp) :: n complex(${rk}$), allocatable :: clambda(:) real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$) real(${rk}$), parameter :: atol = tiny(0.0_${rk}$) n = size(lambda,dim=1,kind=ilp) allocate(clambda(n)) call stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,clambda,right,left, & overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err0) ! Check that no eigenvalues have meaningful imaginary part if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR, & 'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(clambda))) endif ! Return real components only lambda(:n) = real(clambda,kind=${rk}$) call linalg_error_handling(err0,err) end subroutine stdlib_linalg_real_eig_${ep}$_${ri}$ #:endfor #:endfor #:for rk,rt,ri in RC_KINDS_TYPES !> Utility function: Scale generalized eigenvalue elemental complex(${rk}$) function scale_general_eig_${ri}$(alpha,beta) result(lambda) !! A generalized eigenvalue for a pair of matrices (A,B) is a scalar lambda or a ratio !! alpha/beta = lambda, such that A - lambda*B is singular. It is usually represented as the !! pair (alpha,beta), there is a reasonable interpretation for beta=0, and even for both !! being zero. complex(${rk}$), intent(in) :: alpha ${rt}$, intent(in) :: beta real (${rk}$), parameter :: rzero = 0.0_${rk}$ complex(${rk}$), parameter :: czero = (0.0_${rk}$,0.0_${rk}$) if (beta==#{if rt.startswith('real')}#r#{else}#c#{endif}#zero) then if (alpha/=czero) then lambda = cmplx(ieee_value(1.0_${rk}$, ieee_positive_inf), & ieee_value(1.0_${rk}$, ieee_positive_inf), kind=${rk}$) else lambda = ieee_value(1.0_${rk}$, ieee_quiet_nan) end if else lambda = alpha/beta end if end function scale_general_eig_${ri}$ #:endfor end submodule stdlib_linalg_eigenvalues