rpaeig Subroutine

  subroutine rpaeig(ee,vl,vr,apb,amb,scr,tamm_dancoff)

    use precision, only: dp
    use eigen, only: diag_symm_packed

    implicit none

    real(kind=dp), intent(out), dimension(:) :: ee
    real(kind=dp), intent(out), dimension(:,:) :: vl, vr
    real(kind=dp), intent(in), dimension(:,:) :: apb
    real(kind=dp), intent(inout), dimension(:,:) :: amb
    real(kind=dp), intent(inout), dimension(:) :: scr
    logical, intent(in) :: tamm_dancoff

    integer :: ierr, j, nvec

    nvec = ubound(vl, 2)

    if (tamm_dancoff) then

  !   Diagonailze A: VR
      if (nvec==1) then
        ee(1) = vr(1,1)
        vr(1,1) = 1.0_dp
      else
        call dtrttp('u', nvec, apb, nvec, scr, ierr)
        call diag_symm_packed(1,nvec,nvec,nvec,scr,ee,vr,ierr)
      end if
      return
    end if

!   sqrt(A-B)
    if (nvec==1) then
      amb(1,1) = sqrt(amb(1,1))
    else
      call dtrttp('u', nvec, amb, nvec, scr, ierr)
      call diag_symm_packed(1,nvec,nvec,nvec,scr,ee,vr,ierr)
      ee(1:nvec) = sign(sqrt(abs(ee(1:nvec))),ee(1:nvec))
      do j = 1, nvec
        vl(:,j) = vr(:,j)*ee(j)
      end do
      call dgemm('n','t',nvec,nvec,nvec, &
                 1.0_dp,vr,nvec,vl,nvec, &
                 0.0_dp,amb,nvec)
    end if

!   Form sqrt(A-B)*(A+B)*sqrt(A-B)
!   (A+B)*sqrt(A-B) : VL
    call dgemm('n','n',nvec,nvec,nvec, &
               1.0_dp,apb,nvec,amb,nvec, &
               0.0_dp,vl,nvec)
!   sqrt(A-B)*(A+B)*sqrt(A-B) : VR
    call dgemm('n','n',nvec,nvec,nvec, &
               1.0_dp,amb,nvec,vl,nvec, &
               0.0_dp,vr,nvec)

!   Diagonailze sqrt(A-B)*(A+B)*sqrt(A-B) : VR
    if (nvec==1) then
      ee(1) = vr(1,1)
      vl(1,1) = 1.0_dp
    else
      call dtrttp('u', nvec, vr, nvec, scr, ierr)
      call diag_symm_packed(1,nvec,nvec,nvec,scr,ee,vl,ierr)
    end if

!   Current vector VL is sqrt(1/(A-B)))|X+Y>.

!   VL into right eigenvector  VR = |X+Y>
    call dgemm('n','n',nvec,nvec,nvec, &
               1.0_dp,amb,nvec,vl,nvec, &
               0.0_dp,vr,nvec)

!   Left eigenvector   VL = = 1/E (A+B)|X+Y>
    call dgemm('n','n',nvec,nvec,nvec, &
               1.0_dp,apb,nvec,vr,nvec, &
               0.0_dp,vl,nvec)

    do j = 1, nvec
      vl(1:nvec,j) = vl(1:nvec,j)/sign(sqrt(abs(ee(j))),ee(j))
    end do

  end subroutine rpaeig