program f90_GJ
    implicit none
    integer::i,j,n,m,nt
    real(8),allocatable::ary0(:,:),y0(:)
    real(8),allocatable::ary1(:,:),ary2(:,:),x1(:),x2(:),x3(:),x4(:)
    real(8)::fact0=1.0D-30,det
    integer::ib

    nt=20
    allocate(ary0(1:nt,1:nt),ary1(1:nt,1:nt))
    allocate(y0(1:nt),x1(1:nt),x2(1:nt),x3(1:nt),x4(1:nt))
    n=11

!    ary0(1,1)= 2.0D0; ary0(1,2)= 3.0D0; ary0(1,3)= 1.0D0; y0(1)= 4.0D0
!    ary0(2,1)= 4.0D0; ary0(2,2)= 1.0D0; ary0(2,3)=-3.0D0; y0(2)=-2.0D0
!    ary0(3,1)=-1.0D0; ary0(3,2)= 2.0D0; ary0(3,3)= 2.0D0; y0(3)= 2.0D0
    !連立一次方程式[a]{x}={y}の解{x}を求める．ary(i,j)=[a|y]を入力する
    !入力出力結果
    ![a]={ 2, 3, 1} {y}={ 4}{x}={ 2}
    !    { 4, 1,-3}     {-2}    {-1}
    !    {-1, 2, 2}     { 2]    { 3}


    !ヒルベルト行列と解が{1}となる係数ベクトル算出
    do i=1,n
        do j=1,n
            ary0(i,j)=1.0D0/(dble(i)+dble(j)-1.0D0)
        end do
    end do
    !{y}=[ary0]{x}
    do i=1,n
        x1(i)=1.0D0
    end do
    do i=1,n
        y0(i)=0.0D0
        do j=1,n+1
            y0(i)=y0(i)+ary0(i,j)*x1(j)
        end do
    end do


    !連立一次方程式（河西さん）
    do i=1,n
        do j=1,n
            ary1(i,j)=ary0(i,j)
        end do
        ary1(i,n+1)=y0(i)
    end do
    call STEGJ(n,ary1,nt)
    do i=1,n
        x4(i)=ary1(i,n+1)
    end do


    !連立一次方程式（山本さん）
    do i=1,n
        do j=1,n
            ary1(i,j)=ary0(i,j)
        end do
        x1(i)=y0(i)
    end do
    call GJMAT(n,ary1,x1,nt)


    !逆行列(FEMで使用)
    do i=1,n
        do j=1,n
            ary1(i,j)=ary0(i,j)
        end do
    end do
    call MATINV(n,ary1,nt)
    !{x}=[ary0]^(-1){y}
    do i=1,n
        x2(i)=0.0D0
        do j=1,n
            x2(i)=x2(i)+ary1(i,j)*y0(j)
        end do
    end do


    !バンドマトリックスCholesky法(FEMで使用)
    do i=1,n
        do j=1,n
            ary1(i,j)=ary0(i,j)
        end do
        x3(i)=y0(i)
    end do
    ib=IBAND(n,ary1,fact0,nt)
    call BAND(n,ib,ary1,nt)
    call CHBAND10(n,ib,ary1,nt)!三角行列
    call CHBAND11(n,ib,ary1,x3,nt)!前進消去・後退代入


    write(6,*) 'n=',n
    do i=1,n
        write(6,'(4f12.6)') x1(i),x2(i),x3(i),x4(i)
    end do

stop

