program f90_JACOBI
    implicit none
    integer::i,j,n,nt
    real(8),allocatable::ary0(:,:),ary1(:,:)
    real(8),allocatable::ev1(:),vec1(:,:)
    real(8),allocatable::ev2(:),vec2(:,:),b(:,:)
    character::linebuf*1000,dummy*100

    nt=6
    n=nt
    allocate(ary0(1:nt,1:nt))
    allocate(ary1(1:nt,1:nt))
    allocate(ev1(1:nt))
    allocate(vec1(1:nt,1:nt))
    allocate(ev2(1:nt))
    allocate(vec2(1:nt,1:nt))
    allocate(b(1:nt,1:nt))

!    ary0(1, 1) = 4 ; ary0(1, 2) = 2 ; ary0(1, 3) = 3
!    ary0(2, 1) = 2 ; ary0(2, 2) = 6 ; ary0(2, 3) = 1
!    ary0(3, 1) = 3 ; ary0(3, 2) = 1 ; ary0(3, 3) = 5
    !*** Eigen value ***
    !9.000000e+000   4.732051e+000   1.267949e+000
    !*** Eigen vector ***
    !9.227825e-001  -2.633891e-001   1.060632e+000
    !9.227825e-001   9.829816e-001  -2.841954e-001
    !9.227825e-001  -7.195925e-001  -7.764364e-001

    !ヒルベルト行列
    do i=1,n
        do j=1,n
            ary0(i,j)=1.0D0/(dble(i)+dble(j)-1.0D0)
        end do
    end do

    do i=1,n
        do j=1,n
            ary1(i,j)=ary0(i,j)
        end do
    end do
    call JACOBI(n,ary1,ev1,vec1,nt)

    do i=1,n
        do j=1,n
            ary1(i,j)=ary0(i,j)
            b(i,j)=0.0D0
        end do
        b(i,i)=1.0D0
    end do
    call GJACOBI(n,ary1,b,ev2,vec2,nt)

    write(dummy,*) ev1(1)
    linebuf=trim(dummy)
    do j=2,n
        write(dummy,*) ev1(j)
        linebuf=trim(linebuf)//','//trim(dummy)
    end do
    call del_spaces(linebuf)
    write(6,*) trim(linebuf)
    do i=1,n
        write(dummy,*) vec1(i,1)
        linebuf=trim(dummy)
        do j=2,n
            write(dummy,*) vec1(i,j)
            linebuf=trim(linebuf)//','//trim(dummy)
        end do
        call del_spaces(linebuf)
        write(6,*) trim(linebuf)
    end do
    write(6,'(a)') '****'

    write(dummy,*) ev2(1)
    linebuf=trim(dummy)
    do j=2,n
        write(dummy,*) ev2(j)
        linebuf=trim(linebuf)//','//trim(dummy)
    end do
    call del_spaces(linebuf)
    write(6,*) trim(linebuf)
    do i=1,n
        write(dummy,*) vec2(i,1)
        linebuf=trim(dummy)
        do j=2,n
            write(dummy,*) vec2(i,j)
            linebuf=trim(linebuf)//','//trim(dummy)
        end do
        call del_spaces(linebuf)
        write(6,*) trim(linebuf)
    end do

stop

