qp2/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f

subroutine lapack_diag_non_sym(n, A, WR, WI, VL, VR)

  BEGIN_DOC
  ! You enter with a general non hermitian matrix A(n,n) 
  !
  ! You get out with the real WR and imaginary part WI of the eigenvalues 
  !
  ! Eigvalue(n) = WR(n) + i * WI(n)
  !
  ! And the left VL and right VR eigenvectors 
  !
  ! VL(i,j) = <i|Psi_left(j)>  :: projection on the basis element |i> on the jth left  eigenvector 
  !
  ! VR(i,j) = <i|Psi_right(j)> :: projection on the basis element |i> on the jth right eigenvector 
  !
  ! The real part of the matrix A can be written as A = VR D VL^T
  !
  END_DOC

  implicit none

  integer,          intent(in)  :: n
  double precision, intent(in)  :: A(n,n)
  double precision, intent(out) :: WR(n), WI(n), VL(n,n), VR(n,n)

  integer                       :: lda, ldvl, ldvr, LWORK, INFO
  double precision, allocatable :: Atmp(:,:), WORK(:)

  lda  = n
  ldvl = n
  ldvr = n

  allocate( Atmp(n,n) )
  Atmp(1:n,1:n) = A(1:n,1:n)

  allocate(WORK(1))
  LWORK = -1 ! to ask for the optimal size of WORK
  call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
  if(INFO.gt.0)then
    print*,'dgeev failed !!',INFO
    stop
  endif
  LWORK = max(int(WORK(1)), 1) ! this is the optimal size of WORK 
  deallocate(WORK)

  allocate(WORK(LWORK))

  ! Actual diagonalization 
  call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
  if(INFO.ne.0) then
    print*,'dgeev failed !!', INFO
    stop
  endif

  deallocate(Atmp, WORK)

end subroutine lapack_diag_non_sym


subroutine non_sym_diag_inv_right(n,A,leigvec,reigvec,n_real_eigv,eigval)
 implicit none
 BEGIN_DOC
! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors 
!
! of a non hermitian matrix A(n,n)
!
! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n" 
 END_DOC
 integer, intent(in) :: n
 double precision, intent(in) :: A(n,n)
 double precision, intent(out) :: reigvec(n,n),leigvec(n,n),eigval(n)
 double precision, allocatable :: Aw(:,:)
 integer, intent(out) :: n_real_eigv
 print*,'Computing the left/right eigenvectors ...'
 character*1 :: JOBVL,JOBVR
 JOBVL = "V" ! computes the left  eigenvectors 
 JOBVR = "V" ! computes the right eigenvectors 
 double precision, allocatable :: WR(:),WI(:),Vl(:,:),VR(:,:),S(:,:),inv_reigvec(:,:)
 integer :: i,j
 integer :: n_good
 integer, allocatable :: list_good(:), iorder(:)
 double precision :: thr
 thr = 1.d-10
 ! Eigvalue(n) = WR(n) + i * WI(n)
 allocate(WR(n),WI(n),VL(n,n),VR(n,n),Aw(n,n))
 Aw = A
 do i = 1, n
  do j = i+1, n
   if(dabs(Aw(j,j)-Aw(i,i)).lt.thr)then
     Aw(j,j)+= thr
     Aw(i,i)-= thr
!    if(Aw(j,i) * A(i,j) .lt.0d0  )then
!     if(dabs(Aw(j,i) * A(i,j)).lt.thr**(1.5d0))then
!      print*,Aw(j,j),Aw(i,i)
!      print*,Aw(j,i) , A(i,j)
      Aw(j,i) = 0.d0
      Aw(i,j) = Aw(j,i)
!     endif
!    endif
   endif
  enddo
 enddo
 call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
 ! You track the real eigenvalues 
 n_good = 0
