Demo entry 4076518

test

   

Submitted by anonymous on Mar 16, 2016 at 22:07
Language: Fortran. Code size: 14.5 kB.

module zfad_solve
    use libraries
    use zgap_paras
    use zgap_utils
    use zfad_utils
    use zfad_paras
    implicit none

private
    type(quark_loader):: qrk
    type(zero_fadq_para),target:: para
    type(chey_cls):: cheys2,cheyxs,cheyzp,cheyzq,cheyz0

    complex(8),target:: tP(4)
    real(8),allocatable:: s2(:,:),xs(:,:),pz(:,:),qz(:,:),x0(:,:)
    real(8),allocatable:: qk2(:,:),qkz(:,:),qk0(:,:),qkt(:,:)
    type (fadv_mom),allocatable:: moms(:,:,:,:,:)
    complex(8),allocatable,target:: E(:,:)
    complex(8),allocatable:: R(:),FR(:),F3Y(:,:,:,:,:,:,:)
    complex(8),pointer:: GX(:,:,:,:,:,:,:),GY(:,:,:,:,:,:,:,:),GZ(:,:,:,:,:,:,:,:)
    complex(8),pointer:: F3X(:,:,:,:,:,:,:),MatWav(:,:,:,:,:,:)
    type (MPI_shared_array):: GXptr, GYptr, GZptr, F3Xptr, MatWavptr
    complex(8),allocatable:: MatKer(:,:,:,:,:,:,:,:)
    type (sparse_mat_cx):: H1Mat, H2Mat, GMat, Kmat
    integer,allocatable:: H1ja(:),H1ia(:),H2ja(:),H2ia(:),Gja(:),Gia(:),Kja(:),Kia(:)

    real(8):: totaltimer1,totaltimer2,timer1,timer2

    integer:: jdx1,jdx2,jtasks,idx1,idx2,itasks,lflag

    integer:: Z,I_s2,I_xs,I_pz,I_qz,I_q0,I_qk2,I_qkz,I_qk0,I_qkt,I_Nf,I_Nl

type zfad_solve_cls
    type(zero_fadq_para),pointer:: para
    integer:: idxP = 0
    integer:: maxP = 0
    complex(8),pointer:: tP(:)
    complex(8),pointer:: E(:,:)
contains
    procedure:: info => zfad_solve_info
    procedure:: initiate => calculate_in
    procedure:: finalize => calculate_out
    procedure:: solve => zfad_solve_main
end type

public zfad_solve_cls

CONTAINS

subroutine zfad_solve_info(this)
    class(zfad_solve_cls):: this

    if(jrank /= 0) return !only allow main process output

    print*,CG
    print*,"****************************************************************"
    print*,"*                                                              *"
    print*,"*      Solve Faddeev equation at zero T (complex plane)        *"
    print*,"*                                                              *"
    print*,"****************************************************************"
    print*,CR
    write(*,"('kernel : ',I4,' * ',f8.5,' GB')") &
        & jsize,dble(product(shape(MatKer)))*16d0/10d0**9
    write(*,"('baryon,[P2,ImP] : ',a8,'  [',f8.5,',',f8.5,']')") &
        & para%baryon,real(tP(4)**2),imag(tP(4))
    write(*,"('gln,rnm,D,omega,mm : ',2a4,3f10.5)") &
        & para%pgap%gln,para%pgap%rnm,&
        & para%pgap%D,para%pgap%omega,para%pgap%mm
    print*,CX
end subroutine

subroutine zfad_solve_main(this)
    class(zfad_solve_cls):: this

    !#call calculate_in(this)

    if( jrank == 0 ) call cpu_time(totaltimer1)

    call calculate_matrix()     !#calculate kernel matrix

    if( jrank == 0 ) then
        call eigiter(E, calculate_residue, 1, 50)
        lflag=0; call calculate_residue(Z,R,FR)
    else
        do while( lflag /= 0)
            call calculate_residue(Z,R,FR)
        end do
    end if

    if( jrank == 0 ) then
		call cpu_time(totaltimer2)
		print*,CR,'timmer(s): ',totaltimer2-totaltimer1,CX
	end if

    !#call calculate_out(this)
end subroutine

