! ヤコビ法のサブルーチンの使用例
implicit none
real(8),allocatable :: a(:,:)
real(8),allocatable :: v(:,:)
integer :: n=4    !行列の大きさ
integer :: i,j,k,l
allocate(a(n,n),v(n,n))

! 例(実対称行列でないと駄目)
a(1,:)=(/ 1. , 1. , 1. , 1. /)
a(2,:)=(/ 1. , 1. , 1. , 1. /)
a(3,:)=(/ 1. , 1. , 1. , 1. /)
a(4,:)=(/ 1. , 1. , 1. , 1. /)

print*,'与えた実対称行列'
do i=1,n
   print*,a(i,:)
enddo

call eigen_jacobi(n,a,v)

print*,'固有値（対角行列）'
do i=1,n
   print*,a(i,:)
enddo

print*,'対応する固有ベクトル'
do i=1,n
   print*,v(i,:)
enddo

end

! 実対称行列の固有値をヤコビ法で求めるサブルーチン
! 戸川隼人『詳解 数値計算演習』(共立出版)p.163-169 を参考にした．
subroutine eigen_jacobi(n,a,v)
  implicit none
  integer,intent(in)    :: n      ! 行列の大きさ
  real(8),intent(inout) :: a(n,n) ! 対象とする実対称行列
  real(8),intent(out)   :: v(n,n) ! 固有ベクトル
  real(8) :: a1(n,n)              ! a_k+1
  real(8) :: v1(n,n)              ! v_k+1
  real(8) :: err=1.d-8            ! 許容誤差
  real(8) :: amax,th,c,s
  integer :: i,j,k,l

  a1=a
  ! 初期の固有ベクトルの単位行列化
  v=0.d0
  do i=1,n
     v(i,i)=1.d0
  enddo
  v1=v

  ! 反復計算のループ
  do
    ! 非対角行列の中で絶対値が最大なものを探す
    amax=0.d0
    do j=2,n
       do i=1,j-1
          if (abs(a(i,j))>amax) then
             amax=abs(a(i,j))
             k=i
             l=j
          endif
       enddo
    enddo
    if (amax**2<err) exit ! 収束判定(誤差はamax^2のオーダー)

    ! 行列P中のsin(th),cos(th)
    th=atan(2.d0*a(k,l)/(a(k,k)-a(l,l)))/2
    c=cos(th)
    s=sin(th)

    ! PAP^tの計算(効率的にk,l行とk,l列だけを計算する)
    a1(k,:)=a(k,:)*c+a(l,:)*s
    a1(l,:)=a(l,:)*c-a(k,:)*s
    a1(:,k)=a1(k,:)
    a1(:,l)=a1(l,:)
    a1(k,k)=a(k,k)*c**2+2.d0*a(k,l)*c*s+a(l,l)*s**2
    a1(l,l)=a(l,l)*c**2-2.d0*a(k,l)*c*s+a(k,k)*s**2
    a1(k,l)=a(k,l)*c**2+(a(l,l)-a(k,k))*c*s-a(k,l)*s**2
    a1(l,k)=a(k,l)*c**2+(a(l,l)-a(k,k))*c*s-a(k,l)*s**2
    a=a1 ! Aの更新

    ! VP^tの計算(効率的にi,j列だけを計算する)
    v1(:,k)=v(:,k)*c+v(:,l)*s
    v1(:,l)=v(:,l)*c-v(:,k)*s
    v=v1 ! Vの更新 

  enddo
end
