program main
    implicit none
    integer m1,m2,m3
    call main1()
    stop
    end program 
!=======================================================    
    subroutine main1()
    ! 读取不明行数的文件
    implicit none
    integer count,stat
    real h1,h2,h3,h4,h5,h22
    h1 = 0
    h2 = 0
    h3 = 0
    h22 = 0
    open(1,file = "test.dat")
    open(2,file = "test-format.dat")
    count = 0
    do while (.true.)
        count = count + 1
        h22 = h1
        read(1,*,iostat = STAT)h1,h2,h3,h4,h5
        if(h22.ne.h1)write(2,*)""  ! 在这里加空行是为了gnuplot绘制密度图
        write(2,999)h1,h2,h3,h4,h5   ! 数据格式化
        if(stat .ne. 0) exit ! 当这个参数不为零的时候,证明读取到文件结尾
    end do
    ! write(*,*)h1,h2,h3
    ! write(*,*)count
    close(1)
    close(2)
999 format(10f11.6)
    return
    end subroutine main1

    module pub
    implicit none
    integer N,iternum,hn
    real err,eta,dk,domg
    parameter( hn = 2, N = hn,iternum = 200,err = 1e-16,eta = 0.01,dk = 0.01,domg = dk)
    real,parameter::pi = 3.1415926535
    complex,parameter::im = (0.,1.0)  
    complex ones(N,N),GLL(N,N),GRR(N,N),GB(N,N)
    complex H00(N,N)  ! diagonal elementery
    complex H01(N,N)  ! off-diag elementery
    complex g1(hn,hn),g2(hn,hn),g3(hn,hn)
    !---------------------------------------------
    real m0,mu   
    real tx,ty
    real ax,ay
    end module pub
!==================================================================
    program main
    use pub
    !======parameter value setting =====
    m0 = 1.5
    tx = 1.0
    ty = 1.0
    ax = 1.0
    ay = 1.0
    call surfstat()
    stop
    end program main
!============================================================================================
    subroutine surfstat()
    ! surfstat calculates surface state using green's function method---J.Phys.F.Met.Phys.15(1985)851-858     
    ! 利用已经求得的格林函数来计算对应的态密度
    use pub
    implicit none
    real kx,omg,re1,re2,re3
    real t_start,t_end
    integer i1
    open(20,file="dos.dat")
    call cpu_time(t_start)
    !------------------------------------------
    do kx = -pi,pi,dk
        call matset(kx)
        do omg = -3,3,domg
            call itera(omg,kx)
            re1 = log(abs(sum(aimag(GLL))))
            re2 = log(abs(sum(aimag(GRR))))
            re3 = log(abs(sum(aimag(GB))))
            write(20,999)kx/pi,omg,re1,re2,re3
        end do
        write(20,*)" "
    end do
    !------------------------------------------
    call cpu_time(t_end)
    close(20)
    write(*,*)"Timing const is: ",t_end - t_start
999 format(30f16.12)
    return
    end subroutine surfstat
!=================================================================
    subroutine itera(omega,ky)
    use pub
    real omega,real_temp,ky
    integer iter
    complex omegac
    complex g0dem(N,N), g0(N,N) ! Green's Function
    complex epsiloni(N,N),epsilons(N,N),epsilons_t(N,N),alphai(N,N),betai(N,N)   ! 迭代过程变量
    complex GLLdem(N,N),GRRdem(N,N),GBdem(N,N),mat1(N,N),mat2(N,N)
    !----------------------------
    call matset(ky)
    epsiloni = H00
    epsilons = H00
    epsilons_t = H00
    alphai = H01
    betai  = conjg(transpose(H01))
    omegac = omega + eta*im
    !----------------------------------
    do iter = 1, iternum

        g0dem = omegac*ones - epsiloni  !  Green's Function
        call inv(g0dem, g0)

        mat1 = matmul(alphai, g0)
        
        mat2 = matmul(betai, g0)

        g0 = matmul(mat1,betai)

        epsiloni = epsiloni + g0

        epsilons = epsilons + g0

        g0 = matmul(mat2,alphai)

        epsiloni = epsiloni + g0

        epsilons_t = epsilons_t + g0

        g0 = matmul(mat1, alphai)
        alphai = g0

        g0 = matmul(mat2, betai)
        betai = g0
          
        real_temp = sum(abs(alphai))   
        if (real_temp .le. err)then
            exit
        end if

     end do ! end of iteration

     ! calculate surface green's function
     GLLdem = omegac*ones- epsilons
     call inv(GLLdem, GLL)
    !  GLL = epsilons


     GRRdem = omegac*ones- epsilons_t
     call inv(GRRdem, GRR)
    !  GRR = epsilons_t

     GBdem = omegac*ones- epsiloni
     call inv(GBdem, GB)
    !  GB = epsiloni
    end subroutine itera
!================================================================
    subroutine matset(ky)
    ! 构建哈密顿量
    use pub
    real ky
    integer m,l
    call Pauli()
    do m = 1,hn
        do l = 1,hn
            H00(m,l) = (m0-ty*cos(ky))*g1(m,l) + ay*sin(ky)*g3(m,l) 

            H01(m,l) = (-tx*g1(m,l) - im*ax*g2(m,l))/2
        end do
    end do
    !----------------------
    !  初始化单位矩阵
    do ix = 1,N
        ones(ix,ix) = 1.0
    end do
    return
    end subroutine matset
!=======================矩阵求逆====================================
    subroutine inv(matin,matout)
    use pub
    complex,intent(in) :: matin(N,N)  ! N is dimension of matrix which can be readed from pub model (use pub)
    complex:: matout(size(matin,1),size(matin,2))
    real:: work2(size(matin,1))            ! work2 array for LAPACK
    integer::ipiv(size(matin,1))     ! pivot indices
    integer info  ! inv solution information
    ! Store matin in matout to prevent it from being overwritten by LAPACK
    matout = matin
    ! SGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    call CGETRF(N,N,matout,N,ipiv,info)
    if (info.ne.0) stop 'Matrix is numerically singular!'
    ! SGETRI computes the inverse of a matrix using the LU factorization
    ! computed by SGETRF.
    call CGETRI(N,matout,N,ipiv,work2,N,info)
    if (info.ne.0) stop 'Matrix inversion failed!'
    return
    end subroutine inv
!====================================================
    subroutine Pauli()
    use pub
    !----------
    g1(1,1) = 1
    g1(2,2) = -1
    !----------
    g2(1,2) = 1
    g2(2,1) = 1
    !----------
    g3(1,2) = -im
    g3(2,1) = im
    return
    end subroutine Pauli