Compiler Explorer

Source code

program test
  use stdlib_linalg, only: eig
  use stdlib_linalg_lapack, only: ggev
    use stdlib_linalg_constants
    use stdlib_linalg
    use stdlib_linalg_state
    use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_positive_inf, ieee_quiet_nan
  implicit none    
    character(*), parameter :: this="eig"

integer, parameter :: n=2
  double complex, dimension(n,n) :: A_Z,S_Z
  double precision, dimension(n,n) :: A_D,S_D
   double complex, dimension(n,n) :: vecs_r
   double complex,dimension(n) :: eigs
  integer :: i

! Set matrix
  A_Z = reshape( [ [1, 6], &
                 [9, 2] ], [n,n] )

S_Z = 0
  do i= 1,2
    S_Z(i,i) = i
  end do

A_D = real(A_Z)
  S_D = real(S_Z)

call stdlib_linalg_eig_generalized_z(A_Z,S_Z,eigs,right=vecs_r) !Fails
  !call eig(A_Z,S_Z,eigs,right=vecs_r) !Fails
  write(*,*) eigs
  write(*,*) vecs_r

contains

!> 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

!> Request for eigenvector calculation
     elemental character function eigenvectors_task(required)
        logical, intent(in) :: required
        eigenvectors_task = merge('V','N',required)
     end function eigenvectors_task

subroutine stdlib_linalg_eig_generalized_z(a,b,lambda,right,left,&
                                                       overwrite_a,overwrite_b,err)
     !! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, 
     !! and optionally right or left eigenvectors.        
         !> Input matrix A[m,n]
         complex(dp), intent(inout), dimension(:,:), target :: a 
         !> Generalized problem matrix B[n,n]
         complex(dp), intent(inout), dimension(:,:), target :: b         
         !> Array of eigenvalues
         complex(dp), intent(out) :: lambda(:)
         !> [optional] RIGHT eigenvectors of A (as columns)
         complex(dp), optional, intent(out), target :: right(:,:)
         !> [optional] LEFT eigenvectors of A (as columns)
         complex(dp), optional, intent(out), target :: left(:,:)
         !> [optional] Can A data be overwritten and destroyed? (default: no)
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] Can B data be overwritten and destroyed? (default: no)
         logical(lk), optional, intent(in) :: overwrite_b
         !> [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,ldb,nb
         logical(lk) :: copy_a,copy_b
         character :: task_u,task_v
         complex(dp), target :: work_dummy(1),u_dummy(1,1),v_dummy(1,1)
         complex(dp), allocatable :: work(:)
         complex(dp), dimension(:,:), pointer :: amat,umat,vmat,bmat
         real(dp), allocatable :: rwork(:)
         complex(dp), allocatable :: beta(:)
         
         !> 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
         
         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

! 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

! 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))

! Decide if U, V eigenvectors
         task_u = eigenvectors_task(present(left))
         task_v = eigenvectors_task(present(right))         
         
         if (present(right)) then 
                        
            ! For a complex matrix, GEEV returns complex arrays. 
            ! Point directly to output.
            vmat => right
            
            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
            
            ! For a complex matrix, GEEV returns complex arrays. 
            ! Point directly to output.
            umat => left

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_ggev: if (err0%ok()) then

ldu = size(umat,1,kind=ilp)
             ldv = size(vmat,1,kind=ilp)

! Compute workspace size
             allocate(rwork(5*n))

lwork = -1_ilp
            
             call ggev(task_u,task_v,n,amat,lda,&
                       bmat,ldb, &
                       lambda,  &
                       beta, &
                       umat,ldu,vmat,ldv,&
                       work_dummy,lwork,rwork,info)
             call handle_ggev_info(err0,info,shape(amat),shape(bmat))

! Compute eigenvalues
             if (info==0) then

!> Prepare working storage
                lwork = nint(real(work_dummy(1),kind=dp), kind=ilp)
                allocate(work(lwork))

!> Compute eigensystem
                call ggev(task_u,task_v,n,amat,lda,&
                          bmat,ldb, &
                          lambda,  &
                          beta, &
                          umat,ldu,vmat,ldv,&            
                          work,lwork,rwork,info)
                call handle_ggev_info(err0,info,shape(amat),shape(bmat))

endif
             
             ! Finalize storage and process output flag
             
             ! Scale generalized eigenvalues
             lambda(:n) = scale_general_eig(lambda(:n),beta)
             
         
         endif get_ggev
         
         if (copy_a) deallocate(amat)
         if (copy_b) deallocate(bmat)
         call linalg_error_handling(err0,err)

end subroutine stdlib_linalg_eig_generalized_z
     
     !> Utility function: Scale generalized eigenvalue
     elemental complex(dp) function scale_general_eig(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(dp), intent(in) :: alpha
         complex(dp),          intent(in) :: beta
         
         real   (dp), parameter :: rzero = 0.0_dp
         complex(dp), parameter :: czero = (0.0_dp,0.0_dp)
         
         if (beta==czero) then 
            if (alpha/=czero) then 
                lambda = cmplx(ieee_value(1.0_dp, ieee_positive_inf), &
                               ieee_value(1.0_dp, ieee_positive_inf), kind=dp)
            else
                lambda = ieee_value(1.0_dp, ieee_quiet_nan)
            end if            
         else
            lambda = alpha/beta
         end if
         
     end function scale_general_eig

end program test