!
!     2次精度の差分法を用いた
!     ベナール対流の成長率の計算(rigid-rigid, T-T)。
!
program SolveEigenvalueProblemOfZBenardConvection
  
  implicit none
  
  integer, parameter :: m = 41
  
  real(8) al, R, Pr
  integer i,j,ii,jj
  real(8) al0, al1, dal, R0, R1, dR
  complex(8) ca(1:2*m-2, 1:2*m-2), cb(1:2*m-2, 1:2*m-2), eigenvector(1:2*m-2), eigenvalue

  Pr = 1.0d0

  al0 = 0.0d0
  al1 = 8.0d0
  dal = 0.1d0
  R0 = 0.0d0
  R1 = 10000.0d0
  dR  = 100.0d0

  ii = int((al1-al0)/dal)
  jj = int((R1-R0)/dR)

  print *, "# k, Ra, eigenvalue(real), eigenvalue(imag)"
  do j=0,jj
    R = R0 + dR*dble(j)
    do i=0,ii
      al = al0 + dal*dble(i)
      call calculate_matrices(m, al, R, Pr, ca, cb)
      call solve_eigenvalue_problem(2*m-2, ca, cb, eigenvector, eigenvalue)
      print *, al, R, dble(eigenvalue), dimag(eigenvalue)
    end do
  end do
  
  stop
end program SolveEigenvalueProblemOfZBenardConvection

! 
! 一般化固有値問題 $A x = s B x$ に現れる行列 A, Bを計算する
! 
!  ここでそれぞれのブロックはほぼ対角成分のみの行列を想定する。
!  ( Pr*(D^2-a^2)^2  -a^2*R*Pr)(W) = s*(D^2-a^2   0)(W)
!  ( 1               D^2-a^2  )(T)     (0         1)(T)
!
subroutine calculate_matrices(m, al, R, Pr, ca, cb)
  
  implicit none

  integer, intent(in) :: m
  real(8), intent(in) :: al, R, Pr
  complex(8), intent(out) :: ca(1:2*m-2, 1:2*m-2), cb(1:2*m-2, 1:2*m-2)

  integer i
  real(8) dz

  dz = 1.0d0/dble(m)

  ca(1:2*m-2, 1:2*m-2) = 0.0d0
  cb(1:2*m-2, 1:2*m-2) = 0.0d0

  i=1
  ca(i,i)   = ca(i,i) + Pr*1.0d0/dz**4
  ca(i,i)   = ca(i,i) + Pr*(6.0d0/dz**4+4.0d0*al**2/dz**2+al**4)
  ca(i+1,i) = ca(i+1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)
  ca(i+2,i) = ca(i+2,i) + Pr*1.0d0/dz**4

  ca(i,i+m-1) = ca(i,i+m-1) + 1.0d0

  ca(i+m-1,i) = ca(i+m-1,i) -al**2*R*Pr

  ca(i+m-1,i+m-1) = ca(i+m-1,i+m-1) -2.0d0/dz**2-al**2
  ca(i+1+m-1,i+m-1) = ca(i+1+m-1,i+m-1) + 1.0d0/dz**2

  cb(i,i) = cb(i,i) -2.0d0/dz**2 -al**2
  cb(i+1,i) = cb(i+1,i) + 1.0d0/dz**2

  cb(i+m-1,i) = cb(i+m-1,i) + 0.0d0

  cb(i,i+m-1) = cb(i,i+m-1) + 0.0d0

  cb(i+m-1,i+m-1) = cb(i+m-1,i+m-1) + 1.0d0

  i=2
  ca(i-1,i) = ca(i-1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)
  ca(i,i)   = ca(i,i) + Pr*(6.0d0/dz**4+4.0d0*al**2/dz**2+al**4)
  ca(i+1,i) = ca(i+1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)
  ca(i+2,i) = ca(i+2,i) + Pr*1.0d0/dz**4

  ca(i,i+m-1) = ca(i,i+m-1) + 1.0d0

  ca(i+m-1,i) = ca(i+m-1,i) -al**2*R*Pr

  ca(i-1+m-1,i+m-1) = ca(i-1+m-1,i+m-1) + 1.0d0/dz**2
  ca(i+m-1,i+m-1) = ca(i+m-1,i+m-1) -2.0d0/dz**2-al**2
  ca(i+1+m-1,i+m-1) = ca(i+1+m-1,i+m-1) + 1.0d0/dz**2

  cb(i-1,i) = cb(i-1,i) + 1.0d0/dz**2
  cb(i,i) = cb(i,i) -2.0d0/dz**2 -al**2
  cb(i+1,i) = cb(i+1,i) + 1.0d0/dz**2

  cb(i+m-1,i) = cb(i+m-1,i) + 0.0d0

  cb(i,i+m-1) = cb(i,i+m-1) + 0.0d0

  cb(i+m-1,i+m-1) = cb(i+m-1,i+m-1) + 1.0d0

  do i=3,m-3
    ca(i-2,i) = ca(i-2,i) + Pr*1.0d0/dz**4
    ca(i-1,i) = ca(i-1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)
    ca(i,i)   = ca(i,i) + Pr*(6.0d0/dz**4+4.0d0*al**2/dz**2+al**4)
    ca(i+1,i) = ca(i+1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)
    ca(i+2,i) = ca(i+2,i) + Pr*1.0d0/dz**4

    ca(i,i+m-1) = ca(i,i+m-1) + 1.0d0

    ca(i+m-1,i) = ca(i+m-1,i) -al**2*R*Pr

    ca(i-1+m-1,i+m-1) = ca(i-1+m-1,i+m-1) + 1.0d0/dz**2
    ca(i+m-1,i+m-1) = ca(i+m-1,i+m-1) -2.0d0/dz**2-al**2
    ca(i+1+m-1,i+m-1) = ca(i+1+m-1,i+m-1) + 1.0d0/dz**2

    cb(i-1,i) = cb(i-1,i) + 1.0d0/dz**2
    cb(i,i) = cb(i,i) -2.0d0/dz**2 -al**2
    cb(i+1,i) = cb(i+1,i) + 1.0d0/dz**2

    cb(i+m-1,i) = cb(i+m-1,i) + 0.0d0

    cb(i,i+m-1) = cb(i,i+m-1) + 0.0d0

    cb(i+m-1,i+m-1) = cb(i+m-1,i+m-1) + 1.0d0
  end do
  i=m-2
  ca(i-2,i) = ca(i-2,i) + Pr*1.0d0/dz**4
  ca(i-1,i) = ca(i-1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)
  ca(i,i)   = ca(i,i) + Pr*(6.0d0/dz**4+4.0d0*al**2/dz**2+al**4)
  ca(i+1,i) = ca(i+1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)

  ca(i,i+m-1) = ca(i,i+m-1) + 1.0d0

  ca(i+m-1,i) = ca(i+m-1,i) -al**2*R*Pr

  ca(i-1+m-1,i+m-1) = ca(i-1+m-1,i+m-1) + 1.0d0/dz**2
  ca(i+m-1,i+m-1) = ca(i+m-1,i+m-1) -2.0d0/dz**2-al**2
  ca(i+1+m-1,i+m-1) = ca(i+1+m-1,i+m-1) + 1.0d0/dz**2

  cb(i-1,i) = cb(i-1,i) + 1.0d0/dz**2
  cb(i,i) = cb(i,i) -2.0d0/dz**2 -al**2
  cb(i+1,i) = cb(i+1,i) + 1.0d0/dz**2

  cb(i+m-1,i) = cb(i+m-1,i) + 0.0d0

  cb(i,i+m-1) = cb(i,i+m-1) + 0.0d0

  cb(i+m-1,i+m-1) = cb(i+m-1,i+m-1) + 1.0d0

  i=m-1
  ca(i-2,i) = ca(i-2,i) + Pr*1.0d0/dz**4
  ca(i-1,i) = ca(i-1,i) + Pr*(-4.0d0/dz**4-2.0d0*al**2/dz**2)
  ca(i,i)   = ca(i,i) + Pr*(6.0d0/dz**4+4.0d0*al**2/dz**2+al**4)
  ca(i,i) = ca(i,i) + Pr*1.0d0/dz**4

  ca(i,i+m-1) = ca(i,i+m-1) + 1.0d0

  ca(i+m-1,i) = ca(i+m-1,i) -al**2*R*Pr

  ca(i-1+m-1,i+m-1) = ca(i-1+m-1,i+m-1) + 1.0d0/dz**2
  ca(i+m-1,i+m-1) = ca(i+m-1,i+m-1) -2.0d0/dz**2-al**2

  cb(i-1,i) = cb(i-1,i) + 1.0d0/dz**2
  cb(i,i) = cb(i,i) -2.0d0/dz**2 -al**2

  cb(i+m-1,i) = cb(i+m-1,i) + 0.0d0

  cb(i,i+m-1) = cb(i,i+m-1) + 0.0d0

  cb(i+m-1,i+m-1) = cb(i+m-1,i+m-1) + 1.0d0