! do i = 1, n
!  write(*,'(100(F16.12,X))')A(:,i)
! enddo
 do i = 1, n
  print*,'Im part of lambda = ',dabs(WI(i))
  if(dabs(WI(i)).lt.thr)then
   n_good += 1
  else
   print*,'Found an imaginary component to eigenvalue'
   print*,'Re(i) + Im(i)',WR(i),WI(i)
   write(*,'(100(F10.5,X))')VR(:,i)
   write(*,'(100(F10.5,X))')VR(:,i+1)
   write(*,'(100(F10.5,X))')VL(:,i)
   write(*,'(100(F10.5,X))')VL(:,i+1)
  endif
 enddo
 allocate(list_good(n_good),iorder(n_good))
 n_good = 0
 do i = 1, n
  if(dabs(WI(i)).lt.thr)then
   n_good += 1
   list_good(n_good) = i
   eigval(n_good) = WR(i)
  endif
 enddo
 n_real_eigv = n_good 
 do i = 1, n_good
  iorder(i) = i
 enddo
 ! You sort the real eigenvalues 
 call dsort(eigval,iorder,n_good)
 do i = 1, n_real_eigv
  do j = 1, n
   reigvec(j,i) = VR(j,list_good(iorder(i)))
   leigvec(j,i) = VL(j,list_good(iorder(i)))
  enddo
 enddo
 allocate(inv_reigvec(n_real_eigv,n_real_eigv))
! call get_pseudo_inverse(reigvec,n_real_eigv,n_real_eigv,n_real_eigv,inv_reigvec,n_real_eigv,thr)
! do i = 1, n_real_eigv
!  do j = 1, n
!   leigvec(j,i) = inv_reigvec(i,j)
!  enddo
! enddo
 allocate( S(n_real_eigv,n_real_eigv) )

  ! S = VL x VR
  call dgemm( 'T', 'N', n_real_eigv, n_real_eigv, n_real_eigv, 1.d0                              &
            , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
            , 0.d0, S, size(S, 1) )
   do i = 1,n_real_eigv
    write(*,'(100(F10.5,X))')S(:,i)
   enddo
! call lapack_diag_non_sym(n,S,WR,WI,VL,VR)
! print*,'Eigenvalues of S'
! do i = 1, n
!  print*,WR(i),dabs(WI(i))
! enddo
  call dgemm( 'T', 'N', n_real_eigv, n_real_eigv, n_real_eigv, 1.d0                              &
            , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
            , 0.d0, S, size(S, 1) )
! call get_inv_half_svd(S, n_real_eigv, inv_reigvec)

  double precision :: accu_d,accu_nd
  accu_nd = 0.d0
  accu_d = 0.d0
  do i = 1, n_real_eigv
    do j = 1, n_real_eigv
      if(i==j) then
       accu_d += S(j,i) * S(j,i)
      else
       accu_nd = accu_nd + S(j,i) * S(j,i)
      endif
    enddo
  enddo
  accu_nd = dsqrt(accu_nd)

  print*,'accu_nd = ',accu_nd
  if( accu_nd .lt. 1d-10 ) then
    ! L x R is already bi-orthogonal
    !print *, ' L & T bi-orthogonality: ok'
    return
  else
   print*,'PB with bi-orthonormality!!'
   stop
  endif
end