subroutine calculate_in(this)
    class(zfad_solve_cls):: this
    integer:: idxP

    para = zfad_load_para()

    jtasks = para%Cs2 !#size of jobs
    call MPI_distribute(jdx1,jdx2,jtasks,lflag)

    itasks = para%Cs2 !#size of shared jobs
    call MPI_distribute_shared(idx1,idx2,itasks)

    if (jrank == 0 .and. this%idxP == 0) then
        !#call append (trim(para%ofile)//'.txt', para%baryon)
    end if

    this%maxP = ceiling((para%maxP - para%minP)/para%deltaP)
    if (para%tyP == 'RE') then
        tP = (/Zero,Zero,Zero, sqrt(One*( para%maxP - this%idxP * para%deltaP ))/)
    else if (para%tyP == 'IM') then
        tP = (/Zero,Zero,Zero, Im*( para%maxP - this%idxP * para%deltaP )/)
    end if
    this%idxP = this%idxP + 1

    Z = Nf * Nl * para%Cs2 * para%Cxs * para%Cxz * para%Cxz * para%Cx0 !#set the total dimension

    call Mat_H1_IDX(H1Mat%ja, H1Mat%ia); call Mat_H2_IDX(H2Mat%ja, H2Mat%ia)
    call Mat_G_IDX(GMat%ja, GMat%ia); call Mat_K_IDX(KMat%ja, KMat%ia)
    allocate(H1Mat%a(size(H1Mat%ja))); allocate(H2Mat%a(size(H2Mat%ja)))
    allocate(GMat%a(size(GMat%ja))); allocate(KMat%a(size(KMat%ja)))

    allocate(s2(para%Cs2,2),xs(para%Cxs,2),pz(para%Cxz,2),qz(para%Cxz,2),x0(para%Cx0,2))
    allocate(moms(para%Cs2,para%Cxs,para%Cxz,para%Cxz,para%Cx0))
    allocate(R(Z),FR(Z),E(Z+2,1))

    allocate(F3Y(Nl,Nf,para%Cxz,para%Cs2,para%Cxs,para%Cxz,para%Cx0))  !# F(p,q) = F(p2,zp,q2,zq,z0)
    allocate(MatKer(size(KMat%a),jdx1:jdx2,para%Cxs,para%Cxz,para%Cx0,para%Mx2,para%Mxz,para%Mx0)) !# K(q,qk)

    call MPI_shared_alloc (F3Xptr, (/Nl,Nf,para%Cxz,para%Cs2,para%Cxs,para%Cxz,para%Cx0/), F3X)
    call MPI_shared_alloc (GXptr, (/Nl,Nf,para%Cxz,para%Cs2,para%Cxs,para%Cxz,para%Cx0/), GX)
    call MPI_shared_alloc (GZptr, (/Nl,Nf,para%Cxz,para%Cs2,para%Cxs,para%Mx2,para%Cxz,para%Cx0/), GZ)
    call MPI_shared_alloc (GYptr, (/Nl,Nf,para%Cxz,para%Cs2,para%Cxs,para%Mx2,para%Mxz,para%Mx0/), GY)
    call MPI_shared_alloc (MatWavptr, (/size(GMat%a),para%Cxz,para%Cs2,para%Cxs,para%Cxz,para%Cx0/), MatWav)

    call cheys2%init(s2, para%I_F, para%U_V, 'E', 'a')   !!!# s2 = 4q2 + 3p2
    call cheyxs%init(xs, -1d0, 1d0, 'D', 'a')   !!!# x = (4q2-3p2)/s2
    call cheyzp%init(pz, -1d0, 1d0, 'D', 'a')
    call cheyzq%init(qz, -1d0, 1d0, 'D', 'a')
    call cheyz0%init(x0, -1d0, 1d0, 'D', 'a')

    qk2  = intpoint('GS',para%Mx2,para%I_F,para%U_V,mapexp,radial4R)
    qkz  = intpoint('C2',para%Mxz)
    qk0  = intpoint('GS',para%Mx0)
    qkt  = intpoint('GS',para%Mxt,0d0,2d0*pi,maplin)

    call qrk%load(para)

    this%para => para
    this%tP => tP
    this%E => E

end subroutine