end subroutine calculate_matrices


! 
! 固有値問題を解き, 第一不安定モードまたは最小減衰モードと, それに対応する
! 固有関数を求める
! 
! Note:
!    非対称行列の一般化固有値問題を解くLapackのサブルーチン zgegv を用いる.
!
!    zgegvの使い方は, 
!    call zgegv( 'N', 'V', N, A, LDA, B, LDB, ALPHA, BETA,
!                VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
!
!    LDA, LDBは leading dimension of A, B
!    固有値 lambda は lambda = alpha/beta の形で返る.
!    VLは左固有ベクトル, VRは右固有ベクトル
!    詳しい使い方は http://www.netlib.org/lapack/complex16/zgegv.f を参照すること.
!
subroutine solve_eigenvalue_problem(n, ca, cb, eigenvector, eigenvalue)

  implicit none

  integer, intent(in) :: n
  complex(8), intent(inout) :: ca(1:n, 1:n), cb(1:n, 1:n)
  complex(8), intent(out) :: eigenvector(1:n), eigenvalue

  integer info
  complex(8) alpha(1:n), beta(1:n), dummy(1:1, 1:n), rev(1:n, 1:n), work(1:24*n)
  real(8) work2(1:8*n)
  
  complex(8) lambda(1:n)
  integer i_max

  call zgegv( 'N', 'V', n, ca, n, cb, n, alpha, beta, &
       &      dummy, 1, rev, n, work, 24*n, work2, info )
  if( info /= 0 )then
     print *, 'error has occured in zgegv: error_code=', info
     call exit(1)
  end if
  lambda(1:n) = alpha(1:n)/beta(1:n)
  
  i_max = maxloc(dble(lambda), 1)

  eigenvector(1:n) = rev(1:n, i_max)
  eigenvalue = lambda(i_max)

end subroutine solve_eigenvalue_problem