subroutine lapack_diag_non_sym_new(n, A, WR, WI, VL, VR)

  BEGIN_DOC
  !
  ! You enter with a general non hermitian matrix A(n,n) 
  !
  ! You get out with the real WR and imaginary part WI of the eigenvalues 
  !
  ! Eigvalue(n) = WR(n) + i * WI(n)
  !
  ! And the left VL and right VR eigenvectors 
  !
  ! VL(i,j) = <i|Psi_left(j)>  :: projection on the basis element |i> on the jth left  eigenvector 
  !
  ! VR(i,j) = <i|Psi_right(j)> :: projection on the basis element |i> on the jth right eigenvector 
  !
  END_DOC

  implicit none
  integer,          intent(in)  :: n
  double precision, intent(in)  :: A(n,n)
  double precision, intent(out) :: WR(n), WI(n), VL(n,n), VR(n,n)

  character*1                   :: JOBVL,JOBVR,BALANC,SENSE
  integer                       :: ILO, IHI
  integer                       :: lda, ldvl, ldvr, LWORK, INFO
  double precision              :: ABNRM
  integer,          allocatable :: IWORK(:)
  double precision, allocatable :: WORK(:), SCALE_array(:), RCONDE(:), RCONDV(:)
  double precision, allocatable :: Atmp(:,:)

  allocate( Atmp(n,n) )
  Atmp(1:n,1:n) = A(1:n,1:n)

  JOBVL  = "V" ! computes the left  eigenvectors 
  JOBVR  = "V" ! computes the right eigenvectors 
  BALANC = "B" ! Diagonal scaling and Permutation for optimization
  SENSE  = "B"
  lda  = n
  ldvl = n
  ldvr = n
  allocate(WORK(1),SCALE_array(n),RCONDE(n),RCONDV(n),IWORK(2*n-2))
  LWORK = -1 ! to ask for the optimal size of WORK
  call dgeevx(BALANC,JOBVL,JOBVR,SENSE,&  ! CHARACTERS 
              n,Atmp,lda,              &  ! MATRIX TO DIAGONALIZE
              WR,WI,                   &  ! REAL AND IMAGINARY PART OF EIGENVALUES 
              VL,ldvl,VR,ldvr,         &  ! LEFT AND RIGHT EIGENVECTORS 
              ILO,IHI,SCALE_array,ABNRM,RCONDE,RCONDV, & ! OUTPUTS OF OPTIMIZATION
              WORK,LWORK,IWORK,INFO)

  !if(INFO.gt.0)then
  ! print*,'dgeev failed !!',INFO
  if( INFO.ne.0 ) then
    print *, 'dgeevx failed !!', INFO
    stop
  endif

  LWORK = max(int(work(1)), 1) ! this is the optimal size of WORK 
  deallocate(WORK)
  allocate(WORK(LWORK))
  ! Actual dnon_hrmt_real_diag_newiagonalization 
  call dgeevx(BALANC,JOBVL,JOBVR,SENSE,&  ! CHARACTERS 
              n,Atmp,lda,              &  ! MATRIX TO DIAGONALIZE
              WR,WI,                   &  ! REAL AND IMAGINARY PART OF EIGENVALUES 
              VL,ldvl,VR,ldvr,         &  ! LEFT AND RIGHT EIGENVECTORS 
              ILO,IHI,SCALE_array,ABNRM,RCONDE,RCONDV, & ! OUTPUTS OF OPTIMIZATION
              WORK,LWORK,IWORK,INFO)

  !if(INFO.ne.0)then
  ! print*,'dgeev failed !!',INFO
  if( INFO.ne.0 ) then
    print *, 'dgeevx failed !!', INFO
    stop
  endif

  deallocate( Atmp )
  deallocate( WORK, SCALE_array, RCONDE, RCONDV, IWORK )

end subroutine lapack_diag_non_sym_new

! ---

subroutine lapack_diag_non_sym_right(n, A, WR, WI, VR)

  implicit none

  integer,          intent(in)  :: n
  double precision, intent(in)  :: A(n,n)
  double precision, intent(out) :: WR(n), WI(n), VR(n,n)

  integer                       :: i, lda, ldvl, ldvr, LWORK, INFO
  double precision, allocatable :: Atmp(:,:), WORK(:), VL(:,:)

  lda  = n
  ldvl = 1
  ldvr = n

  allocate( Atmp(n,n), VL(1,1) )
  Atmp(1:n,1:n) = A(1:n,1:n)

  allocate(WORK(1))
  LWORK = -1
  call dgeev('N', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
  if(INFO.gt.0)then
    print*,'dgeev failed !!',INFO
    stop
  endif

  LWORK = max(int(WORK(1)), 1) ! this is the optimal size of WORK 
  deallocate(WORK)

  allocate(WORK(LWORK))

  ! Actual diagonalization 
  call dgeev('N', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
  if(INFO.ne.0) then
    print*,'dgeev failed !!', INFO
    stop
  endif

  deallocate(Atmp, WORK, VL)

! print *, ' JOBL = F'
! print *, ' eigenvalues'
! do i = 1, n
!   write(*, '(1000(F16.10,X))') WR(i), WI(i)
! enddo
! print *, ' right eigenvect' 
! do i = 1, n
!   write(*, '(1000(F16.10,X))') VR(:,i)
! enddo

end subroutine lapack_diag_non_sym_right

! ---

subroutine non_hrmt_real_diag(n, A, leigvec, reigvec, n_real_eigv, eigval)

  BEGIN_DOC
  !
  ! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors 
  !
  ! of a non hermitian matrix A(n,n)
  !
  ! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n" 
  !
  END_DOC

  implicit none
  integer,          intent(in)  :: n
  double precision, intent(in)  :: A(n,n)
  integer,          intent(out) :: n_real_eigv
  double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)

  integer                       :: i, j, n_good
  double precision              :: thr, threshold, accu_d, accu_nd
  integer,          allocatable :: list_good(:), iorder(:)
  double precision, allocatable :: Aw(:,:)
  double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:), S(:,:), S_inv_half_tmp(:,:)

  print*, ' Computing the left/right eigenvectors with lapack ...'

  ! Eigvalue(n) = WR(n) + i * WI(n)
  allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n))
  Aw = A
  !print *, ' matrix to diagonalize', Aw
  call lapack_diag_non_sym(n, Aw, WR, WI, VL, VR)

  ! ---
  ! You track the real eigenvalues 

  thr = 1d-15

  n_good = 0
  do i = 1, n
    if(dabs(WI(i)).lt.thr) then
      n_good += 1
    else
      print*, ' Found an imaginary component to eigenvalue'
      print*, ' Re(i) + Im(i)', WR(i), WI(i)
    endif
  enddo

  allocate(list_good(n_good), iorder(n_good))
  n_good = 0
  do i = 1, n
    if(dabs(WI(i)).lt.thr) then
      n_good += 1
      list_good(n_good) = i
      eigval(n_good) = WR(i)
    endif
  enddo
  n_real_eigv = n_good
  do i = 1, n_good
   iorder(i) = i
  enddo

  ! You sort the real eigenvalues 
  call dsort(eigval, iorder, n_good)
  do i = 1, n_real_eigv
    do j = 1, n
      reigvec(j,i) = VR(j,list_good(iorder(i)))
      leigvec(j,i) = Vl(j,list_good(iorder(i)))
    enddo
  enddo