subroutine calculate_out(this)
    class(zfad_solve_cls):: this
    complex(8):: row(5)

    if(jrank == 0)then !#only allow main process output

    print*,CG,E(Z+1,:),CX

    call append (trim(para%ofile)//'.txt',(/real(E(Z+1,1)),imag(E(Z+1,1))/),"","Re|Im")

    end if

    deallocate( s2,xs,pz,qz,x0,qk2,qkz,qk0,qkt,E,R,FR,F3Y,MatKer )
    call MPI_shared_free(GXptr)
    call MPI_shared_free(GZptr)
    call MPI_shared_free(GYptr)
    call MPI_shared_free(F3Xptr)
    call MPI_shared_free(MatWavptr)
end subroutine

subroutine calculate_matrix()
    complex(8):: FF(4)
    real(8):: p(4),q(4),qk(4),psqu,qsqu
    type (fadv_mom):: mom

    do I_s2=1,para%Cs2
        do I_xs=1,para%Cxs
            do I_pz=1,para%Cxz
                do I_qz=1,para%Cxz
                    do I_q0=1,para%Cx0
                        mom%s2 = s2(I_s2,1)
                        mom%xs = xs(I_xs,1)
                        mom%zp = pz(I_pz,1)
                        mom%zq = qz(I_qz,1)
                        mom%z0 = x0(I_q0,1)
                        call mom_permutation(mom)
                        moms(I_s2,I_xs,I_pz,I_qz,I_q0) = mom
                    end do
                end do
            end do
        end do
    end do

    do I_q0=1,para%Cx0
        do I_qz=1,para%Cxz
            do I_xs=1,para%Cxs
                do I_s2=idx1,idx2 !#1,para%Cs2
                    do I_pz=1,para%Cxz
                        mom = moms(I_s2,I_xs,I_pz,I_qz,I_q0)
                        !mom%p = 0d0 !*************
                        call qrk%calc (tP, mom%p, mom%q, FF, para)
                        !FF(1) = 1d0; FF(2) = 0.5d0; FF(3) = 1d0; FF(4) = 0.5d0
                        call Proj_G(tP, mom%p, mom%q, FF, MatWav(:,I_pz,I_s2,I_xs,I_qz,I_q0))
                    end do
                end do
            end do
        end do
    end do
    call MPI_shared_sync (MatWavptr)

    MatKer = 0d0
    do I_xs=1,para%Cxs
        do I_q0=1,para%Cx0
            do I_qz=1,para%Cxz
                do I_s2=jdx1,jdx2 !#1,para%Cs2, cross-node parallel
                    qsqu = s2(I_s2,1) * (1 + xs(I_xs,1))/2d0
                    q = (/0d0,sqrt((1d0-x0(I_q0,1)**2)*(1d0-qz(I_qz,1)**2)),&
                      & x0(I_q0,1)*sqrt(1d0-qz(I_qz,1)**2),qz(I_qz,1)/) * sqrt(qsqu)
                    do I_qk0=1,para%Mx0
                        do I_qkz=1,para%Mxz
                            do I_qk2=1,para%Mx2
                                qk = (/0d0,sqrt((1d0-qk0(I_qk0,1)**2)*(1d0-qkz(I_qkz,1)**2)),&
                                   & qk0(I_qk0,1)*sqrt(1d0-qkz(I_qkz,1)**2),qkz(I_qkz,1)/) * sqrt(qk2(I_qk2,1))
                                do I_qkt=1,para%Mxt
                                    call Kerl_K(q, qk, qkt(I_qkt,1), para, KMat%a)
                                    MatKer(:,I_s2,I_xs,I_qz,I_q0,I_qk2,I_qkz,I_qk0) &
                                        & = MatKer(:,I_s2,I_xs,I_qz,I_q0,I_qk2,I_qkz,I_qk0) + KMat%a &
                                        & * (qk2(I_qk2,2) * qkz(I_qkz,2) * qk0(I_qk0,2) * qkt(I_qkt,2))
                                end do
                            end do
                        end do
                    end do
                end do
            end do
        end do
        if(jrank==0) write (6,"(A,I3,'%')",advance="no") char(13),int(I_xs*100d0/(para%Cxs)); call flush(6)
    end do
end subroutine

subroutine calculate_residue(Z,X,F)
    integer:: Z
    complex(8):: X(Z),F_local(Z),F(Z)

    if( jrank == 0 ) call cpu_time(timer1)

    call MPI_broadcast(lflag)
    if( lflag /= 0 )then
        call MPI_broadcast(X)

        F_local=0.d0; call fad_equation(X,F_local)

        call MPI_reducesum(F_local,F)
    end if

    if( jrank == 0 ) then
		call cpu_time(timer2)
        write (6,"(A,'[',f4.1,'s]')",advance="no") char(13),timer2-timer1; call flush(6)
	end if

end subroutine

subroutine fad_equation(X,F)
    complex(8):: X(:),F(:)
    real(8):: ts2,txs,psqu,qsqu
    complex(8):: amp1(Nl,Nf),amp2(Nl,Nf),wav(Nl,Nf)

    F3X = reshape (X, shape(F3X)) !I_Nl,I_Nf,I_pz,I_s2,I_xs,I_qz,I_q0

    do I_q0=1,para%Cx0
        do I_qz=1,para%Cxz
            do I_xs=1,para%Cxs
                do I_s2=idx1,idx2 !#1,para%Cs2, in-node parallel
                    do I_pz=1,para%Cxz
                        call amp_permutation(cheyxs,cheyzp,cheyzq,cheyz0, &
                                & moms(I_s2,I_xs,I_pz,I_qz,I_q0),&
                                & H1Mat,H2Mat,F3X(:,:,:,I_s2,:,:,:),amp1,amp2)
                        wav = F3X(:,:,I_pz,I_s2,I_xs,I_qz,I_q0) !+ amp1 + amp2
                        !if(jrank == 0) print*,"x",amp1(1,1)
                        GMat%a = MatWav(:,I_pz,I_s2,I_xs,I_qz,I_q0)
                        do I_Nf=1,Nf !************
                            wav(:,I_Nf) = matmul2(GMat, wav(:,I_Nf))
                            !wav(1,I_Nf) = GMat%a(1) * wav(1,I_Nf)
                        end do

                        GX(:,:,I_pz,I_s2,I_xs,I_qz,I_q0) = wav
                    end do
                end do
            end do
        end do
    end do
    call MPI_shared_sync (GXptr)

    do I_qk2=1,para%Mx2
        do I_xs=1,para%Cxs
            do I_s2=idx1,idx2 !#1,para%Cs2, in-node parallel
                psqu = s2(I_s2,1) * (1d0 - xs(I_xs,1))*2d0/3d0 !*********
                ts2 = qk2(I_qk2,1) + 3d0/4d0*psqu; txs = (qk2(I_qk2,1) - 3d0/4d0*psqu)/ts2
                if (ts2 > cheys2%ux) ts2 = cheys2%ux; if (ts2 < cheys2%lx) ts2 = cheys2%lx
                call cheys2%tran(ts2); call cheyxs%tran(txs)
                do I_qz=1,para%Cxz
                    do I_q0=1,para%Cx0
                        do I_pz=1,para%Cxz
                            do I_Nf=1,Nf
                                do I_Nl=1,Nl
                                    GZ(I_Nl,I_Nf,I_pz,I_s2,I_xs,I_qk2,I_qz,I_q0) &
                                        & = mulchey(cheys2%cheytran(1,:),cheyxs%cheytran(1,:),&
                                        & GX(I_Nl,I_Nf,I_pz,:,:,I_qz,I_q0))
                                end do
                            end do
                        end do
                    end do
                end do
            end do
        end do
    end do
    call MPI_shared_sync (GZptr)

    do I_qk0=1,para%Mx0
        do I_qkz=1,para%Mxz
            call cheyzq%tran(qkz(I_qkz,1)); call cheyz0%tran(qk0(I_qk0,1))
            do I_qk2=1,para%Mx2
                do I_xs=1,para%Cxs
                    do I_s2=idx1,idx2 !#1,para%Cs2, in-node parallel
                        do I_pz=1,para%Cxz
                            do I_Nf=1,Nf
                                do I_Nl=1,Nl
                                    GY(I_Nl,I_Nf,I_pz,I_s2,I_xs,I_qk2,I_qkz,I_qk0) &
                                        & = mulchey(cheyzq%cheytran(1,:),cheyz0%cheytran(1,:),&
                                        & GZ(I_Nl,I_Nf,I_pz,I_s2,I_xs,I_qk2,:,:))
                                end do
                            end do
                        end do
                    end do
                end do
            end do
        end do
    end do
    call MPI_shared_sync (GYptr)

    F3Y = 0d0
    do I_qk0=1,para%Mx0
        do I_qkz=1,para%Mxz
            do I_qk2=1,para%Mx2
                do I_q0=1,para%Cx0
                    do I_qz=1,para%Cxz
                        do I_xs=1,para%Cxs
                            do I_s2=jdx1,jdx2 !#1,para%Cs2, cross-node parallel
                                KMat%a = MatKer(:,I_s2,I_xs,I_qz,I_q0,I_qk2,I_qkz,I_qk0)
                                do I_pz=1,para%Cxz
                                    do I_Nf=1,Nf !*********
                                        F3Y(:,I_Nf,I_pz,I_s2,I_xs,I_qz,I_q0) &
                                            & = F3Y(:,I_Nf,I_pz,I_s2,I_xs,I_qz,I_q0) &
                                            & + matmul2(KMat, GY(:,I_Nf,I_pz,I_s2,I_xs,I_qk2,I_qkz,I_qk0) )
                                        !F3Y(1,I_Nf,I_pz,I_s2,I_xs,I_qz,I_q0) &
                                        !    & = F3Y(1,I_Nf,I_pz,I_s2,I_xs,I_qz,I_q0) &
                                        !    & + KMat%a(1)*GY(1,I_Nf,I_pz,I_s2,I_xs,I_qk2,I_qkz,I_qk0)
                                    end do
                                end do
                            end do
                        end do
                    end do
                end do
            end do
        end do
    end do

    F = reshape (F3Y, shape(F))
end subroutine

end module

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).