Skip to content

Commit

Permalink
Minor changes. refs #4
Browse files Browse the repository at this point in the history
  • Loading branch information
gha3mi committed Jul 22, 2023
1 parent 7d65a67 commit 179d356
Showing 1 changed file with 14 additions and 8 deletions.
22 changes: 14 additions & 8 deletions src/foreig.f90
Original file line number Diff line number Diff line change
Expand Up @@ -119,17 +119,17 @@ pure subroutine dgeev(fjobvl, fjobvr, fn, fA, flda, fwr, fwi, fvl, fldvl, fvr, f
else
call dgeev('N', 'V', m, A, m, wr, wi, vl, m, vr, m, work1, -1, info)
end if

lwork = nint(work1(1))
allocate(work(lwork))

if (present(eig_vecl)) then
call dgeev('V', 'V', m, A, m, wr, wi, vl, m, vr, m, work, lwork, info)
eig_vecl = vl
else
call dgeev('N', 'V', m, A, m, wr, wi, vl, m, vr, m, work, lwork, info)
end if

eig_val = wr
eig_vecr = vr

Expand All @@ -138,7 +138,7 @@ pure subroutine dgeev(fjobvl, fjobvr, fn, fA, flda, fwr, fwi, fvl, fldvl, fvr, f
end subroutine dgeev_rel
!===============================================================================


!===============================================================================
!> author: Seyed Ali Ghasemi
pure subroutine dggev_rel(matrix, eig_vecr, eig_val, eig_vecl)
Expand All @@ -151,6 +151,7 @@ pure subroutine dggev_rel(matrix, eig_vecr, eig_val, eig_vecl)
real(rk), dimension(:), allocatable :: work
integer :: m, lwork, info
real(rk) :: work1(1)
real(rk), dimension(size(matrix,1),size(matrix,1)) :: A

interface
pure subroutine dggev(fjobvl, fjobvr, fn, fa, flda, fb, ldb, &
Expand All @@ -168,21 +169,26 @@ pure subroutine dggev(fjobvl, fjobvr, fn, fa, flda, fb, ldb, &
m = size(matrix, 1)
allocate(eig_val(m), eig_vecr(m, m))
eig_vecr = matrix
A = matrix

if (present(eig_vecl)) then
eig_vecl = matrix
call dggev('V', 'V', m, eig_vecl, m, eig_vecr, m, alphar, alphai, beta, vl, m, vr, m, work1, -1, info)
call dggev('V', 'V', m, eig_vecl, m,&
eig_vecr, m, alphar, alphai, beta, vl, m, vr, m, work1, -1, info)
else
call dggev('N', 'V', m, eig_vecr, m, matrix, m, alphar, alphai, beta, vl, m, vr, m, work1, -1, info)
call dggev('N', 'V', m, eig_vecr, m,&
A, m, alphar, alphai, beta, vl, m, vr, m, work1, -1, info)
end if

lwork = nint(work1(1))
allocate(work(lwork))

if (present(eig_vecl)) then
call dggev('V', 'V', m, eig_vecl, m, eig_vecr, m, alphar, alphai, beta, vl, m, vr, m, work, lwork, info)
call dggev('V', 'V', m, eig_vecl, m, &
eig_vecr, m, alphar, alphai, beta, vl, m, vr, m, work, lwork, info)
else
call dggev('N', 'V', m, eig_vecr, m, matrix, m, alphar, alphai, beta, vl, m, vr, m, work, lwork, info)
call dggev('N', 'V', m, eig_vecr, m,&
A, m, alphar, alphai, beta, vl, m, vr, m, work, lwork, info)
end if

eig_val = alphar / beta
Expand Down

0 comments on commit 179d356

Please sign in to comment.