! print *, ' ordered eigenvalues'
! print *, ' right eigenvect' 
! do i = 1, n
!   print *, i, eigval(i)
!   write(*, '(1000(F16.10,X))') reigvec(:,i)
! enddo

  ! ---

  allocate( S(n_real_eigv,n_real_eigv), S_inv_half_tmp(n_real_eigv,n_real_eigv) )

  ! S = VL x VR
  call dgemm( 'T', 'N', n_real_eigv, n_real_eigv, n_real_eigv, 1.d0 &
            , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1)  &
            , 0.d0, S, size(S, 1) )

  accu_nd = 0.d0
  accu_d  = 0.d0
  do i = 1, n_real_eigv
    do j = 1, n_real_eigv
      if(i==j) then
        accu_d += S(j,i)
      else
        accu_nd = accu_nd + S(j,i) * S(j,i)
      endif
    enddo
  enddo
  accu_nd = dsqrt(accu_nd)

  threshold = 1.d-15
  if( (accu_nd .gt. threshold) .or. (dabs(accu_d-dble(n_real_eigv)) .gt. threshold) ) then

    print*, ' sum of off-diag S elements = ', accu_nd
    print*, ' Should be zero '
    print*, ' sum of     diag S elements = ', accu_d
    print*, ' Should be ',n
    print*, ' Not bi-orthonormal !!'
    print*, ' Notice that if you are interested in ground state it is not a problem :)'
  endif

end subroutine non_hrmt_real_diag

! ---