contains

    subroutine JACOBI(n,array,ev,vec,nt)
        !Jacobi法による固有値・固有ベクトル計算（奥村先生）
        !固有値ev(k),固有ベクトルvec(k,i)

        integer,intent(in)::n,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        real(8),intent(out)::ev(1:nt),vec(1:nt,1:nt)

        integer::i,j,k,im
        real(8)::s,c,aoff,eps,t
        real(8)::aa,aj,ak,vj,vk,vmax,bb

        do i=1,n
            do j=1,n
                vec(i,j)=0.0
            end do
            vec(i,i)=1.0
        end do

        aoff=0.0
        do j=1,n-1
            do k=j+1,n
                aoff=aoff+array(j,k)*array(j,k)
            end do
        end do
        eps=aoff*2.0
        do j=1,n
            eps=eps+array(j,j)*array(j,j)
        end do
        eps=eps*0.00000000005

        do while(aoff>eps)
            do j=1,n-1
                do k=j+1,n
                    if(array(j,k)==0.0)then
                        t=0.0
                    else
                        aa=(array(k,k)-array(j,j))/(2.0*array(j,k))
                    end if
                    if(aa>=0)then
                        t=1.0/(aa+sqrt(1.0+aa*aa))
                    else
                        t=1.0/(aa-sqrt(1.0+aa*aa))
                    end if
                    c=1.0/sqrt(1.0+t*t)
                    s=t*c
                    array(j,j)=array(j,j)-t*array(j,k)
                    array(k,k)=array(k,k)+t*array(j,k)
                    array(j,k)=0.0
                    do i=1,j-1
                        aj=array(i,j)
                        ak=array(i,k)
                        array(i,j)=aj*c-ak*s
                        array(i,k)=aj*s+ak*c
                    end do
                    do i=j+1,k-1
                        aj=array(j,i)
                        ak=array(i,k)
                        array(j,i)=aj*c-ak*s
                        array(i,k)=aj*s+ak*c
                    end do
                    do i=k+1,n
                        aj=array(j,i)
                        ak=array(k,i)
                        array(j,i)=aj*c-ak*s
                        array(k,i)=aj*s+ak*c
                    end do
                    do i=1,n
                        vj=vec(i,j)
                        vk=vec(i,k)
                        vec(i,j)=vj*c-vk*s
                        vec(i,k)=vj*s+vk*c
                    end do
                end do
            end do
            aoff=0.0
            do j=1,n-1
                do k=j+1,n
                    aoff=aoff+array(j,k)*array(j,k)
                end do
            end do
        end do

        do k=1,n
            ev(k)=array(k,k)
        end do

        do k=1,n-1
            vmax=-1.0E+30
            do i=k,n
                if(ev(i)>=vmax)then
                    im=i
                    vmax=ev(i)
                end if
            end do
            aa=ev(k)
            ev(k)=ev(im)
            ev(im)=aa
            do i=1,n
                bb=vec(i,k)
                vec(i,k)=vec(i,im)
                vec(i,im)=bb
            end do
        end do
    endsubroutine JACOBI


    subroutine GJACOBI(nf,a,b,ev,vec,nt)
        !************************************
        !一般化JACOBI法による固有値計算
        !(パソコン構造力学，宮澤健二著)
        !************************************
        integer,intent(in)::nf,nt
        real(8),intent(inout)::a(1:nt,1:nt),b(1:nt,1:nt)
        real(8),intent(out)::ev(1:nt),vec(1:nt,1:nt)

        integer::i,j,k
        real(8)::fmk,fmg,epsk,epsg
        real(8)::fact=1.0D-30
        integer::kmax,ip,iq,im
        real(8)::amax,dam,de1,de2,de3,dd,ss,x,gg,ww
        real(8)::amin

        fmk=0.0D0
        fmg=0.0D0
        do i=1,nf
            do j=1,nf
                if(abs(a(i,j))>fmk)fmk=abs(a(i,j))
                if(abs(b(i,j))>fmg)fmg=abs(b(i,j))
            end do
        end do
        epsk=fmk*fact
        epsg=fmg*fact
        do i=1,nf
            do j=1,nf
                vec(i,j)=0.0
            end do
            vec(i,i)=1.0
        end do
        kmax=5*nf*nf
        do k=1,kmax
            ip=1
            iq=2
            amax=a(1,2)**2.0D0/(a(1,1)**2.0D0+a(2,2)**2.0D0)
            do i=1,nf-1
                do j=i+1,nf
                    dam=a(i,j)**2.0D0/(a(i,i)**2.0D0+a(j,j)**2.0D0)
                    if(amax<dam)then
                        ip=i
                        iq=j
                        amax=dam
                    end if
                end do
            end do
            do i=1,nf-1
                do j=i+1,nf
                    dam=b(i,j)**2.0D0/(b(i,i)**2.0D0+b(j,j)**2.0D0)
                    if(amax<dam)then
                        ip=i
                        iq=j
                        amax=dam
                    end if
                end do
            end do
            de1=a(ip,ip)*b(ip,iq)-a(ip,iq)*b(ip,ip)
            de2=a(iq,iq)*b(ip,iq)-a(ip,iq)*b(iq,iq)
            de3=a(ip,ip)*b(iq,iq)-a(iq,iq)*b(ip,ip)
            dd=de3*de3+4.0D0*de1*de2
            ss=sqrt(dd)
            if(de3<0.0D0)then
                x=0.5D0*(de3-ss)
            else
                x=0.5D0*(de3+ss)
            end if
            if(x==0.0D0)then
                dam=0.0D0
                gg=-a(ip,iq)/a(iq,iq)
            else
                dam=de2/x
                gg=-de1/x
            end if
            do j=1,nf
                ww=a(ip,j)
                a(ip,j)=a(ip,j)+gg*a(iq,j)
                a(iq,j)=a(iq,j)+dam*ww
                ww=b(ip,j)
                b(ip,j)=b(ip,j)+gg*b(iq,j)
                b(iq,j)=b(iq,j)+dam*ww
            end do
            do i=1,nf
                ww=a(i,ip)
                a(i,ip)=a(i,ip)+gg*a(i,iq)
                a(i,iq)=a(i,iq)+dam*ww
                ww=b(i,ip)
                b(i,ip)=b(i,ip)+gg*b(i,iq)
                b(i,iq)=b(i,iq)+dam*ww
                ww=vec(i,ip)
                vec(i,ip)=vec(i,ip)+gg*vec(i,iq)
                vec(i,iq)=vec(i,iq)+dam*ww
            end do
            if(amax<epsk.and.amax<epsg)exit
        end do
        do i=1,nf
            ev(i)=a(i,i)/b(i,i)
        end do

        !並び替え
        !ev(k)：最小固有値k=1，最大固有値k=nf
        !vec(i,k)：k番目の固有値に対応する固有ベクトル
        do k=1,nf-1
            amin=1.0E+30
            do i=k,nf
                if(ev(i)<=amin)then
                    im=i
                    amin=ev(i)
                end if
            end do
            ww=ev(k)
            ev(k)=ev(im)
            ev(im)=ww
            do i=1,nf
                ww=vec(i,k)
                vec(i,k)=vec(i,im)
                vec(i,im)=ww
            end do
        end do
    end subroutine GJACOBI


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


end program f90_JACOBI