contains


    subroutine STEGJ(n,ary,nt)
        !Gauss-Jordan法による連立一次方程式の解法
        !河西朝雄著「C 言語によるはじめてのアルゴリズム入門」
        integer,intent(in)::n,nt
        real(8),intent(inout)::ary(1:nt,1:nt)

        integer::i,j,k,s
        real(8)::p,d,max,dumy

        do k=1,n
            max=0.0D0
            s=k
            do j=k,n
                if(abs(ary(j,k))>max)then
                    max=abs(ary(j,k))
                    s=j
                end if
            end do
            do j=1,n+1
                dumy=ary(k,j)
                ary(k,j)=ary(s,j)
                ary(s,j)=dumy
            end do
            p=ary(k,k)
            do j=k,n+1
                ary(k,j)=ary(k,j)/p
            end do
            do i=1,n
                if(i/=k)then
                    d=ary(i,k)
                    do j=k,n+1
                        ary(i,j)=ary(i,j)-d*ary(k,j)
                    end do
                end if
            end do
        end do
    end subroutine STEGJ


    subroutine GJMAT(n,ary,b,nt)
        !ガウス-ジョルダン法(逆行列＋解)
        !http://akita-nct.jp/yamamoto/lecture/2006/5E/Linear_eauations/gaussj_html/
        !著者: 山本昌志
        integer,intent(in)::n,nt
        real(8),intent(inout)::ary(1:nt,1:nt),b(1:nt)

        integer::ipv,i,j
        real(8)::inv_pivot,temp
        real(8)::big
        integer::pivot_row,row(1:nt)

        do ipv=1,n
            !----最大値探索----------------------------
            big=0.0D0
            do i=ipv,n
                if(abs(ary(i,ipv))>big)then
                    big=abs(ary(i,ipv))
                    pivot_row=i
                end if
            end do
            if(big==0.0D0)return !異常終了
            row(ipv)=pivot_row

            !----行の入れ替え--------------------------
            if(ipv/=pivot_row)then
                do i=1,n
                    temp=ary(ipv,i)
                    ary(ipv,i)=ary(pivot_row,i)
                    ary(pivot_row,i)=temp
                end do
                temp=b(ipv)
                b(ipv)=b(pivot_row)
                b(pivot_row)=temp
            end if

            !----対角成分=1(ピボット行の処理)----------
            inv_pivot=1.0D0/ary(ipv,ipv)
            ary(ipv,ipv)=1.0D0
            do j=1,n
                ary(ipv,j)=ary(ipv,j)*inv_pivot
            end do
            b(ipv)=b(ipv)*inv_pivot

            !----ピボット列=0(ピボット行以外の処理)----
            do i=1,n
                if(i/=ipv)then
                    temp=ary(i,ipv)
                    ary(i,ipv)=0.0D0
                    do j=1,n
                        ary(i,j)=ary(i,j)-temp*ary(ipv,j)
                    end do
                    b(i)=b(i)-temp*b(ipv)
                end if
            end do
        end do

        !----列の入れ替え(逆行列)--------------------------
        do j=n,1,-1
            if(j/=row(j))then
                do i=1,n
                    temp=ary(i,j)
                    ary(i,j)=ary(i,row(j))
                    ary(i,row(j))=temp
                end do
            end if
        end do
    end subroutine GJMAT


    subroutine MATINV(n,ary,nt)
        !*************************************
        !逆行列計算
        !山田嘉昭著：マトリックス法材料力学
        !*************************************
        integer,intent(in)::n,nt
        real(8),intent(inout)::ary(1:nt,1:nt)

        integer::i,j,k,l,l1,irow,icol,jrow,jcol
        real(8)::aamax,s,t
        real(8)::factor=1.0D-30
        integer::ipivot(n)
        real(8)::pivot(n)
        integer::nindex(n,2)

        do j=1,n
            ipivot(j)=0
        end do
        do i=1,n
            aamax=0.0D0
            do j=1,n
                if(ipivot(j)-1/=0)then
                    do k=1,n
                        if(0<ipivot(k)-1)return
                        if(ipivot(k)-1<0)then
                            if(abs(aamax)-abs(ary(j,k))<0)then
                                irow=j
                                icol=k
                                aamax=ary(j,k)
                            end if
                        end if
                    end do
                end if
            end do
            ipivot(icol)=ipivot(icol)+1
            if(irow-icol/=0)then
                do l=1,n
                    s=ary(irow,l)
                    ary(irow,l)=ary(icol,l)
                    ary(icol,l)=s
                end do
            end if
            nindex(i,1)=irow
            nindex(i,2)=icol
            pivot(i)=ary(icol,icol)
            if(abs(pivot(i))<factor)then
                write(6,*) '#Singular'
                stop
            end if
            ary(icol,icol)=1.0D0
            do l=1,n
                ary(icol,l)=ary(icol,l)/pivot(i)
            end do
            do l1=1,n
                if(l1-icol/=0)then
                    t=ary(l1,icol)
                    ary(l1,icol)=0.0D0
                    do l=1,n
                        ary(l1,l)=ary(l1,l)-ary(icol,l)*t
                    end do
                end if
            end do
        end do
        do i=1,n
            l=n+1-i
            if(nindex(l,1)-nindex(l,2)/=0)then
                jrow=nindex(l,1)
                jcol=nindex(l,2)
                do k=1,n
                    s=ary(k,jrow)
                    ary(k,jrow)=ary(k,jcol)
                    ary(k,jcol)=s
                end do
            end if
        end do
    end subroutine MATINV


    subroutine CHBAND10(n1,m1,array,nt)
        integer,intent(in)::n1,m1,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        integer :: i,j,k,mm,mn
        real(8) :: z

        array(1,1)=sqrt(array(1,1))
        do j=2,m1
            array(1,j)=array(1,j)/array(1,1)
        end do
        do i=2,n1
            if(n1-i+1>m1)then
                mm=m1
            else
                mm=n1-i+1
            end if
            do j=1,mm
                z=array(i,j)
                if(j/=m1)then
                    if(m1-j+1<i)then
                        mn=m1-j+1
                    else
                        mn=i
                    end if
                    do k=2,mn
                        z=z-array(i-k+1,k)*array(i-k+1,j+k-1)
                    end do
                    if(j>1)then
                        array(i,j)=z/array(i,1)
                    else
                        array(i,1)=sqrt(z)
                    end if
                else
                    array(i,j)=z/array(i,1)
                end if
            end do
        end do
    end subroutine CHBAND10


    subroutine CHBAND11(n1,m1,array,vector,nt)
        integer,intent(in)::n1,m1,nt
        real(8),intent(in)::array(1:nt,1:nt)
        real(8),intent(out)::vector(nt)
        integer :: i,j,mm
        real(8) :: z

        vector(1)=vector(1)/array(1,1)
        do i=2,n1
            z=vector(i)
            if(i>m1)then
                mm=m1
            else
                mm=i
            end if
            do j=2,mm
                z=z-array(i-j+1,j)*vector(i-j+1)
            end do
            vector(i)=z/array(i,1)
        end do

        vector(n1)=vector(n1)/array(n1,1)
        do i=n1-1,1,-1
            z=vector(i)
            do j=2,m1
                if(i+j-1>n1)exit
                z=z-array(i,j)*vector(i+j-1)
            end do
            vector(i)=z/array(i,1)
        end do
    end subroutine CHBAND11


    integer function IBAND(mm,array,fact0,nt)
        integer,intent(in)::mm,nt
        real(8),intent(in)::array(1:nt,1:nt),fact0
        integer :: i,j
        !-----------------------------------------------------
        !バンド幅計算
        !-----------------------------------------------------
        IBAND=1
        do i=1,mm-1
            do j=mm,i+1,-1
                if(fact0<=abs(array(i,j)))exit
            end do
            if(IBAND<=j-i+1)IBAND=j-i+1
        end do
        if(mm<IBAND)IBAND=mm
    end function IBAND


    subroutine BAND(mm,ib,array,nt)
        integer,intent(in)::mm,ib,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        integer :: i,j,k,mn
        real(8) :: work(1:nt)
        !-----------------------------------------------------
        !フルマトリックスのバンド化記憶
        !-----------------------------------------------------
        do i=2,mm
            if(i+ib-1<mm)then
                mn=ib
            else
                mn=mm-i+1
            end if
            do j=1,ib
                work(j)=0.0D0
            end do
            do k=1,mn
                work(k)=array(i,i+k-1)
            end do
            do j=1,ib
                array(i,j)=work(j)
            end do
        end do
    end subroutine BAND


end program f90_GJ