subroutine lapack_diag_general_non_sym(n, A, B, WR, WI, VL, VR)

  BEGIN_DOC
  ! You enter with a general non hermitian matrix A(n,n) and another B(n,n)
  !
  ! You get out with the real WR and imaginary part WI of the eigenvalues 
  !
  ! Eigvalue(n) = (WR(n) + i * WI(n))
  !
  ! And the left VL and right VR eigenvectors 
  !
  ! VL(i,j) = <i|Psi_left(j)>  :: projection on the basis element |i> on the jth left  eigenvector 
  !
  ! VR(i,j) = <i|Psi_right(j)> :: projection on the basis element |i> on the jth right eigenvector 
  END_DOC

  implicit none
  integer,          intent(in)  :: n
  double precision, intent(in)  :: A(n,n), B(n,n)
  double precision, intent(out) :: WR(n), WI(n), VL(n,n), VR(n,n)

  integer                       :: lda, ldvl, ldvr, LWORK, INFO
  integer                       :: n_good
  double precision, allocatable :: WORK(:)
  double precision, allocatable :: Atmp(:,:)

  lda  = n
  ldvl = n
  ldvr = n

  allocate( Atmp(n,n) )
  Atmp(1:n,1:n) = A(1:n,1:n)

  allocate(WORK(1))
  LWORK = -1 
  call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
  if(INFO.gt.0) then
    print*,'dgeev failed !!',INFO
    stop
  endif

  LWORK = max(int(WORK(1)), 1) 
  deallocate(WORK)

  allocate(WORK(LWORK))

  call dgeev('V', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
  if(INFO.ne.0) then
    print*,'dgeev failed !!', INFO
    stop
  endif

  deallocate( WORK, Atmp )

end subroutine lapack_diag_general_non_sym

! ---

subroutine non_hrmt_general_real_diag(n, A, B, reigvec, leigvec, n_real_eigv, eigval)

  BEGIN_DOC
  ! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors 
  !
  ! of a non hermitian matrix A(n,n) and B(n,n) 
  !
  ! A reigvec = eigval * B * reigvec
  !
  ! (A)^\dagger leigvec = eigval * B * leigvec
  !
  ! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n" 
  END_DOC

  implicit none
  integer,          intent(in)  :: n
  double precision, intent(in)  :: A(n,n), B(n,n)
  integer,          intent(out) :: n_real_eigv
  double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)

  integer                       :: i, j
  integer                       :: n_good
  integer, allocatable          :: list_good(:), iorder(:)
  double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:)
  double precision, allocatable :: Aw(:,:), Bw(:,:)

  print*,'Computing the left/right eigenvectors ...'

  allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n), Bw(n,n))
  Aw = A
  Bw = B

  call lapack_diag_general_non_sym(n, A, B, WR, WI, VL, VR)

  ! You track the real eigenvalues 
  n_good = 0
  do i = 1, n
    if(dabs(WI(i)) .lt. 1.d-12) then
      n_good += 1
    else
      print*,'Found an imaginary component to eigenvalue'
      print*,'Re(i) + Im(i)',WR(i),WI(i)
    endif
  enddo

  allocate(list_good(n_good), iorder(n_good))
  n_good = 0
  do i = 1, n
    if(dabs(WI(i)).lt.1.d-12)then
      n_good += 1
      list_good(n_good) = i
      eigval(n_good) = WR(i)
    endif
  enddo
  n_real_eigv = n_good 
  do i = 1, n_good
   iorder(i) = i
  enddo

  ! You sort the real eigenvalues 
  call dsort(eigval, iorder, n_good)
  print*,'n_real_eigv = ', n_real_eigv
  print*,'n           = ', n
  do i = 1, n_real_eigv
    print*,i,'eigval(i) = ', eigval(i) 
    do j = 1, n
      reigvec(j,i) = VR(j,list_good(iorder(i)))
      leigvec(j,i) = Vl(j,list_good(iorder(i)))
    enddo
  enddo

end subroutine non_hrmt_general_real_diag

! ---

