Add dggev. refs #4

gha3mi · Jul 22, 2023 · f464851 · f464851
1 parent 332d283
commit f464851
Showing 1 changed file with 55 additions and 0 deletions.
diff --git a/src/foreig.f90 b/src/foreig.f90
@@ -33,6 +33,8 @@ pure subroutine eig_rel(matrix, eig_vecr, eig_val, eig_vecl, method)
             call dsyev_rel(matrix, eig_vecr, eig_val)
           case ('geev')
             call dgeev_rel(matrix, eig_vecr, eig_val, eig_vecl)
+          case ('ggev')
+            call dggev_rel(matrix, eig_vecr, eig_val, eig_vecl)
          end select
 
       end if
@@ -125,4 +127,57 @@ 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)
+      real(rk), dimension(:,:), intent(in)                         :: matrix
+      real(rk), dimension(:),   allocatable, intent(out)           :: eig_val
+      real(rk), dimension(:,:), allocatable, intent(out)           :: eig_vecr
+      real(rk), dimension(:,:), allocatable, intent(out), optional :: eig_vecl
+      real(rk), dimension(size(matrix,1))                          :: alphar, alphai, beta
+      real(rk), dimension(size(matrix,1),size(matrix,1))           :: vl, vr
+      real(rk), dimension(:), allocatable                          :: work
+      integer                                                      :: m, lwork, info
+      real(rk)                                                     :: work1(1)
+
+      interface
+         pure subroutine dggev(fjobvl, fjobvr, fn, fa, flda, fb, ldb, &
+            falphar, falphai, fbeta, fvl, fldvl, fvr, fldvr, fwork, flwork, finfo)
+            import :: rk
+            character, intent(in)    :: fjobvl, fjobvr
+            integer,   intent(in)    :: fn, flda, ldb, fldvl, fldvr, flwork, finfo
+            real(rk),  intent(inout) :: fa(flda, *), fb(ldb, *)
+            real(rk),  intent(out)   :: falphar(fn), falphai(fn), fbeta(fn)
+            real(rk),  intent(inout) :: fvl(fldvl, *), fvr(fldvr, *)
+            real(rk),  intent(out)   :: fwork(flwork)
+         end subroutine
+      end interface
+
+      m = size(matrix, 1)
+      allocate(eig_val(m), eig_vecr(m, m))
+      eig_vecr = 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)
+      else
+         call dggev('N', 'V', m, eig_vecr, m, matrix, 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)
+      else
+         call dggev('N', 'V', m, eig_vecr, m, matrix, m, alphar, alphai, beta, vl, m, vr, m, work, lwork, info)
+      end if
+
+      eig_val = alphar / beta
+
+      deallocate(work)
+
+   end subroutine dggev_rel
+!===============================================================================
+
 end module foreig