subroutine impose_biorthog_qr(m, n, thr_d, thr_nd, Vl, Vr)

  implicit none 
  integer,          intent(in)    :: m, n
  double precision, intent(in)    :: thr_d, thr_nd
  double precision, intent(inout) :: Vl(m,n), Vr(m,n)

  integer                         :: i, j
  integer                         :: LWORK, INFO
  double precision                :: accu_nd, accu_d
  double precision, allocatable   :: TAU(:), WORK(:)
  double precision, allocatable   :: S(:,:), R(:,:), tmp(:,:)

  ! ---

  call check_biorthog_binormalize(m, n, Vl, Vr, thr_d, thr_nd, .false.)
  
  ! ---
  
  allocate(S(n,n))
  call dgemm( 'T', 'N', n, n, m, 1.d0          &
            , Vl, size(Vl, 1), Vr, size(Vr, 1) &
            , 0.d0, S, size(S, 1) )

  accu_nd = 0.d0
  accu_d  = 0.d0
  do i = 1, n
    do j = 1, n
      if(i==j) then
        accu_d += S(j,i)
      else
        accu_nd = accu_nd + S(j,i) * S(j,i)
      endif
    enddo
  enddo
  accu_nd = dsqrt(accu_nd)

  if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n))/dble(n) .lt. thr_d)) then
    print *, ' bi-orthogonal vectors without QR !'
    deallocate(S)
    return
  endif

  ! -------------------------------------------------------------------------------------
  !                           QR factorization of S: S = Q x R


  print *, ' apply QR decomposition ...'

  allocate( TAU(n), WORK(1) )

  LWORK = -1
  call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO)
  if(INFO .ne. 0) then
    print*,'dgeqrf failed !!', INFO
    stop
  endif

  LWORK = max(n, int(WORK(1)))
  deallocate(WORK)

  allocate( WORK(LWORK) )
  call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO)
  if(INFO .ne. 0) then
    print*,'dgeqrf failed !!', INFO
    stop
  endif

  ! save the upper triangular R
  allocate( R(n,n) )
  R(:,:) = S(:,:)

  ! get Q
  LWORK = -1
  call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO)
  if(INFO .ne. 0) then
    print*,'dorgqr failed !!', INFO
    stop
  endif

  LWORK = max(n, int(WORK(1)))
  deallocate(WORK)

  allocate( WORK(LWORK) )
  call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO)
  if(INFO .ne. 0) then
    print*,'dorgqr failed !!', INFO
    stop
  endif

  deallocate( WORK, TAU )

  !
  ! -------------------------------------------------------------------------------------

  ! ---

  ! -------------------------------------------------------------------------------------
  !                               get bi-orhtog left & right vectors:
  !                                           Vr' = Vr x inv(R) 
  !                                           Vl' = inv(Q) x Vl =  Q.T   x Vl 

  ! Q.T x Vl, where Q = S

  allocate( tmp(n,m) )
  call dgemm( 'T', 'T', n, m, n, 1.d0        &
            , S, size(S, 1), Vl, size(Vl, 1) &
            , 0.d0, tmp, size(tmp, 1) )

  do i = 1, n
    do j = 1, m
      Vl(j,i) = tmp(i,j)
    enddo
  enddo
  deallocate(tmp)

  ! ---

  ! inv(R) 
  !print *, ' inversing upper triangular matrix ...'
  call dtrtri("U", "N", n, R, n, INFO)
  if(INFO .ne. 0) then
    print*,'dtrtri failed !!', INFO
    stop
  endif
  !print *, ' inversing upper triangular matrix OK' 

  do i = 1, n-1
    do j = i+1, n
      R(j,i) = 0.d0
    enddo
  enddo

  !print *, ' inv(R):'
  !do i = 1, n
  !  write(*, '(1000(F16.10,X))') R(i,:)
  !enddo

  ! Vr x inv(R) 
  allocate( tmp(m,n) )
  call dgemm( 'N', 'N', m, n, n, 1.d0        &
            , Vr, size(Vr, 1), R, size(R, 1) &
            , 0.d0, tmp, size(tmp, 1) )
  deallocate( R )

  do i = 1, n
    do j = 1, m
      Vr(j,i) = tmp(j,i)
    enddo
  enddo
  deallocate(tmp)

  return
end subroutine impose_biorthog_qr

! ---

subroutine impose_biorthog_lu(m, n, Vl, Vr, S)

  implicit none 
  integer, intent(in)             :: m, n
  double precision, intent(inout) :: Vl(m,n), Vr(m,n), S(n,n)

  integer                         :: i, j
  integer                         :: INFO
  double precision                :: nrm
  integer,          allocatable   :: IPIV(:)
  double precision, allocatable   :: L(:,:), tmp(:,:), vectmp(:)
  !double precision, allocatable   :: T(:,:), ll(:,:), rr(:,:), tt(:,:)

  !allocate( T(n,n) )
  !T(:,:) = S(:,:)

  print *, ' apply LU decomposition ...'

  ! -------------------------------------------------------------------------------------
  !                           LU factorization of S: S = P x L x U

  allocate( IPIV(n) )

  call dgetrf(n, n, S, n, IPIV, INFO)
  if(INFO .ne. 0) then
    print*, 'dgetrf failed !!', INFO
    stop
  endif

  ! check | S - P x L x U |
  !allocate( ll(n,n), rr(n,n), tmp(n,n) )
  !ll = S
  !rr = S
  !do i = 1, n-1
  !  ll(i,i) = 1.d0
  !  do j = i+1, n
  !    ll(i,j) = 0.d0
  !    rr(j,i) = 0.d0
  !  enddo
  !enddo
  !ll(n,n) = 1.d0
  !call dgemm( 'N', 'N', n, n, n, 1.d0          &
  !          , ll, size(ll, 1), rr, size(rr, 1) &
  !          , 0.d0, tmp, size(tmp, 1) )
  ! deallocate(ll, rr)
  !allocate( vectmp(n) )
  !do j = n-1, 1, -1
  !  i = IPIV(j)
  !  if(i.ne.j) then
  !    print *, j, i
  !    vectmp(:) = tmp(i,:)
  !    tmp(i,:)  = tmp(j,:)
  !    tmp(j,:)  = vectmp(:)
  !  endif
  !enddo
  !deallocate( vectmp )
  !nrm = 0.d0
  !do i = 1, n
  !  do j = 1, n
  !    nrm += dabs(tmp(j,i) - T(j,i))
  !  enddo
  !enddo
  !deallocate( tmp )
  !print*, '|L.T x R - S| =', nrm
  !stop

  ! ------
  ! inv(L) 
  ! ------

  allocate( L(n,n) )
  L(:,:) = S(:,:)

  call dtrtri("L", "U", n, L, n, INFO)
  if(INFO .ne. 0) then
    print*,  'dtrtri failed !!', INFO
    stop
  endif
  do i = 1, n-1
    L(i,i) = 1.d0
    do j = i+1, n
      L(i,j) = 0.d0
    enddo
  enddo
  L(n,n) = 1.d0

  ! ------
  ! inv(U) 
  ! ------
  
  call dtrtri("U", "N", n, S, n, INFO)
  if(INFO .ne. 0) then
    print*,  'dtrtri failed !!', INFO
    stop
  endif

  do i = 1, n-1
    do j = i+1, n
      S(j,i) = 0.d0
    enddo
  enddo

  !
  ! -------------------------------------------------------------------------------------

  ! ---

  ! -------------------------------------------------------------------------------------
  !                               get bi-orhtog left & right vectors:
  !                                           Vr' = Vr x inv(U) 
  !                                           Vl' = inv(L) x inv(P) x Vl

  ! inv(P) x Vl
  allocate( vectmp(n) )
  do j = n-1, 1, -1
    i = IPIV(j)
    if(i.ne.j) then
      vectmp(:) = L(:,j)
      L(:,j)    = L(:,i)
      L(:,i)    = vectmp(:)
    endif
  enddo
  deallocate( vectmp )

  ! Vl'
  allocate( tmp(m,n) )
  call dgemm( 'N', 'T', m, n, n, 1.d0        &
            , Vl, size(Vl, 1), L, size(L, 1) &
            , 0.d0, tmp, size(tmp, 1) )
  deallocate(L)

  Vl = tmp
  deallocate(tmp)

  ! ---

  ! Vr x inv(U) 
  allocate( tmp(m,n) )
  call dgemm( 'N', 'N', m, n, n, 1.d0        &
            , Vr, size(Vr, 1), S, size(S, 1) &
            , 0.d0, tmp, size(tmp, 1) )
  Vr = tmp
  deallocate(tmp)

  !allocate( tmp(n,n) )
  !call dgemm( 'T', 'N', n, n, m, 1.d0          &
  !          , Vl, size(Vl, 1), Vr, size(Vr, 1) &
  !          , 0.d0, tmp, size(tmp, 1) )
  !nrm = 0.d0
  !do i = 1, n
  !  do j = 1, n
  !    nrm += dabs(tmp(j,i))
  !  enddo
  !enddo
  !deallocate( tmp )
  !print*, '|L.T x R| =', nrm
  !stop

  return
end subroutine impose_biorthog_lu

! ---

subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, stop_ifnot)

  implicit none
  integer,          intent(in)  :: n, m
  logical,          intent(in)  :: stop_ifnot
  double precision, intent(in)  :: A(n,n), eigval(m), leigvec(n,m), reigvec(n,m), thr_diag, thr_norm
 
  integer                       :: i, j
  double precision              :: tmp, tmp_abs, tmp_nrm, tmp_rel, tmp_dif
  double precision              :: V_nrm, U_nrm
  double precision, allocatable :: Mtmp(:,:)

  allocate( Mtmp(n,m) )
  
  ! ---

  Mtmp = 0.d0
  call dgemm( 'N', 'N', n, m, n, 1.d0                  &
            , A, size(A, 1), reigvec, size(reigvec, 1) &
            , 0.d0, Mtmp, size(Mtmp, 1) )

  V_nrm   = 0.d0
  tmp_nrm = 0.d0
  tmp_abs = 0.d0
  do j = 1, m
    
    tmp   = 0.d0
    U_nrm = 0.d0
    do i = 1, n
      tmp     = tmp     + dabs(Mtmp(i,j) - eigval(j) * reigvec(i,j))
      tmp_nrm = tmp_nrm + dabs(Mtmp(i,j))
      U_nrm   = U_nrm   + reigvec(i,j) * reigvec(i,j)
    enddo

    tmp_abs = tmp_abs + tmp
    V_nrm   = V_nrm   + U_nrm 
    !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm

  enddo

  if(tmp_abs.lt.10.d-10)then
   tmp_rel = thr_diag/10.d0
  else
   tmp_rel = tmp_abs / tmp_nrm
  endif
  tmp_dif = dabs(V_nrm - dble(m))

  if( stop_ifnot .and. ((tmp_rel .gt. thr_diag) .or. (tmp_dif .gt. thr_norm)) ) then
    print *, ' error in right-eigenvectors'
    print *, ' err tol   = ',thr_diag, thr_norm
    print *, '(tmp_rel .gt. thr_diag) = ',(tmp_rel .gt. thr_diag)
    print *, '(tmp_dif .gt. thr_norm) = ',(tmp_dif .gt. thr_norm)
    print *, ' err estim = ', tmp_abs, tmp_rel
    print *, ' CR norm   = ', V_nrm 
    stop
  endif

  ! ---

  Mtmp = 0.d0
  call dgemm( 'T', 'N', n, m, n, 1.d0                  &
            , A, size(A, 1), leigvec, size(leigvec, 1) &
            , 0.d0, Mtmp, size(Mtmp, 1) )

  V_nrm   = 0.d0
  tmp_nrm = 0.d0
  tmp_abs = 0.d0
  do j = 1, m

    tmp   = 0.d0
    U_nrm = 0.d0
    do i = 1, n
      tmp     = tmp     + dabs(Mtmp(i,j) - eigval(j) * leigvec(i,j))
      tmp_nrm = tmp_nrm + dabs(Mtmp(i,j))
      U_nrm   = U_nrm   + leigvec(i,j) * leigvec(i,j)
    enddo

    tmp_abs = tmp_abs + tmp
    V_nrm   = V_nrm   + U_nrm 
    !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm

  enddo

  if(tmp_abs.lt.10.d-10)then
   tmp_rel = thr_diag/10.d0
  else
   tmp_rel = tmp_abs / tmp_nrm
  endif
  if( stop_ifnot .and. ((tmp_rel .gt. thr_diag) .or. (tmp_dif .gt. thr_norm)) ) then
    print *, ' error in left-eigenvectors'
    print *, ' err tol   = ',thr_diag, thr_norm
    print *, '(tmp_rel .gt. thr_diag) = ',(tmp_rel .gt. thr_diag)
    print *, '(tmp_dif .gt. thr_norm) = ',(tmp_dif .gt. thr_norm)
    print *, ' err estim = ', tmp_abs, tmp_rel
    print *, ' CR norm   = ', V_nrm 
    stop
  endif

  ! ---

  deallocate( Mtmp )

end subroutine check_EIGVEC

! ---

subroutine check_degen(n, m, eigval, leigvec, reigvec)

  implicit none
  integer,          intent(in)    :: n, m
  double precision, intent(in)    :: eigval(m)
  double precision, intent(inout) :: leigvec(n,m), reigvec(n,m)
 
  integer                         :: i, j
  double precision                :: ei, ej, de, de_thr, accu_nd
  double precision, allocatable   :: S(:,:)

  de_thr = 1d-6

  do i = 1, m-1
    ei = eigval(i)

    do j = i+1, m
      ej = eigval(j)
      de = dabs(ei - ej)

      if(de .lt. de_thr) then

        leigvec(:,i) = 0.d0
        leigvec(:,j) = 0.d0
        leigvec(i,i) = 1.d0
        leigvec(j,j) = 1.d0

        reigvec(:,i) = 0.d0
        reigvec(:,j) = 0.d0
        reigvec(i,i) = 1.d0
        reigvec(j,j) = 1.d0

      endif

    enddo
  enddo

  ! ---

  allocate( S(m,m) )

  ! S = VL x VR
  call dgemm( 'T', 'N', m, m, n, 1.d0                              &
            , leigvec, size(leigvec, 1), reigvec, size(reigvec, 1) &
            , 0.d0, S, size(S, 1) )

  accu_nd = 0.d0
  do i = 1, m
    do j = 1, m
      if(i==j) cycle
      accu_nd = accu_nd + S(j,i) * S(j,i)
    enddo
  enddo
  accu_nd = dsqrt(accu_nd)

  deallocate( S )

  print *, ' check_degen: L & T bi-orthogonality: ok'
  print *, ' accu_nd = ', accu_nd